TY - JOUR
AU1 - Roberts,, Lewis
AU2 - Griffith,, Thom
AU3 - Champneys,, Alan
AU4 - Piano,, Martina
AU5 - Kiely,, Janice
AU6 - Luxton,, Richard
AB - Abstract A mathematical model is developed to describe the action of a novel form of fluidic biosensor that uses paramagnetic particles (PMPs) that have been pre-coated with target-specific antibodies. In an initial phase the particles are introduced to a sample solution containing the target which then binds to the particles via antigen–antibody reactions. During the test phase a magnet is used to draw the PMPs to the sensor surface which is similarly coated with specific antibodies. During this process, cross-links are formed by the antigens thereby binding the PMPs to the sensor surface. After the magnetic field is removed, a voltage change across an inductor below the sensor surface is recorded, which is deemed to depend on the number of magnetic particles that have been bound to the sensor surface. The fundamental question addressed is to explain the range of experimentally observed dose–response curves, and how this depends on the various parameters of the problem. In particular, observations have shown both rising and falling dose–response curves, as well as ‘hooked’ dose–response curves possessing local maxima. Initially a particle-dynamics computational model is produced to determine the time scales of the key processes involved, but is shown to be unable to produce differently shaped dose–response curves. The computational model suggests spatio-temporal effects are unimportant, therefore a homogenized rate-equation model is developed for each of the key phases of the immunoassay process. Binding rates are shown to depend on various geometric factors related to the diameter of the PMPs and the size of the sensor surface. The dose–response is shown to depend crucially on various saturation effects during each phase, and conditions can be derived, in some cases analytically, for each of the three qualitatively different curve types. Furthermore, non-dimensionalization reveals 5 key dimensionless parameters and the dependence of these curve shapes on each is revealed. The results point to future quantitative approaches to sensor design and calibration. 1. Introduction Immunoassay techniques exploiting the high specificity of antigen–antibody pair binding are regularly used as biosensors to determine the amount of target protein present in a sample. In a magnetic immunoassay, reactive magnetic particles (particles that have been coated with antibodies) provide the means for separating the target antigen from the sample as well as the signal by which the target antigen is quantified. The targeted protein may represent an important biomarker in the early detection of disease, therefore the potential for fast separation and quantification of target antigen makes the development of magnetic immunoassays important in the field of biomedical diagnostics (Koh & Josephson, 2009). Many different types of biosensor technology have been developed for the detection and measurement of magnetic particles (see Kokalj et al., 2015; Jamshaid et al., 2016 for recent reviews). In immunoassay design, to ensure the assay is sensitive to the required range of target antigen concentrations, there are several experimental parameters that need to be optimized. Such parameters include the number of magnetic particles per assay, the capture-antibody density on the magnetic particle surface, and, if a sensor surface is employed, the density of capture antibodies on the sensor surface. The key is to ensure that the output signal has maximal sensitivity in the range of target concentrations that is clinically important (see e.g. Xu et al., 2006; Barnett et al., 2014). Also important is to eliminate signal interferences such as the high-dose hook effect where large analyte concentrations block the immobilization of particles at the sensor surface. The hook effect is particularly problematic for the practical application of immunoassay technology in biomedical diagnostics since it can cause false negatives and hence potential misdiagnoses. However, the hook effect can be avoided through careful design and optimization of the experimental parameters in the assay (Cole et al., 2001). As is common in commercial assay development, the typical approach used to optimize experimental parameters is painstaking and involves many trial-and-error repeats. However, unlike high-throughput techniques, the design of an assay for a new biomarker, is typically carried out by hand in a laboratory, and therefore each repeat is both time consuming and costly. So far, mathematical modelling techniques have been used only sparingly within the magneto-immunoassay field. For example, the coating of magnetic particles used in such assays have recently been modelled using random sequential adsorption theory (Mackey et al., 2016). Also, simple Langmuir kinetics have been applied in order to understand the process of antigen capture (Saha et al., 2014). At the other scale of complexity, computational fluid dynamics (CFD) models have been used to model the motion of the magnetic particles in an assay when under the influence of external magnetic fields (see, e.g. Munir et al., 2011). The purpose of this article is to describe the construction of an end-to-end mathematical model of a magnetic immunoassay. Specifically, we shall describe a novel assay, developed at the University of the West of England (Richardson et al., 2001) and further developed in Sharif et al. (2013) and Barnett et al. (2014), that uses inductive sensors to determine the amount of magnetic particles that have been cross-linked to the sensor surface through antibodies attached to both the particles and the sensor surface. We seek a model that is simple enough to provide rapid insight into how the output signal depends on all the design parameters, yet is sufficiently complex that all the key processes are captured. We shall proceed in two steps. First a stochastic, computational model of the dynamics of the magnetic particles will be constructed. This takes the form of a particle dynamics model; essentially a continuous-time-and-space agent-based model. The purpose of that model is to test the hypothesis that a spatially homogenized approximation is likely to describe the physics with sufficient accuracy. Then we proceed to derive rate equations in the form of systems of ordinary differential equations (ODEs). Analysis of these ODEs then leads to predictions for the signal dependency on all design parameters. Our overall objective is to demonstrate the power of the modelling approach to produce a quantitative tool that can rapidly circumvent the current time-consuming approach to immunoassay design and calibration. The rest of the article is outlined as follows. Section 2 describes the particular magnetic immunoassay we seek to model, and presents typical laboratory output. Section 3 then describes and presents results from the computational model. Section 4 develops the deterministic model and performs non-dimensionalization to identify several dimensionless parameters. Section 5 presents results from the deterministic model in which under certain simplifying assumptions, an asymptotic analysis leads us to an analytical approximation of the signal output. These results are then compared with numerical simulations of the same model. The corresponding dose–response curves are computed and their dependence on each of the dimensionless parameters is revealed. Finally, Section 6 provides concluding remarks, interpreting the results from the model, along with suggestions for future work. 2. The magnetic immunoassay The immunoassay under consideration is depicted schematically in Fig. 1 and, as described in Richardson et al. (2001), Sharif et al. (2013) and Barnett et al. (2014), uses antibody-coated paramagnetic particles (PMPs) and a reaction vessel whose surface is also coated with antibodies. The target antigen must bind to both the PMP and the sensor surface in order for a signal to be recorded. Such a modality requires that the target antigen has multiple epitopes (parts of the molecule to which antibodies may become attached), which is typical of most complex proteins. Fig. 1. Open in new tabDownload slide A schematic depiction of the magnetic immunoassay. Antibody-coated PMPs are added to the sample containing the target antigen |$T$|⁠. The target antigen binds to the antibody-coated PMPs fixing it to the particle surface. When the magnetic field is applied, the PMPs are attracted to the surface. Once the PMPs are close enough to the sensor surface, target antigen molecules bound to the PMPs can cross-link with free antibodies on the sensor surface |$S$|⁠. A blocking agent is used to prevent any target binding to the particle itself or the sensor surface. If enough cross-links are made, the particle is bound to the surface. When the magnetic field is removed, the free particles |$P_F$| drift away but the bound particles |$P$| remain, changing the self-inductance of the sensor coil from its initial state. At any time, a free binding site on the sensor surface can also be occupied by a free target molecule |$T$| inhibiting the binding of active PMPs. This is more likely to occur when the target analyte concentration is high. To accommodate the sensor coil, the magnet stops at a minimum distance |$z_0$| from the sensor surface. The radius of the sensor surface is |$a$|⁠. Fig. 1. Open in new tabDownload slide A schematic depiction of the magnetic immunoassay. Antibody-coated PMPs are added to the sample containing the target antigen |$T$|⁠. The target antigen binds to the antibody-coated PMPs fixing it to the particle surface. When the magnetic field is applied, the PMPs are attracted to the surface. Once the PMPs are close enough to the sensor surface, target antigen molecules bound to the PMPs can cross-link with free antibodies on the sensor surface |$S$|⁠. A blocking agent is used to prevent any target binding to the particle itself or the sensor surface. If enough cross-links are made, the particle is bound to the surface. When the magnetic field is removed, the free particles |$P_F$| drift away but the bound particles |$P$| remain, changing the self-inductance of the sensor coil from its initial state. At any time, a free binding site on the sensor surface can also be occupied by a free target molecule |$T$| inhibiting the binding of active PMPs. This is more likely to occur when the target analyte concentration is high. To accommodate the sensor coil, the magnet stops at a minimum distance |$z_0$| from the sensor surface. The radius of the sensor surface is |$a$|⁠. The procedure for carrying out the assay is as follows. First, antibody-coated PMPs are added to the test solution containing the target antigen, as illustrated schematically in Fig. 1. This solution is then left for an incubation period of approximately 1 min to allow a sufficient amount of target to fix onto the antibodies on the PMP surface. We shall refer to a PMP that has become bound by a target protein as ‘active’. The solution containing the target antigen and antibody-coated PMPs is then introduced to the antibody-coated sensor surface. A PMP becomes bound to the sensor surface if enough cross-links are made between itself and the sensor surface. A cross-link is formed when a target antigen that is bound to a PMP binds to an antibody on the sensor surface via an antibody–antigen reaction. The number of such bound PMPs is measured by a magnetometer, which consists of a inductor coil and a permanent magnet that is an adjustable vertical distance |$z$| from the sensor surface, as indicated in Fig. 1. In order to bring the PMPs in close proximity to the sensor surface, the permanent magnet is moved close to the sensor surface (at what we shall call the on position) to attract the PMPs. When the magnetic field is removed (i.e. the magnet is moved away from the sensor surface to the off position such that the magnetic field at the sensor surface becomes negligible), the unbound particles drift away from the surface via Brownian motion. Once a sufficient number of PMPs have bound to the sensor surface, they exhibit a net magnetization that will change the inductance of the coil and hence the output voltage of the magnetometer. For the particular magnetometer used, the voltage signal is proportional to the number of PMPs bound to the sensor surface (see Section 2.1 below). The underlying measurement principle then is that the number of immobilized PMPs on the surface will depend on the initial concentrations of the target molecule. In essence, the target is acting as the glue that binds the two forms of antibody together. The more glue, the more binding and hence the more magnetic particles that will be captured by the sensor surface once the magnetic field is turned off. In turn, the more bound PMPs the higher the signal in the inductor. Hence, we should expect to see a dose–response curve by plotting the signal from the magnetometer against the concentration of the target molecule in the original solution. It is important though to note that the relationship between the magnetometer signal and the target-antigen concentration need not necessarily be described by a positive monotonic function. If the amount of target introduced to the assay is greater than the number of binding sites or the incubation time is too short to allow the number of binding sites on the PMPs to saturate, then there will be some free target |$T$| in the assay when the solution containing the PMPs is added to the reaction vessel. These target molecules can bind to a free antibody site |$S$| on the sensor surface and thus no cross-link can form there, reducing the chance of immobilizing a PMP on the sensor surface. This can result in a type of interference in the dose–response of the immunoassay known as the hook effect where high concentrations of analyte produce lower signals (Tate & Ward, 2004). This phenomenon is frequently observed in one-sided immunoassays, such as the one considered here, and is due to the competition for available sensor binding-sites between active PMPs and target-antigen molecules. 2.1. The bio-sensor system The magnetometer-based biosensor system we are modelling is essentially an LC circuit and was described by Richardson et al. (2001) who used circuit theory to show that the output signal of the magnetometer is linearly dependent on the number of immobilized PMPs. Their argument can be summarized as follows. Consider a planar inductor, with self-inductance |$L_0$|⁠; if there are |$n$| immobilized PMPs on the sensor then its inductance will be given by $$ L = L_0 + kn, $$ where |$k$| is a constant of proportionality. The change in self-inductance of the inductor is due to the net magnetic field from the immobilized PMPs which interacts with the magnetic field of the inductor. From circuit theory, an LC circuit has a natural frequency that is dependent upon the inductance via a relationship of the form $$ f_0 = (L_0 C)^{-\frac{1}{2}}. $$ Therefore, if there are |$n$| bound particles, the frequency of the LC circuit can be written as $$ f_n = (L C)^{-\frac{1}{2}} = ((L_0 + kn) C)^{-\frac{1}{2}} = \left(L_0 C\right)^{-\frac{1}{2}} \left(1 + \frac{kn}{L_0}\right)^{-\frac{1}{2}}. $$ Assuming that |$kn/L_0 \ll 1$|⁠, we obtain $$ f_n \approx f_0 \left(1 - \frac{1}{2} \frac{kn}{L_0}\right)\!, $$ and so the change in the natural frequency is given by \begin{equation} {\it{\Delta}} f = f_n - f_0 \approx - f_0\frac{1}{2} \frac{kn}{L_0} . \label{eq:delf} \end{equation} (2.1) The magnetometer has been designed such that, over the active frequency range for the experiment, there is a positive linear relationship between the change in the LC circuit’s natural frequency and the change in magnetometer output voltage. If all experiments are conducted within this frequency range, then the change in voltage is proportional to the change in the natural frequency, \begin{equation} |{\it{\Delta}} V| \propto |{\it{\Delta}} f| \label{eq:delv}. \end{equation} (2.2) Using (2.1) and (2.2), we can say that the change in voltage output has a linear dependence upon the number of immobilized PMPs: $$ |{\it{\Delta}} V| \propto n. $$ 2.2. Magnetometer output signal Figure 2 shows a typical time profile for the magnetometer output signal during the measurement process in the immunoassay. After a lead time |$t_L$|⁠, which allows the assay to settle, the magnetic field is applied for a fixed time time |$t_\mathrm{ON}$|⁠, resulting in a step change in the magnetometer output signal. The magnetic field is then turned off for a time |$t_\mathrm{OFF}$|⁠. This process is repeated several times (typically about three) before a final relaxation time |$t_\mathrm{R}$|⁠. The magnetic field is ‘turned off’ by moving a permanent magnet between two positions, one in close proximity to the sensor surface |$z_0$| and the other at a set distance |$z_0+z$| from the sensor surface such that there is negligible effect from its magnetic field on the PMP dynamics. Fig. 2. Open in new tabDownload slide Schematic time profile for the magnetic field at the sensor surface during an assay. Fig. 2. Open in new tabDownload slide Schematic time profile for the magnetic field at the sensor surface during an assay. Figure 3(a) shows the observed time course of the magnetometer output signal during a particular experimental immunoassay measurement. Note how this output shares some of the signal features illustrated in Fig. 2. The large step changes in the output signal are produced by the movement of the magnet. However, following application of the magnetic field, there are additional, small but measurable changes in the voltage output. It is these changes that are used to measure the number of immobilized PMPs for a given set of assay parameters and target concentration. Fig. 3. Open in new tabDownload slide (a) Typical observed voltage output signal over time from the magnetometer during an assay. (b, c) The lower and upper parts of the signal in (a) plotted on a finer y-axis scale. Fig. 3. Open in new tabDownload slide (a) Typical observed voltage output signal over time from the magnetometer during an assay. (b, c) The lower and upper parts of the signal in (a) plotted on a finer y-axis scale. In Fig. 3(b,c), we have plotted the lower and upper parts of the signal on a finer scale, thus ignoring the large jump in the signal due to the change in the external magnetic field. For the lower part of the signal in (b), notice that the voltage varies nonlinearly with time when the magnetic field is present. Despite the potential interest of these dynamics, similar nonlinearities are observed during assays both with and without PMPs. Therefore, we cannot attribute this trend to the dynamics of binding particles. For the upper part of the signal in (c), notice that the voltage monotonically decreases with time when the magnetic field is not present. In addition, note that the gradient of both these upper and lower portions of the signal decrease with each pull of the magnet. In order to produce a dose–response curve from such a signal, an appropriate difference in voltage must be computed from the signal to measure the number of bound particles for a particular concentration of target. For a given set of assay parameters, the signal (response) will vary due to the change in concentration of the target (dose). By obtaining measurements for multiple different concentrations of target, we can plot a dose–response curve. There are several time instances within the signal that could potentially be used to measure a change in voltage. As in Sharif et al. (2013), the voltage change |${\it{\Delta}} V$| is taken as the difference between the initial voltage and the voltage immediately preceding the second application of the magnetic field, i.e. |${\it{\Delta}} V = V_4 - V_0$| (Fig. 4). Fig. 4. Open in new tabDownload slide Schematic of a typical output signal illustrating the measurement of the voltage change |${\it{\Delta}} V$| for a given concentration. Fig. 4. Open in new tabDownload slide Schematic of a typical output signal illustrating the measurement of the voltage change |${\it{\Delta}} V$| for a given concentration. Figure 5 shows three such experimental dose–response curves. These represent a combined analysis of historical and recently collected data and are produced here in order to be illustrative of typical responses; a systematic quantitative analysis is left for future work. The |$y$|-axis of each point on this plot represents a measure of the signal |${\it{\Delta}} V = V_4 - V_0$|⁠. Fig. 5. Open in new tabDownload slide Typical observed dose–response curves: (a) dose response from an immunoassay of C-reactive protein (CRP) marker |$(n=2$|⁠, curve is double exponential fit), (b) dose response from an immunoassay of CRP marker demonstrating the hook effect |$(n=1$|⁠, curve is double exponential fit) and (c) dose response from immunoassay of neutrophil elastase (NE) marker |$(n=2$|⁠, curve is linear fit). All error bars are |$\pm$| SEM. Fig. 5. Open in new tabDownload slide Typical observed dose–response curves: (a) dose response from an immunoassay of C-reactive protein (CRP) marker |$(n=2$|⁠, curve is double exponential fit), (b) dose response from an immunoassay of CRP marker demonstrating the hook effect |$(n=1$|⁠, curve is double exponential fit) and (c) dose response from immunoassay of neutrophil elastase (NE) marker |$(n=2$|⁠, curve is linear fit). All error bars are |$\pm$| SEM. Note that the dose–response in panel Fig. 5(a) shows a large jump in the voltage signal from the control condition (no antigen present in sample) followed by an approximately linear increase as the target concentration is increased. By contrast, the dose–response in panel Fig. 5(b) (from the same analyte as that in panel Fig. 5(a) but over a wider concentration range) illustrates the hook effect. Specifically, there is a continuous increase in response with dose, up to a maximum value followed by a decrease at larger analyte concentrations. The dose–response in panel Fig. 5(c), is for a different analyte and demonstrates a linear response over the range of tested analyte concentrations. In summary, experiments show that observed magnetometer dose–response curves can take various forms and are dependent upon the particular analyte and analyte concentration range to be measured. The purpose of constructing a mathematical model of this magnetic immunoassay is to inform the design and optimization of future assays for new biomarkers. Most importantly, in any design we must ensure that any calibrated dose–response is sensitive to the clinically relevant concentration range for that biomarker. We would also like to be able to explain and quantify the hook effect. 3. A computational model As an initial model, we have constructed a particle-dynamics simulation that considers the physical motion of the PMPs only and uses a simplified probabilistic binding model for the immobilization of the PMPs on the sensor surface. We aim to show that this approach is able to produce output signals and dose–response curves characteristic of those found experimentally in Figs 3 and 5. 3.1. PMP dynamics We consider the physical motion of the PMPs in the presence of a changing magnetic field. We assume that each PMP undergoes Brownian motion and so, in the absence of an external magnetic field, the motion of an individual PMP can be modelled using the stochastic ODE \begin{equation} m\ddot{\mathbf{x}} = - \lambda\dot{\mathbf{x}} + \epsilon \boldsymbol{\zeta}(t), \end{equation} (3.1) where |$\mathbf{x} = [x,y,z]^{\top}$|⁠, |$m$| is the particle mass, |$\lambda$| is the coefficient of the viscous force, |$\boldsymbol{\zeta}(t)$| is a 3D stationary Gaussian process where |$\langle \zeta_i(t) \rangle = 0$| and |$\langle \zeta_i(t) \zeta_i(t')\rangle = \delta(t-t')$| for |$i = 1,2,3$| where |$\delta(t)$| is the Dirac-delta function. According to Shevkoplyas et al. (2007), the additional force on a magnetic particle in an external magnetic field in the |$z$|-direction is |$\mathbf{F} = [0,0,F_z(z)]$| where \begin{equation} F_\mathrm{z}(z) = \frac{V_\mathrm{P} {\it{\Delta}} \chi}{\mu_0} B_\mathrm{z}(z) \frac{\partial}{\partial z} B_\mathrm{z}(z). \label{eq:fz} \end{equation} (3.2) Here, |$V_\mathrm{P}$| is the volume of an individual PMP, |${\it{\Delta}}\chi = \chi_\mathrm{P} - \chi_\mathrm{m}$| is the difference between the magnetic susceptibility of the particle |$\chi_\mathrm{P}$| and the medium |$\chi_\mathrm{m}$|⁠, |$\mu_0$| is the permeability of free space and |$B_\mathrm{z}(z)$| is the magnetic field along the axis of a permanent magnet as a function of height |$z$|⁠. From Biot–Savart theory (Griffiths, 1998), the latter can be modelled by the magnetic field produced by an equivalent circular current ring given by \begin{equation} B_\mathrm{z}(z) = \frac{B_0}{ [(z + z_0)^2 + a^2]^\frac{3}{2}}, \label{eq:Bz} \end{equation} (3.3) where |$z_0$| is the distance from the sensor surface to the centre point of the magnet, |$a$| is the radius of the magnet and |$B_0$| is a constant of proportionality. Differentiating this function with respect to |$z$|⁠, we get \begin{equation} \frac{\partial}{\partial z}B_\mathrm{z}(z) = - \frac{B_0 (z+z_0)}{ [(z + z_0)^2 + a^2]^\frac{5}{2}}. \label{eq:diffB} \end{equation} (3.4) Substituting (3.3) and (3.4) into (3.2), we can write the force on a PMP explicitly as a function of |$z$|⁠: $$ F_\mathrm{z}(z) = -c \frac{(z+z_0)}{ [(z + z_0)^2 + a^2]^4}, $$ where |$c = \frac{V_P {\it{\Delta}} \chi B_0^2}{\mu_0}$| is a constant of proportionality. In order to simplify this model, we consider the motion of a PMP in two dimensions, |$x$| and |$z$|⁠, such that its dynamics are modelled by \begin{equation*} \begin{aligned} m\ddot{x} & = - \lambda\dot{x} + \epsilon \zeta(t), \\ m\ddot{z} & = - \lambda\dot{z} + \epsilon \zeta(t) + \kappa(t) F_\mathrm{z}(z), \end{aligned} \end{equation*} where \begin{equation*} \kappa(t) = \begin{cases} 0 \qquad \mbox{when the magnet is}\,\, {\small{\text{OFF}}}, \\ 1 \qquad \mbox{when the magnet is}\,\, {\small{\text{ON}}}. \end{cases} \end{equation*} 3.2. Magnetization dynamics The magnetization of a PMP is dependent on its height |$z$| from the sensor surface (⁠|$z=0$| at the sensor surface). For simplicity, the magnetization of a PMP is assumed to be in one of three magnetic field-dependent states. These three states correspond to the magnetic field being any of: on, off after an on period, or off having never previously been applied. The magnetization of the particle is written as a vector |$\mathbf{M}(t) = |\mathbf{M}|(t)\hat{\mathbf{M}}(t) = M(t)\hat{\mathbf{M}}(t)$| where |$|\mathbf{M}|(t) = M(t)$|⁠. We present a heuristic model for the state of the magnetization by considering the strength of the magnetic field at the particle’s position |$z$|⁠, and how much rotation the particle experiences while in the assay. We assume that a PMP far from the magnet rotates more freely than a PMP close to the magnet and hence the random direction of that PMP’s magnetic dipole is less constrained at greater heights, |$z$|⁠. When the magnet is on, we assume that a PMP gains magnetization instantly and the subsequent dynamics for the magnetization are given by \begin{equation} \hat{\mathbf{M}}_\mathrm{on}(t) = \begin{pmatrix} \sin \theta_\mathrm{on}(t) \\ \cos \theta_\mathrm{on}(t) \end{pmatrix}, \qquad M_\mathrm{on}(t) \propto B_z(z(t)), \end{equation} (3.5) where \begin{equation} \theta_\mathrm{on}(t) = \xi_\mathrm{on} \left[1 - q\frac{B_z(z(t))}{B_z(0)}\right]\!, \end{equation} (3.6) and |$\xi_\mathrm{on}$| is a random number taken from the uniform distribution |$\mathrm{unif}\left(-\frac{\pi}{2},\frac{\pi}{2}\right)$| at each time-step. The choice of the value of |$q$| is arbitrary, but theoretically, it should be in the range |$0.5 < q < 1$|⁠. We choose |$q = 0.7$|⁠. The magnitude of the magnetization can then be written \begin{equation} M_\mathrm{on}(t) = M_0 \frac{B_z(z(t))}{B_z(0)}. \end{equation} (3.7) When the magnetic field is removed, the magnitude of the particle’s magnetization will decay over time. The orientation of the magnetization with respect to the |$z$|-axis is affected by the natural decay of the magnetization over time but will be affected more by the rotational motion of the particle. Another heuristic model can be proposed to model the dynamics when the magnetic field is not present. The magnitude of the magnetization is assumed to be \begin{equation} M_\mathrm{off}(t) = M_{t_r}\,{\rm e}^{-t/\tau}, \label{eq:magdec} \end{equation} (3.8) where |$M_{t_r}=M_\mathrm{on}(t_r)$|⁠, |$\tau$| is a decay constant and |$t$| is the time from when the magnetic field was removed (⁠|$t=t_r$|⁠). For the direction of the magnetization of a particle, \begin{equation} \hat{\mathbf{M}}_\mathrm{off}(t) = \begin{pmatrix} \sin \theta_\mathrm{off}(t) \\ \cos \theta_\mathrm{off}(t) \end{pmatrix} \end{equation} (3.9) with \begin{equation} \theta_\mathrm{off}(t) = \xi_\mathrm{off} \left(1-q \,{\rm e}^{-t/\tau}\right) \end{equation} (3.10) and |$\xi_\mathrm{off}$| is a random number taken from the distribution |$\mathrm{unif}(-\pi,\pi)$| at each time step. The values of all the parameters of the model can be found in Table 1. Table 1 Parameters used in the computational model Property . Symbol . Default value . Number of PMPs in assay |$P_T$| |$10^6$| Assay radius |$a$| |$5$||$\mathrm{mm}$| Sample diameter |$D$| |$0.05$||$\mathrm{mm}$| Assay height |$L$| |$2$||$\mathrm{mm}$| Viscous constant |$\lambda$| |$0.01$| Stochastic constant |$\epsilon$| |$0.01$| Magnet strength |$B_0$| |$10$| PMP mass |$m$| |$1$| Minimum distance from centre of magnet to assay surface |$z_0$| |$1$||$\mathrm{mm}$| Magnetization decay constant |$\tau$| |$10$||$\mathrm{s}$| Minimum normalized alignment of particle magnetization with field |$q$| |$0.7$| Property . Symbol . Default value . Number of PMPs in assay |$P_T$| |$10^6$| Assay radius |$a$| |$5$||$\mathrm{mm}$| Sample diameter |$D$| |$0.05$||$\mathrm{mm}$| Assay height |$L$| |$2$||$\mathrm{mm}$| Viscous constant |$\lambda$| |$0.01$| Stochastic constant |$\epsilon$| |$0.01$| Magnet strength |$B_0$| |$10$| PMP mass |$m$| |$1$| Minimum distance from centre of magnet to assay surface |$z_0$| |$1$||$\mathrm{mm}$| Magnetization decay constant |$\tau$| |$10$||$\mathrm{s}$| Minimum normalized alignment of particle magnetization with field |$q$| |$0.7$| Table 1 Parameters used in the computational model Property . Symbol . Default value . Number of PMPs in assay |$P_T$| |$10^6$| Assay radius |$a$| |$5$||$\mathrm{mm}$| Sample diameter |$D$| |$0.05$||$\mathrm{mm}$| Assay height |$L$| |$2$||$\mathrm{mm}$| Viscous constant |$\lambda$| |$0.01$| Stochastic constant |$\epsilon$| |$0.01$| Magnet strength |$B_0$| |$10$| PMP mass |$m$| |$1$| Minimum distance from centre of magnet to assay surface |$z_0$| |$1$||$\mathrm{mm}$| Magnetization decay constant |$\tau$| |$10$||$\mathrm{s}$| Minimum normalized alignment of particle magnetization with field |$q$| |$0.7$| Property . Symbol . Default value . Number of PMPs in assay |$P_T$| |$10^6$| Assay radius |$a$| |$5$||$\mathrm{mm}$| Sample diameter |$D$| |$0.05$||$\mathrm{mm}$| Assay height |$L$| |$2$||$\mathrm{mm}$| Viscous constant |$\lambda$| |$0.01$| Stochastic constant |$\epsilon$| |$0.01$| Magnet strength |$B_0$| |$10$| PMP mass |$m$| |$1$| Minimum distance from centre of magnet to assay surface |$z_0$| |$1$||$\mathrm{mm}$| Magnetization decay constant |$\tau$| |$10$||$\mathrm{s}$| Minimum normalized alignment of particle magnetization with field |$q$| |$0.7$| 3.3. Immobilization of PMPs A PMP is captured on the sensor surface through the cross-linking of bound antigens on the particle and antibodies on the sensor surface. Once a PMP is captured it is stationary and remains in the same orientation. If a particle is within a feasible distance |$\ell$| of the sensor surface such that it is physically possible for it to bind (i.e. the length of the protein can stretch from the surface of the PMP to the sensor surface) then we can attribute a non-zero probability for the particle binding to the sensor surface. This probability will be a function of many things including the target concentration, distribution of the target on the PMPs and density of the antibody on the assay surface. For the purposes of this simple model though, we consider that a particle has a fixed probability |$p$| of being immobilized on the assay surface if it is within a distance |$\ell$| of the assay surface. Therefore, at each time step, a random number |$x_\mathrm{rand}$| from the uniform distribution |$\mathrm{unif}(0,1)$| is generated and if |$x_\mathrm{rand} < p$| then the particle is immobilized, otherwise it is not immobilized during that time step. If a particle has been immobilized it maintains its current |$x$|-position for the remainder of the simulation, with its |$z$|-position fixed at |$z=0$|⁠. We assume that the direction of magnetization of an immobilized particle is \begin{equation} \hat{\mathbf{M}} = \begin{pmatrix} 0 \\ 1 \end{pmatrix} \end{equation} (3.11) and the magnitude of the magnetization is always |$M_0$|⁠, given by (3.8). 3.4. Simulation details Restricting ourselves to a computationally feasible number of particles, we consider a thin cylindrical region of the assay with radius |$D \ll a$| and height |$L$| where the radius and height of the reaction vessel are |$a$| and |$h$| respectively. The average number of PMPs in this cylinder |$n_p$| is given by \begin{equation} n_p = P_T \left(\frac{\pi D^2 h}{\pi a^2 h}\right) = P_T\left(\frac{D}{a} \right)^2, \end{equation} (3.12) where |$P_T$| is the total number of particles in the assay. For the dynamics of the PMPs in this thin cylinder, we assume, for simplicity, periodic boundary conditions in the |$x$|-direction and inelastic collisions with the assay surface at |$z=0$|⁠. The simulation was implemented in MATLAB and carried out on Blue Crystal, the University of Bristol supercomputer. 3.5. Results Figure 6 depicts a typical output signal from the computational model using the parameter values in Table 1. The plotted signal is an average of 2000 simulations. There are several features in this signal that are comparable to the observed signal in Fig. 3. For example, when the magnet is turned on or off, the signal undergoes an immediate large jump and during intervals when the magnet is on the signal increases, whereas when it is off the signal decreases. Another feature of the simulated signal that is similar to the observed signal is that for the times when the magnet is on the signal time profile displays significant nonlinearity and the slope decreases between subsequent on intervals. This is due to the gradual fixing of magnetic dipole orientations in the |$z$|-direction as PMPs become immobilized. Moreover, during the off periods, the initial gradient becomes more negative during subsequent intervals (compare e.g. |$50< t< 60$| and |$80< t< 90$|⁠). This is due to the Brownian motion of unbound particles drifting away from the sensor surface as well as the exponential decay of the particles’ magnetization. Fig. 6. Open in new tabDownload slide Typical time profile of the output signal from the computational model. In this case, the number of particles |$P_T=10^6$| and the probability of binding is |$p=0.019$|⁠. The output is a linear function of the net magnetic field |$B_z$| at |$z=0$|⁠, chosen so that the model output resembles the observed magnetometer voltage signal. Fig. 6. Open in new tabDownload slide Typical time profile of the output signal from the computational model. In this case, the number of particles |$P_T=10^6$| and the probability of binding is |$p=0.019$|⁠. The output is a linear function of the net magnetic field |$B_z$| at |$z=0$|⁠, chosen so that the model output resembles the observed magnetometer voltage signal. In addition, it can be observed in the simulation (data not shown) that the PMPs are attracted to the sensor surface over a shorter time scale (⁠|$<$|1 s) than displayed by the nonlinearities. Since the non-linearities are due to binding processes at the sensor surface, a key conclusion then is that the spatio-temporal effects of the magnet pull between on and off would seem to be rapid and so—assuming no aggregation of the PMPs—we can in principle make well-mixed assumptions during each of the magnet-on and magnet-off time intervals. 3.6. Dose–response curves We can use the model to produce dose–response curves by measuring the number of bound PMP particles on the surface after the first pull of the magnet, which serves as a proxy for the change in voltage which would be obtained experimentally. The number of PMPs immobilized at this time can be measured for different values of the probability |$p$|⁠, where the change in probability is a proxy for the change in the concentration of the target given the same assay conditions. In Fig. 7, we show how a dose–response curve can be produced as we vary the binding probability |$p$|⁠. There are two dose–response curves for two different assay conditions; where the number of PMPs in the assay is either |$P_T = 10^5$| or |$P_T = 10^6$|⁠. The model successfully produces a monotonically increasing signal, whose slope decreases at higher |$p$| values, but does not replicate the hook effect. The lack of a hook effect is because our modelling of the binding process as a single probabilistic coefficient is too simplistic to capture the confounding effect of the free antigen binding to the sensor surface at higher initial antigen concentration rates. Fig. 7. Open in new tabDownload slide Dose responses from the computational model. Fig. 7. Open in new tabDownload slide Dose responses from the computational model. Nevertheless, the overall effect of the computational model is to produce a Boltzman-like or sigmoidal response function, which is exactly what one would expect from a simple spatially lumped model that ignored the angular motion of the PMPs and any variation in the |$x$|-direction. Essentially it would seem that a spatially lumped model under a well-mixed assumption could easily produce the same results. 4. Rate-equation model We now consider separately the biochemical binding process of the target antigen to antibodies on the PMPs and the sensor surface, as well as the process by which PMPs become bound to the sensor surface. To do this, we use the law of mass action to derive rate equations of the key species within the assay, namely free PMPs, antigen-bound PMPs, surface-bound PMPs, free antigen, surface-bound antigen and free surface antibodies. Applying simple mass balance to each of these species leads to a coupled system of ODEs that can be simulated over the three experimental time periods. In the model, all antibody concentrations should be considered as effective antibody concentrations. This is because in reality, only a small fraction of antibodies immobilized on a particle or sensor surface will be capable of antigen capture. The majority of immobilized antibodies will be oriented such that their antigen binding sites are unavailable—it has been reported that only up to 4% of antibodies immobilized on a particle surface will actually be functional (Saha et al., 2014). 4.1. Dynamics during each experimental step We consider three distinct time intervals during the testing procedure; the incubation period |$0 \leq t< t_\mathrm{inc}$|⁠, the lead period |$t_\mathrm{inc} \leq t< t_2$|⁠, and the period that the magnet is applied for the first time |$t_2 \leq t< t_3$|⁠. We write |$t_2 = t_\mathrm{inc} + t_\mathrm{lead}$| and |$t_3 = t_\mathrm{inc} + t_\mathrm{lead} + t_\mathrm{ON}$| where |$t_\mathrm{inc}$|⁠, |$t_\mathrm{lead}$| and |$t_\mathrm{ON}$| are the durations of the incubation, lead and magnet-on periods, respectively. The total amount of simulated time is thus $$ t_\mathrm{end} = t_\mathrm{inc} + t_\mathrm{lead} + t_\mathrm{ON}. $$ In principle, we can then apply further periods of the magnet being off and on again, but we shall assume in what follows that the signal is proportional to the number of bound PMPs at the end of the first magnet on period. Antigen–antibody interactions can be described using the Langmuir kinetic model (Saha et al., 2014), which assumes binding occurs via a simple reversible reaction scheme. That is $$ {\rm A} + {\rm L} \mathop {\rightleftharpoons}\limits_{k^+}^{k^-} {\rm AL}, $$ where A and L are the antibody and ligand antigen respectively and |$k^+$| and |$k^-$| are the forward and reverse rate constants for the reaction. Antigen–antibody binding generally occurs with high affinity (⁠|$k^-/k^+ \ll 1$|⁠), therefore to simplify analysis of the model we consider antigen–antibody binding to be irreversible and use the reaction scheme \begin{equation} {\rm A} + {\rm L} \xrightarrow {k_0} {\rm AL}, \label{eq:chem} \end{equation} (4.1) which proceeds according to some rate constant |$k_0$|⁠. We define the model variable |$b$| to be the average concentration of antigen-bound PMP antibodies in the assay if only one PMP were present. The total concentration of PMP binding sites in the assay is therefore given by |$P_\mathrm{T}b$| where |$P_\mathrm{T}$| is the total number of PMPs in the assay. According to the reaction scheme (4.1), and using mass-action kinetics, the rate of change of |$b$| and the free target concentration |$T$| during the incubation period are given by \begin{equation} \begin{aligned} \dot{b}&= k_0 T(A - b), \\ \dot{T}&= -k_0 {P_\mathrm{T}} T(A - b) , \end{aligned} \label{eq:modp1} \end{equation} (4.2) where |$A >0$| is the average capture-antibody concentration if only one PMP were present and |${P_\mathrm{T}}>0$| is the total number of PMPs in the assay. For initial conditions |$T(0) = T_0$| and |$b(0) = 0$|⁠, (4.2) have the solution $$ \begin{aligned} b(t)&= \frac{A{T_0}(\exp([{T_0} - A{P_\mathrm{T}}]k_0 t]) - 1)}{{T_0} \exp([{T_0} - A{P_\mathrm{T}}]k_0 t) - A{P_\mathrm{T}}},\\ T(t)&= {T_0} -{P_\mathrm{T}} b(t), \end{aligned} $$ During the lead time, the solution containing the PMPs and target analyte is introduced to the sensor surface. In addition to a continuation of the dynamics occurring in the incubation period, we must consider the change of concentration of unbound antibody |$S$| on the sensor surface. We thus get a coupled ODE system of the form $$ \begin{aligned} \dot{b}&= k_0 T(A - b), \\ \dot{T}&= -k_0 P_T T(A - b) - k_1 S T, \\ \dot{S}&= - k_1 S T, \end{aligned} \label{eq:leadeq} $$ where |$k_1$| is a rate constant determining the velocity of the target-antigen/sensor-antibody binding reaction. For |$t > t_2$|⁠, when the magnetic field is on we must additionally consider the change in the number of PMPs |$P$| that become immobilized on the sensor surface. We also introduce a function |$f$| that places a threshold on the number of cross-links formed per PMP before PMPs start to immobilize. We explain the form of this function in more detail below. The ODE system now becomes \begin{equation} \begin{aligned} \dot{b}&= k_0 T (A - b), \\ \dot{T}&= -k_0 P_T T(A - b) - k_1 S T, \\ \dot{S}&= - k_1 S T - r k_2 (P_T - P) f(\alpha_1 b) f(\alpha_2 S),\\ \dot{P}&= k_2 (P_T - P) f(\alpha_1 b) f(\alpha_2 S), \end{aligned} \label{eq:magONode} \end{equation} (4.3) where |$k_2$| is a rate constant determining the velocity of the cross-linking reaction between an individual PMP-bound antigen and a surface antibody and |$r$|⁠, |$\alpha_1$| and |$\alpha_2$| are geometric constants that will be defined shortly. The function |$f$| takes the form of a linear relationship with a deadband for small argument: \begin{equation} f(x) := x H(x - Q), \label{eq:efffun} \end{equation} (4.4) where |$H$| is the Heaviside function $$ H(x - Q) = \begin{cases} 1&\mbox{if } x \geq Q \\ 0&\mbox{otherwise} . \end{cases} $$ Here the constant |$Q$| is the threshold value for |$\alpha_1b$| and |$\alpha_2S$| required for a PMP to become immobilized. This threshold is due to there being a minimum number of cross links per PMP for the PMP to bind securely to the sensor surface. |$Q$| will depend on a number of geometric factors. The key idea in defining |$r$|⁠, |$\alpha_1$|⁠, |$\alpha_2$| and |$Q$| is to recognize that the diameter of a PMP is an important factor in determining how many PMP binding-sites will be in physical proximity to the sensor antibodies once the PMP has moved into contact with the sensor surface. Binding sites on the opposite side of the PMP to that which is in contact, will not in general be able to bind to sensor antibodies. The parameters |$\alpha_1$| and |$\alpha_2$| are therefore geometric factors that determine which proportion of the PMP binding sites and the sensor antibodies can contribute to the rate of PMP immobilization (see Fig. 8). We also assume that the PMP mass density is uniform and constant across all PMPs. Therefore, the minimum force required to immobilize a PMP is proportional to its volume. In addition, we also assume that the force required to hold the PMP to the surface is proportional to the number of cross-links |$Q$|⁠. Fig. 8. Open in new tabDownload slide Schematic of a PMP immobilized on the sensor surface. Due to the PMP’s shape, binding is restricted to a binding zone of radius R (filled region at the PMP’s base). Fig. 8. Open in new tabDownload slide Schematic of a PMP immobilized on the sensor surface. Due to the PMP’s shape, binding is restricted to a binding zone of radius R (filled region at the PMP’s base). Given these assumptions we obtain, \begin{equation} Q = Q_0 d^3, \end{equation} (4.5) where |$d$| is the diameter of the PMP and |$Q_0$| is a constant of proportionality which depends on the density of the PMPs. From experimental experience the number of cross-links required to immobilize a typical PMP is |$\sim 10^2$| (Barnett et al., 2014; Sharif et al., 2013) which, given the other parameters in Table 3, corresponds to |$Q \sim \frac{10^2}{\mathrm{N_A}V_\mathrm{a}}\sim 10^{-12}\mu$|M, where N|$_\mathrm{A}$| and |$V_\mathrm{a}$| are Avogadro’s constant and the assay volume, respectively. Therefore the expected order of magnitude for |$Q_0 \sim 10^{-12} / (10^{-6})^3 \sim 10^{6}$|⁠. Table 3 Parameters in the mathematical model with the values used in this article. The values used are typical to an order of magnitude Parameter . Meaning . Value . Units . |$k_0$| Rate constant for target binding to PMPs |$1$| |$\mu$|M|$^{-1}$|s|$^{-1}$| |$k_1$| Rate constant for target binding to Ab on sensor surface 0.1 |$\mu$|M|$^{-1}$|s|$^{-1}$| |$k_2$| Rate constant for PMP binding to Ab on sensor surface (magnet on) |$10^{19}$| |$\mu$|M|$^{-2}$|s|$^{-1}$| |$A$| Average Ab concentration per PMP |$4 \times 10^{-10}$| |$\mu$|M |$P_T$| Total number of PMPs in assay |$10^{6}$| — |$h$| Binding zone depth |$10^{-8}$| m |$d$| Diameter of particle |$10^{-6}$| m |$a$| Radius of assay |$5\times 10^{-3}$| m |$Q_0$| Constant of proportionality |$2\times 10^{6}$| |$\mu$|M|$^{-1}$|m|$^{-3}$| |$t_\mathrm{inc}$| Incubation time |$60$| s |$t_\mathrm{lead}$| Lead time |$60$| s |$t_\mathrm{ON}$| Time magnet is on |$20$| s Parameter . Meaning . Value . Units . |$k_0$| Rate constant for target binding to PMPs |$1$| |$\mu$|M|$^{-1}$|s|$^{-1}$| |$k_1$| Rate constant for target binding to Ab on sensor surface 0.1 |$\mu$|M|$^{-1}$|s|$^{-1}$| |$k_2$| Rate constant for PMP binding to Ab on sensor surface (magnet on) |$10^{19}$| |$\mu$|M|$^{-2}$|s|$^{-1}$| |$A$| Average Ab concentration per PMP |$4 \times 10^{-10}$| |$\mu$|M |$P_T$| Total number of PMPs in assay |$10^{6}$| — |$h$| Binding zone depth |$10^{-8}$| m |$d$| Diameter of particle |$10^{-6}$| m |$a$| Radius of assay |$5\times 10^{-3}$| m |$Q_0$| Constant of proportionality |$2\times 10^{6}$| |$\mu$|M|$^{-1}$|m|$^{-3}$| |$t_\mathrm{inc}$| Incubation time |$60$| s |$t_\mathrm{lead}$| Lead time |$60$| s |$t_\mathrm{ON}$| Time magnet is on |$20$| s Table 3 Parameters in the mathematical model with the values used in this article. The values used are typical to an order of magnitude Parameter . Meaning . Value . Units . |$k_0$| Rate constant for target binding to PMPs |$1$| |$\mu$|M|$^{-1}$|s|$^{-1}$| |$k_1$| Rate constant for target binding to Ab on sensor surface 0.1 |$\mu$|M|$^{-1}$|s|$^{-1}$| |$k_2$| Rate constant for PMP binding to Ab on sensor surface (magnet on) |$10^{19}$| |$\mu$|M|$^{-2}$|s|$^{-1}$| |$A$| Average Ab concentration per PMP |$4 \times 10^{-10}$| |$\mu$|M |$P_T$| Total number of PMPs in assay |$10^{6}$| — |$h$| Binding zone depth |$10^{-8}$| m |$d$| Diameter of particle |$10^{-6}$| m |$a$| Radius of assay |$5\times 10^{-3}$| m |$Q_0$| Constant of proportionality |$2\times 10^{6}$| |$\mu$|M|$^{-1}$|m|$^{-3}$| |$t_\mathrm{inc}$| Incubation time |$60$| s |$t_\mathrm{lead}$| Lead time |$60$| s |$t_\mathrm{ON}$| Time magnet is on |$20$| s Parameter . Meaning . Value . Units . |$k_0$| Rate constant for target binding to PMPs |$1$| |$\mu$|M|$^{-1}$|s|$^{-1}$| |$k_1$| Rate constant for target binding to Ab on sensor surface 0.1 |$\mu$|M|$^{-1}$|s|$^{-1}$| |$k_2$| Rate constant for PMP binding to Ab on sensor surface (magnet on) |$10^{19}$| |$\mu$|M|$^{-2}$|s|$^{-1}$| |$A$| Average Ab concentration per PMP |$4 \times 10^{-10}$| |$\mu$|M |$P_T$| Total number of PMPs in assay |$10^{6}$| — |$h$| Binding zone depth |$10^{-8}$| m |$d$| Diameter of particle |$10^{-6}$| m |$a$| Radius of assay |$5\times 10^{-3}$| m |$Q_0$| Constant of proportionality |$2\times 10^{6}$| |$\mu$|M|$^{-1}$|m|$^{-3}$| |$t_\mathrm{inc}$| Incubation time |$60$| s |$t_\mathrm{lead}$| Lead time |$60$| s |$t_\mathrm{ON}$| Time magnet is on |$20$| s The constant |$\alpha_1$| represents the proportion of the binding sites on a PMP that contribute to PMP immobilization. This corresponds to the ratio of the surface area of a region of a PMP we shall call the binding zone (illustrated by the filled region at the PMP’s base in Fig. 8) and the total surface area of the PMP. Geometrically, the binding zone is a spherical cap with depth |$h$| on a PMP with diameter |$d$|⁠, thus the area of the binding zone is given by |$\pi h d$|⁠. Hence we can write $$ \alpha_1 = \frac{\pi d h}{4\pi (d/2)^2} =\frac{h}{d}. $$ The constant |$\alpha_2$|⁠, represents the proportion of sensor-surface antibodies that contribute to the PMP immobilization. It is the ratio of the circular area |$\pi R^2$| to the area of the assay surface (see Fig. 8). Therefore, we have \begin{equation} \alpha_2 = \frac{\pi R^2}{\pi a^2}, \label{eq:alpha2} \end{equation} (4.6) where |$a$| is the radius of the vessel. From Pythagoras’s theorem, |$R^2 = h(d-h)$|⁠, and substituting this expression into (4.6) we get $$ \alpha_2 = \frac{h(d-h)}{a^2}. $$ Additionally, when a PMP has attached to the sensor surface, the average concentration of binding sites on the sensor surface that are no longer active is given by |$r = S_2 \frac{\left(d/2\right)^2}{a^2}$| where |$S_2 = S(t_2)$|⁠. The variables and parameters used in this model are summarized in Tables 2 and 3. Table 2 Variables in the mathematical model and their initial values Variable . Meaning . Initial value . Units . |$b$| Volume concentration of antigen-bound antibodies on one PMP |$b_0 = b(t_0) = 0$| |$\mu$|M |$T$| Volume concentration of free Target in the assay |$T_0 = T(t_0) = 0.3$| |$\mu$|M |$S$| Free Ab binding sites on assay surface as volume concentration |$S_0 = S(t_0) = 0.1$| |$\mu$|M |$P$| Number of bound particles on surface |$P_0 = P(t_0) = 0$| — Variable . Meaning . Initial value . Units . |$b$| Volume concentration of antigen-bound antibodies on one PMP |$b_0 = b(t_0) = 0$| |$\mu$|M |$T$| Volume concentration of free Target in the assay |$T_0 = T(t_0) = 0.3$| |$\mu$|M |$S$| Free Ab binding sites on assay surface as volume concentration |$S_0 = S(t_0) = 0.1$| |$\mu$|M |$P$| Number of bound particles on surface |$P_0 = P(t_0) = 0$| — Table 2 Variables in the mathematical model and their initial values Variable . Meaning . Initial value . Units . |$b$| Volume concentration of antigen-bound antibodies on one PMP |$b_0 = b(t_0) = 0$| |$\mu$|M |$T$| Volume concentration of free Target in the assay |$T_0 = T(t_0) = 0.3$| |$\mu$|M |$S$| Free Ab binding sites on assay surface as volume concentration |$S_0 = S(t_0) = 0.1$| |$\mu$|M |$P$| Number of bound particles on surface |$P_0 = P(t_0) = 0$| — Variable . Meaning . Initial value . Units . |$b$| Volume concentration of antigen-bound antibodies on one PMP |$b_0 = b(t_0) = 0$| |$\mu$|M |$T$| Volume concentration of free Target in the assay |$T_0 = T(t_0) = 0.3$| |$\mu$|M |$S$| Free Ab binding sites on assay surface as volume concentration |$S_0 = S(t_0) = 0.1$| |$\mu$|M |$P$| Number of bound particles on surface |$P_0 = P(t_0) = 0$| — 4.2. Non-dimensionalization The rate-equation model can be non-dimensionalized through selection of appropriate scales as follows: \begin{equation} b^* = \frac{b}{A} \qquad T^* = \frac{T}{T_0} \qquad S^* = \frac{S}{S_0} \qquad P^* = \frac{P}{P_T} \qquad t^* = \frac{1}{k_1S_0}. \label{eq:og_sys} \end{equation} (4.7) The characteristic scale factors we have chosen for the dependent variables are the maximum values determined by the relevant model parameter. Substituting the scales into the original model and dropping the asterisk notation for the remainder of this article, yields the following non-dimensionalized model. During the incubation period, (⁠|$0 \leq t < t_1$|⁠), where |$t_1 = t_\mathrm{inc}$| we have \begin{align} \dot{b} &= \kappa_0\tau_0T(1-b), \nonumber \\ \dot{T} &= -\kappa_0\nu T(1-b). \label{eq:incnd} \end{align} (4.8) During the lead period, (⁠|$t_1 \leq t < t_2$|⁠), where |$t_2-t_1 = t_\mathrm{lead}$| we have \begin{align} \dot{b} &= \kappa_0\tau_0T(1-b), \nonumber \\ \dot{T} &= -\kappa_0\nu T(1-b) - ST, \nonumber \\ \dot{S} &= -\tau_0ST . \label{eq:lead} \end{align} (4.9) Finally, during the magnet-on period, (⁠|$t_2 \leq t < t_\mathrm{end}$|⁠), where |$t_\mathrm{end}-t_2 = t_\mathrm{ON}$| we have \begin{align} \dot{b} &= \kappa_0\tau_0T(1-b), \nonumber \\ \dot{T} &= -\kappa_0\nu T(1-b) - ST, \nonumber \\ \dot{S} &= -\tau_0ST - \kappa_2\eta(1-P)\hat{f}_1(b)\hat{f}_2(S), \nonumber \\ \dot{P} &= \kappa_2(1-P)\hat{f}_1(b)\hat{f}_2(S). \label{eq:magon} \end{align} (4.10) The function |$f$| from the unscaled model is replaced by |$\hat{f}_1$| and |$\hat{f}_2$| which are defined \begin{equation} \hat{f}_1(x) := x H(x - \chi_1), \qquad \hat{f}_2(x) := x H(x - \chi_2) \label{eq:efffun1} \end{equation} (4.11) where \begin{equation} \chi_1 = \frac{Q_0d^3}{\alpha_1A}, \qquad \chi_2 = \frac{Q_0d^3}{\alpha_2S_0}. \end{equation} (4.12) Finally, dimensionless parameter groups are given by \begin{equation} \kappa_0 = \frac{k_0}{k_1}, \qquad \tau_0 = \frac{T_0}{S_0}, \qquad \nu = \frac{P_TA}{S_0}, \qquad \kappa_2 = \frac{k_2\alpha_1\alpha_2A}{k_1}, \qquad \eta = \frac{rP_T}{S_0}. \end{equation} (4.13) The values of the dimensionless parameters according to the default parameter values in Table 3 are given in Table 4. Here |$\kappa_0$| is a parameter that determines the rate of antigen/PMP binding relative to antigen/sensor binding during the incubation period. The parameter |$\tau_0$| determines the initial amount of target antigen relative to the initial amount of sensor antibody. The ratio |$\nu$| is the total amount of PMP antibody in the assay divided by the total amount of sensor antibody. The parameter |$\kappa_2$| determines the rate of PMP immobilization at the sensor surface relative to the antigen/sensor binding rate. Finally, |$\eta$| determines the proportion of sensor antibody (at |$t=t_2$|⁠) which is obscured if all PMPs are immobilized. The parameter |$\eta$| is a slightly complicated parameter group since it contains the parameter |$r = S_2\frac{\left(d/2\right)^2}{a^2}$| and so can only be calculated after the lead period has completed. For physically realistic parameter values, |$\eta$| will be very small. Table 4 Dimensionless parameters in the mathematical model with the default values used in this article Parameter . Physical interpretation . Value . |$\kappa_0$| Target antigen/PMP Ab binding rate relative to target/sensor binding rate |$10$| |$\tau_0$| Initial target antigen relative to initial sensor antibody 0.33 |$\nu$| Total PMP antibodies relative to total sensor antibodies |$0.0014$| |$\kappa_2$| Rate of PMP immobilization relative to rate of target to sensor binding |$0.064$| |$\eta$| Ratio of sensor occupancy for complete PMP immobilization to initial sensor |$0.04\times S_2$| |$\chi_1$|⁠, |$\chi_2$| Scaled versions of threshold value for PMP immobilization |$Q$| 0.5, 0.013 Parameter . Physical interpretation . Value . |$\kappa_0$| Target antigen/PMP Ab binding rate relative to target/sensor binding rate |$10$| |$\tau_0$| Initial target antigen relative to initial sensor antibody 0.33 |$\nu$| Total PMP antibodies relative to total sensor antibodies |$0.0014$| |$\kappa_2$| Rate of PMP immobilization relative to rate of target to sensor binding |$0.064$| |$\eta$| Ratio of sensor occupancy for complete PMP immobilization to initial sensor |$0.04\times S_2$| |$\chi_1$|⁠, |$\chi_2$| Scaled versions of threshold value for PMP immobilization |$Q$| 0.5, 0.013 Table 4 Dimensionless parameters in the mathematical model with the default values used in this article Parameter . Physical interpretation . Value . |$\kappa_0$| Target antigen/PMP Ab binding rate relative to target/sensor binding rate |$10$| |$\tau_0$| Initial target antigen relative to initial sensor antibody 0.33 |$\nu$| Total PMP antibodies relative to total sensor antibodies |$0.0014$| |$\kappa_2$| Rate of PMP immobilization relative to rate of target to sensor binding |$0.064$| |$\eta$| Ratio of sensor occupancy for complete PMP immobilization to initial sensor |$0.04\times S_2$| |$\chi_1$|⁠, |$\chi_2$| Scaled versions of threshold value for PMP immobilization |$Q$| 0.5, 0.013 Parameter . Physical interpretation . Value . |$\kappa_0$| Target antigen/PMP Ab binding rate relative to target/sensor binding rate |$10$| |$\tau_0$| Initial target antigen relative to initial sensor antibody 0.33 |$\nu$| Total PMP antibodies relative to total sensor antibodies |$0.0014$| |$\kappa_2$| Rate of PMP immobilization relative to rate of target to sensor binding |$0.064$| |$\eta$| Ratio of sensor occupancy for complete PMP immobilization to initial sensor |$0.04\times S_2$| |$\chi_1$|⁠, |$\chi_2$| Scaled versions of threshold value for PMP immobilization |$Q$| 0.5, 0.013 From Section 2.1 we know that the voltage response from the magnetometer is linearly related to the number of immobilized PMPs at time |$t=t_\mathrm{end}$|⁠. Therefore, we can treat the predicted number of PMPs immobilized on the surface after time |$t=t_\mathrm{end}$| as a proxy for the observed change in voltage in a real experiment (see Fig. 4, data point |$(t_4,V_4)$|⁠). The change in the voltage signal from the magnetometer is due only to the net magnetic field from the immobilized PMPs since the free PMPs diffuse away from the sensor surface and therefore have negligible contribution to the change in self-inductance of the inductor. 5. Results We can simulate a dose–response curve for the immunoassay by numerically integrating equations (4.8)–(4.10) to determine |$P(t_\mathrm{end})$| over a range of values for |$\tau_0$|⁠. The results are shown in Fig. 9. These and all subsequent simulations were performed in MATLAB. Fig. 9. Open in new tabDownload slide Initial simulation results. (a–c): Simulated time courses of scaled model variables for different values of |$\tau_0$|⁠. Bottom: |$P(t_\mathrm{end})$| plotted against |$\tau_0$| to give a simulated dose–response curve (2000 data points). Fig. 9. Open in new tabDownload slide Initial simulation results. (a–c): Simulated time courses of scaled model variables for different values of |$\tau_0$|⁠. Bottom: |$P(t_\mathrm{end})$| plotted against |$\tau_0$| to give a simulated dose–response curve (2000 data points). The results show that with increasing |$\tau_0$| the model produces no response from |$P$| until a critical value |$\tau_0 = Q_0 d^4 P_T/(hS_0$|⁠) is reached. This is followed by an abrupt rise in |$P(t_\mathrm{end})$| and a continued increase as |$\tau_0$| increases until the amplitude of the response reaches a peak value. As the amount of initial target relative to sensor antibody becomes greater, |$P(t_\mathrm{end})$| decreases. This decrease is because at high |$\tau_0$|⁠, the number of free sensor binding sites at |$t=t_2$| will be small, reducing the likelihood of PMPs becoming immobilized. If the value of |$\tau_0$| is such that at |$t=t_2$|⁠, the number of binding sites on the sensor surface |$S(t_2) < Q/\alpha_2$| then it is not possible for PMPs to bind to the sensor surface and the response is zero. These observations demonstrate that the model is capable of replicating the hook effect, as illustrated in the lower panel of Fig. 9. 5.1. Asymptotic analysis We can derive an analytical approximation of the dose–response curve that is produced by the full model (4.8–4.10). We do this by letting |$t_1 \rightarrow \infty$| and approximating |$b$| and |$T$| at |$t=t_1$| by their asymptotic values. For example, Fig. 9(c) shows the average number of binding sites per particle |$b$| swiftly reaching saturation (⁠|$b=A$|⁠) within the incubation time |$t_\mathrm{inc}$|⁠. However, Fig. 9(a) shows that with lower values of |$\tau_0$|⁠, |$b$| does not saturate. Despite this, we will show that our analytical expression for a dose–response curve can provide an adequate approximation of the full model’s dose–response given certain conditions for the model parameter values. Our goal is to derive an expression for the number of bound particles on the sensor surface at the conclusion of the assay |$P(t=t_\mathrm{end})$|⁠. We proceed by considering the three time periods we are modelling in turn. Incubation period ($0 \leq t < t_1$): If we assume that the incubation period, |$t_\mathrm{inc}$|⁠, is long enough for |$b$| to saturate, then after the incubation time, we can use asymptotic values of these variables given by \begin{equation} \begin{aligned} \lim_{t_1 \to\infty} b(t_1) & \rightarrow \begin{cases} \tau_0/\nu &\mbox{if } \tau_0 < \nu, \\ 1 & \mbox{if } \tau_0 \geq \nu, \end{cases} \end{aligned} \label{eq:aymB} \end{equation} (5.1) and \begin{equation} \begin{aligned} \lim_{t_1 \to\infty} T(t_1) & \rightarrow \begin{cases} 0 &\mbox{if } \tau_0 < \nu ,\\ 1 - \nu/\tau_0 & \mbox{if } \tau_0 \geq \nu . \end{cases} \end{aligned} \label{eq:aymT} \end{equation} (5.2) Lead period ($t_1 \leq t < t_2$): During the lead time, we use the asymptotic values found in (5.1) and (5.2) as the initial conditions for (4.9). If |$\tau_0 < \nu$|⁠, then |$T(t_1)=0$| and thus there is no change in variables during the lead time. In this case, the initial conditions for the magnet-on period are the asymptotic values from the incubation period, |$[b(t_2), T(t_2), S(t_2), P(t_2)] = [\tau_0/\nu, 0, 1, 0]$|⁠. However, if |$\tau_0 \geq \nu$|⁠, then |$T$| and |$S$| will depend upon time. In this case, since |$b$| is saturated, (4.9) reduces to \begin{align} \dot{T} &= -ST, \end{align} (5.3) \begin{align} \dot{S} &= -\tau_0 ST, \end{align} (5.4) with initial conditions |$[T(t_1),S(t_1)]=[1-\nu/\tau_0, 1]$|⁠. The solution for (5.4) can be written in the form \begin{align} S(t) &= \frac{\delta}{\left(\delta + 1\right)\,{\rm e}^{\delta t} - 1}, \nonumber \\ T(t) &= \frac{1}{\tau_0}\left(S\left(t\right) + \delta\right)\!, \label{eq:sol_lead1} \end{align} (5.5) where |$\delta$| is a constant, which is given by $$ \delta = \tau_0 - \nu - 1. $$ Magnet-on period ($t_2 \leq t \leq t_\mathrm{end}$): The Heaviside functions present in (4.10) render them unsolvable in closed form. However, if they are removed (i.e. we write (4.11) as |$\hat{f}_i(x) = x$|⁠) and we consider the asymptotic approximations (4.3), the equations in (4.10) can be made analytically tractable. The ODE system for this period depends on relative size of |$\tau_0$| to |$\nu$|⁠. If |$\tau_0 < \nu$|⁠, then the initial conditions are |$[b(t_2), T(t_2), S(t_2), P(t_2)] = [\tau_0/\nu,0,1,0]$| and (4.10) reduces to \begin{align} \begin{split} \dot{S} &= -\kappa_2\eta\left(1-P\right)S, \\ \dot{P} &= \kappa_2(1-P)S. \label{eq:red_mag1} \end{split} \end{align} (5.6) We note that (5.6) is analytically solvable and that the solution can be written as \begin{align*} S(t) &= 1 - \eta P(t), \\ P(t) &= \frac{{\rm e}^{(\eta-1)\kappa_2\tau_0 t/\nu}-1}{\eta \,{\rm e}^{(\eta-1)\kappa_2\tau_0t/\nu}-1}, \end{align*} which gives us the expression for |$P(t)$| we wish to derive. However, if |$\tau_0 \geq \nu$|⁠, then initial conditions for |$T$| and |$S$| are calculated from (5.5) and the derivation of an expression for |$P(t)$| becomes more involved. In this case, the initial conditions are \begin{equation} [b(t_2), T(t_2), S(t_2), P(t_2)] = [1,T_2,S_2,0], \end{equation} (5.7) and the ODE system (4.10) reduces to \begin{align} \begin{split} \dot{T} &= -ST, \\ \dot{S} &= -\tau_0ST - \kappa_2\eta(1-P)S, \\ \dot{P} &= \kappa_2(1-P)S. \label{eq:sys3} \end{split} \end{align} (5.8) A solution for (5.8) cannot be found analytically. However, we can still find an expression of the form |$\dot{P} = f(P)$| by noting that \begin{equation} \dot{S} = \tau_0\dot{T} - \eta\dot{P}, \end{equation} (5.9) and therefore \begin{equation} S(t) = \tau_0T(t) - \eta P(t) + c, \label{eq:sub_S} \end{equation} (5.10) where |$c$| is a constant. Using the initial conditions at |$t = t_2$|⁠, we can write |$c = S_2 - \tau_0 T_2$|⁠. Now, dividing |$\dot S$| by |$\dot P$| from (5.8) gives \begin{equation} \frac{\rm dS}{\rm dP} = \frac{-\tau_0T - \kappa_2 \eta(1-P)}{\kappa_2 (1-P)}, \label{eq:sp} \end{equation} (5.11) and then eliminating |$T$| from (5.11) by substitution from (5.10) gives \begin{equation} \frac{{\rm d}S}{{\rm d}P} = \frac{c - S - \eta P - \kappa_2 \eta (1-P)}{\kappa_2 (1-P)}. \label{eq:odeSbyP} \end{equation} (5.12) Rearranging (5.12), we recognize it as a linear first-order ODE with initial condition |$S(0) = S_2$|⁠, \begin{equation} \frac{{\rm d}S}{{\rm d}P} + \frac{1}{\kappa_2 (1-P)}S = \frac{c - \eta P - \kappa_2 \eta (1-P)}{\kappa_2 (1-P)}. \label{eq:dSdP} \end{equation} (5.13) The solution for |$S$| as a function of |$P$| can be written \begin{equation} S(P) = c - \eta P + \tau_0T_2(1-P)^{1/\kappa_2} \label{eq:SofP}. \end{equation} (5.14) Substituting (5.14) into the expression for |$\dot P$| from (5.8) gives the time derivative of |$P$| written entirely as a function of |$P$| \begin{equation} \dot{P} = \kappa_2(1-P)\left(c - \eta P + \tau_0T_2(1-P)^{1/\kappa_2}\right)\!. \label{eq:nl} \end{equation} (5.15) The nonlinearities in (5.15) prevent it from being solvable in closed form. However, considering the full expression for |$c$| and using an alternative expression for |$\eta = S_2r_0P_T$| (where |$r_0 = \frac{(d/2)^2}{a^2}$|⁠) we can rewrite (5.15) as \begin{equation} \dot{P} = \kappa_2(1-P)\left(\tau_0T_2(1-P)^{1/\kappa_2} - \tau_0T_2 + S_2(1-r_0P_TP)\right)\!. \label{eq:lin} \end{equation} (5.16) Note that |$r_0P_TP \leq r_0P_T \sim 10^{-2}$| (where |$P_T \sim 10^6$|⁠) and thus |$1-r_0P_TP \approx 1$|⁠. Therefore, we can rewrite (5.16) as \begin{equation} \dot{P} \approx \kappa_2(1-P)\left(\tau_0T_2(1-P)^{1/\kappa_2} - \tau_0T_2 + S_2\right)\!, \end{equation} (5.17) which is analytically solvable; the solution is \begin{equation} P(t) \approx 1 - \left(\frac{S_2 - \tau_0T_2}{S_2\,{\rm e}^{ct} - \tau_0T_2}\right)^{\kappa_2}. \label{eq:anaPt} \end{equation} (5.18) Equation (5.18) is an analytical approximation for |$P$| as a function of time during the first pull of the magnet (⁠|$t_2 \leq t < t_3$|⁠). At the end of this period |$t = t_\mathrm{end}$|⁠, the value of |$P$| gives us the predicted scaled number of PMPs bound on the sensor surface |$P(t_\mathrm{end})$|⁠. Summary: An approximation for a dose–response curve from the full model equations (4.8)–(4.10) can be calculated as follows: When |$\tau_0 < \nu$|⁠: \begin{equation} P(t_\mathrm{ON}) = \frac{{\rm e}^{(\eta-1)\kappa_2\tau_0 t_\mathrm{ON}/\nu}-1}{\eta \,{\rm e}^{(\eta-1)\kappa_2\tau_0t_\mathrm{ON}/\nu}-1}, \label{eq:P1} \end{equation} (5.19) where |$t_\mathrm{ON}$| is the magnet-on time. When |$\tau_0 \geq \nu$|⁠: \begin{equation} P(T_2, S_2, t_\mathrm{ON})=P(T_2(t_\mathrm{lead}),S_2(t_\mathrm{lead}),t_\mathrm{ON})\approx 1 - \left(\frac{S_2 - \tau_0T_2}{S_2\,{\rm e}^{ct_\mathrm{ON}} - \tau_0T_2}\right)^{\kappa_2}, \label{eq:P2} \end{equation} (5.20) where \begin{align} \begin{split} S_2(t_\mathrm{lead}) &= \frac{\delta}{\left(\delta + 1\right)\,{\rm e}^{\delta t_\mathrm{lead}} - 1}, \nonumber \\ T_2(t_\mathrm{lead}) &= \frac{1}{\tau_0}\left(S_2\left(t_\mathrm{lead}\right) + \delta\right)\!, \end{split} \end{align} where |$t_\mathrm{lead}$| is the lead time and |$\delta = \tau_0 - \nu - 1$|⁠. 5.2. Dose–response curves By calculating |$P(t_\mathrm{end}$|⁠) for different values of |$\tau_0$| (using parameter values in Table 4 unless otherwise indicated) we are able to obtain dose–response curves for the model using either direct simulation or the asymptotic formulae. Figure 10 shows the dose–response curves over the range |$10^{-4} < \tau_0 < 10^{2}$| for the full model equations, the analytical approximation and a number of other modified versions of the full model for different values of |$\kappa_0$|⁠. The dose–response curves from the variations of the full model are included to show the effect of the assumptions made in the analytical approximation on the dose–response relationship. Fig. 10. Open in new tabDownload slide Simulated dose responses from the model. The non-dimensional model variable |$P$| plotted against the dimensionless parameter group |$\tau_0$| for variations of the model including the approximate analytical expression for |$P$|⁠, derived in Section 5.1. Model variants are explained in the main body text. Fig. 10. Open in new tabDownload slide Simulated dose responses from the model. The non-dimensional model variable |$P$| plotted against the dimensionless parameter group |$\tau_0$| for variations of the model including the approximate analytical expression for |$P$|⁠, derived in Section 5.1. Model variants are explained in the main body text. The ‘Full’ curve (plus sign markers) shows the dose–response relationship obtained by numerically integrating the full model (4.8–4.10) over the time period |$0 < t < t_\mathrm{end}$| and the ‘Analytical’ curve (square markers) shows the dose–response relationship obtained by calculating |$P(t_\mathrm{end})$| from (5.19) and (5.20). The ‘No HS’ curve (circle markers) shows the dose–response relationship when the Heaviside function in (4.10) is relaxed such that (4.11) becomes |$f(x)=x$|⁠. Comparing the ‘Full’ and ‘No HS’ curves shows that inclusion of the Heaviside function results in an abrupt increase in |$P$| as |$\tau_0$| is increased beyond a critical value. By contrast, relaxing the Heaviside function results in a gradual increase in |$P$| as |$\tau_0$| is increased. The analytical dose–response does not have the sensitive response from |$P$| around the critical value of |$\tau_0$| since the Heaviside functions are not included in its derivation. The ‘Full’ curve and ‘No HS’ curve approach the analytical approximation for the dose–response as |$\kappa_0$| is increased. This is because a larger value for |$\kappa_0$| corresponds to a larger value for the antigen to PMP antibody binding rate relative to the antigen to sensor antibody binding rate which therefore means that the PMP antibody saturation assumption in the derivation of the approximation becomes increasingly valid. The ‘Sat’ curve (asterisk markers) represents the dose–response of the full model equations but where we have assumed that the incubation time is sufficient for the average concentration of binding sites on a single PMP |$b$| and the amount of free target |$T$| is close to the limits they would reach as |$t \rightarrow \infty$|⁠. The ‘Sat, no HS’ model (cross markers) is a combination of the Sat and No HS models. Both curves illustrate how over a wide range of |$\tau_0$| the analytic approximation is pretty good given a long incubation time. As would be expected, the Sat, no HS model and the analytical model are in good agreement, at least over the range of model parameters explored. 5.3. Parameter sensitivity To investigate the sensitivity of the simulated dose–response to the other dimensionless parameter groups in the model, we ran simulations of the full model equations (corresponding to ‘Full’ in Fig. 10) across a range of parameter values. The results are shown in Fig. 11 (solid lines) as well as the corresponding analytical approximation (dashed lines). Fig. 11. Open in new tabDownload slide Simulated dose responses from numerical integration of the full equations (solid lines) and analytical approximation from Section 5.1 (dashed lines). The non-dimensional model variable |$P$| is plotted against the dimensionless parameter group |$\tau_0$| for different parameter values. In all simulations other than those in panel (a), |$\kappa_0=10^4$|⁠. (a) shows the convergence of the full model solution with the analytical approximation as |$\kappa_0$| is increased. Only one approximate solution is visible as it is not dependent upon |$\kappa_0$| and therefore solutions overlap. (b) Increasing the value of |$\nu$| reduces the sensitivity of the model response but does not effect the decreasing portion of the dose–response curve. (c) Increasing |$\kappa_2$| increases the amplitude of the model response. (d) |$\eta$| is a dimensionless parameter group which includes model variable |$S_2=S(t_2)$|⁠. Therefore, we systematically vary the value of |$\eta/S_2=\frac{(d/2)^2P_T}{a^2S_0}$| to investigate the effect on the model response. Within parameter values that would practically be used in experiment, there is no effect on the model response. Fig. 11. Open in new tabDownload slide Simulated dose responses from numerical integration of the full equations (solid lines) and analytical approximation from Section 5.1 (dashed lines). The non-dimensional model variable |$P$| is plotted against the dimensionless parameter group |$\tau_0$| for different parameter values. In all simulations other than those in panel (a), |$\kappa_0=10^4$|⁠. (a) shows the convergence of the full model solution with the analytical approximation as |$\kappa_0$| is increased. Only one approximate solution is visible as it is not dependent upon |$\kappa_0$| and therefore solutions overlap. (b) Increasing the value of |$\nu$| reduces the sensitivity of the model response but does not effect the decreasing portion of the dose–response curve. (c) Increasing |$\kappa_2$| increases the amplitude of the model response. (d) |$\eta$| is a dimensionless parameter group which includes model variable |$S_2=S(t_2)$|⁠. Therefore, we systematically vary the value of |$\eta/S_2=\frac{(d/2)^2P_T}{a^2S_0}$| to investigate the effect on the model response. Within parameter values that would practically be used in experiment, there is no effect on the model response. Figure 11(a) illustrates how the dose–response curve varies with the parameter |$\kappa_0$|⁠, which represents the ratio between target binding affinities of the PMP antibodies to that of the sensor surface antibodies. Similarly to Fig. 10, the simulated dose–response is in good agreement with the analytical approximation at large values of |$\kappa_0$|⁠. This is because as |$\kappa_0$| gets larger, the validity of the saturation assumption in the approximation is increased. Only one dashed line is visible since the analytical formula is not dependent on the parameter |$\kappa_0$|⁠. Note also that at higher values of |$\kappa_0$|⁠, there is more dose sensitivity at low |$\tau_0$|⁠; for lower |$\kappa_0$|-values, the dose–response curve represents a sudden switch from zero response to maximum. This would suggest that high |$\kappa_0$| would be a good regime for sensing, especially at low-dose levels. Figure 11(b) illustrates the how the dose–response curve varies with the parameter |$\nu$|⁠, which represents the ratio of PMP antibodies to sensor antibodies. The results show that the range of initial target concentration |$\tau_0$| that produces a response from the model increases as the value for |$\nu$| decreases. The switch-on point of the immunoassay increases as the number of PMP antibodies relative to sensor antibodies increases. This suggests that |$\nu$| could be a good parameter to vary in order to control the assay’s limit of detection. Figure 11(c) shows how the dose–response curve varies with the parameter |$\kappa_2$|⁠, which measures the ratio of the rate of immobilization to the rate of antigen to sensor antibody binding. Note that this dimensionless parameter is somewhat complex in that it depends on the geometric constants |$\alpha_1$| and |$\alpha_2$| and the density of the antibodies on each PMP. It would seem that for a given antibody–antigen binding affinity and PMP antibody concentration, this parameter can be treated as a proxy for PMP size, since |$\kappa_2$| increases linearly with PMP volume. The results show that as |$\kappa_2$| increases the magnitude of the model response increases, while the range of |$\tau_0$| to which the response is sensitive does not change. Nevertheless, treating this result as representing what happens as particle size increases is problematic, because the variables |$\chi_1$| and |$\chi_2$| would also increase with particle size, which means that there would presumably also be a higher value of |$\tau_0$| at which the sensor is first switched-on for larger particle sizes. It may therefore be appropriate to fix the particle size for a given antibody–antigen binding affinity and manipulate the amplitude of the model response by changing the average antibody concentration for a single PMP, parameter |$A$|⁠. Figure 11(d) shows that over the range of parameter values simulated (which incorporates the full range one would expect due to physical limitations on the real system) the dose response is insensitive to the somewhat complex parameter |$\eta$|⁠. In short, this suggests that if the number of PMPs is kept within a suitable range relative to the size of the sensor surface, then the obscuring of sensor antibodies by immobilized PMPs should not be a significant factor in the efficiency or sensitivity of the assay. 5.4. Comparison with experimental data In practical usage, for a particular target antigen, the clinically relevant range of |$\tau_0$| will likely be over one or two orders of magnitude. Here, we look at the form of simulated dose–response curve over narrower windows within the full range of values for |$\tau_0$| simulated. Figure 12(a) shows a simulated dose response along a linear |$x$|-axis and the different panels within Fig. 12(b) show narrower ranges of |$\tau_0$| at different points along the curve. Fig. 12. Open in new tabDownload slide (a) Full and No HS model simulated dose response plotted on a linear-scale |$x$|-axis. (b) Panels A–D each show detail of the dose–response curve plotted in (a) over different value ranges for |$\tau_0$|⁠. Fig. 12. Open in new tabDownload slide (a) Full and No HS model simulated dose response plotted on a linear-scale |$x$|-axis. (b) Panels A–D each show detail of the dose–response curve plotted in (a) over different value ranges for |$\tau_0$|⁠. Figure 12(b)(A) shows a region of the dose–response curve that demonstrates the model’s ability to replicate experimental dose–response curves such as that shown in Fig. 5(a). Here at very low concentration of target antigen there is zero response from the model, followed, as |$\tau_0$| increases, by a sudden switch-like, large increase in the response amplitude. This is then followed by a slow decrease in the signal for higher doses. Such a dose response would be good for an assay, such as pregnancy testing, where a binary response is required for the presence of antigen in significant quantity. Figure 12(b)(B) shows that the model is capable of generating a dose–response curve demonstrating the so-called hook effect—as target antigen concentrations become large, the competition between target antigen and PMP binding at the sensor antibody surface results in a decrease in the measured response. Such responses are generally to be avoided when designing immunoassays because there is no one-to-one mapping between output signal and target concentration. In general, we found that the model does not readily generate the type of rising linear response seen in experimental data from Fig. 5(c). Nevertheless, Fig. 12(b)(C) shows how over a narrow range of |$\tau_0$| the model could replicate a roughly linear response, but only if the Heaviside function is ignored. This would correspond in the model to choosing low values of |$\chi_1$| and |$\chi_2$|⁠, which would represent particularly ‘strong’ antibody–antigen pair binding meaning that fewer cross-links would be required to immobilize a PMP. Interestingly, the data in Fig. 5(c) contains arguably the largest uncertainty, as evidenced in the large error bar at 60 ng|$/$|ml. Figures 10 and 11 show that the most robust region of the dose–response curve produced by the full model is the decreasing region following the peak response to |$\tau_0$|⁠. In this region, the result is a near-linear, inverse relationship between dose and signal response, as shown by Fig. 12(b)(D). This region of the curve would therefore be the best region to use for the immunoassay if the main concern is to understand in a single assay the effect of dosage strength beyond an initial signal that represents presence. Such assays would be good for quantifying biomarkers indicative of an infection. The slightly perverse nature of this kind of assay though is that the strongest signal is obtained with the weakest detectable dose. Such a feature might not be desirable however, because a weak response could represent either a strong dose or simply the noise floor in the absence of any detectable target. Nevertheless, this form of assay is not susceptible to the hook effect and the curve of this shape seems resistant to calibration other than by manipulating the value of the initial sensor antibody concentration |$S_0$|⁠. 6. Conclusion This article has constructed for the first time an end-to-end mathematical model of a particular biosensor that is typical of well-mixed particle-based immunoassays. The need for modelling is clear in this case because of the multi-physics involved in the assay so that it is not obvious which effect dominates the dose response. We have accounted for particle diffusion, magnetization and immobilization, as well as different affinities of antigen-antibody binding and the geometric effects that occur when a particle is cross-linked. We first developed a computational particle-dynamics model. The results from this indicated that the movement of particles due to the magnetic field is rapid and dominates any diffusion effects. The model is able to replicate the shape of the voltage output signal over time during the assay (compare Figs 3 and 6). Moreover the model is able to offer an explanation for the growth and decay of the signal during the magnet on and off periods, respectively. The disadvantage of the particle-based model though is that it does not capture the chemistry of binding, other than via a simple probability of a particle being bound, and so it is not readily possible to explain the sensitivity of dose–response curves to assay parameters, nor to capture the hook effect. One key conclusion of the model though is that there are limited spatial effects, so that a well-mixed hypothesis can be used. To capture the varied shapes of observed dose–response curves and how they depend on assay parameters, we secondly produced a spatially-lumped rate-equation model based on the law of mass action. This ODE model was able to produce dose–response curves that emulate all shapes of dose–response curve seen in the experimental results in Fig. 5(b), including the hook effect. Another key contribution has been non-dimensionalization, which enables us to identify dimensionless groupings that can enable a scientific approach to parameter modification during experimental design. Moreover, under the conditions where |$r_0 P_T \ll 1$| and the dimensionless parameter |$\kappa_0$| is sufficiently large, we have shown how the output of the nonlinear ODE model can be well approximated by an analytical expression for an appropriate range of target concentration. Analytic expressions have the advantage over computational results in that they enable direct experimental design without the need for simulation. These results show promise because the qualitative effect on the immunoassay’s dose–response curve can be easily estimated by changing key assay parameters as demonstrated in Figs 10 and 11. We should stress though that the main focus of this paper has not been to precisely match experimental data. Rather we have merely sought to illustrate that the model is capable of reproducing the full range of qualitatively distinct dose–response curves that have been observed in the laboratory and to identify how assay parameters can affect these. In essence, we have merely shown that the model has descriptive power; whether or not it has predictive power remains to be seen. A detailed comparison between our model and a carefully designed experimental campaign, will form the subject of future work. We also would like to adapt the modelling principles developed here to seek to understand particle-based immunoassays that utilize microfluidic channels and lateral flow through porous media. In these cases, spatio-temporal reaction kinetics and fluid mechanic effects are likely to be significant. References Barnett J. , Wraith P., Kiely J., Persad R., Hurley K., Hawkins P. & Luxton R. ( 2014 ) An inexpensive, fast and sensitive quantitative lateral flow magneto-immunoassay for total prostate specific antigen. Biosensors, 4 , 204 – 220 Google Scholar Crossref Search ADS WorldCat Cole L. A. , Shahabi S., Butler S. A., Mitchell H., Newlands E. S., Behrman H. R. & Verrill H. L. ( 2001 ) Utility of commonly used commercial human chorionic gonadotropin immunoassays in the diagnosis and management of trophoblastic diseases. Clin. Chem., 47 , 308 – 315 . OpenURL Placeholder Text WorldCat Griffiths D. J. ( 1998 ) Introduction to Electrodynamics., 3rd edn. New Jersey : Prentice Hall . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Jamshaid T. , Neto E. T. T., Eissa M. M., Zine N., Kunita M. H., El-Salhi A. E. & Elaissari A. ( 2016 ) Magnetic particles: from preparation to lab-on-a-chip, biosensors, microsystems and microfluidics applications. Trends Anal. Chem., 79 , 344 – 362 . Google Scholar Crossref Search ADS WorldCat Kokalj T. , Pérez-Ruiz E. & Lammerytn J. ( 2015 ) Building bio-assays with magnetic particles on a digital microfluidic platform. New Biotechnol., 32 , 485 – 503 . Google Scholar Crossref Search ADS WorldCat Koh I. & Josephson L. ( 2009 ) Magnetic nanoparticle sensors. Sensors, 9 , 8130 – 8145 . Google Scholar Crossref Search ADS WorldCat Mackey D. , Kelly E. & Nooney R. ( 2016 ) Modelling random antibody adsorption and immunoassay activity. Math. Biosci. Eng., 13 , 1159 – 1168 . Google Scholar Crossref Search ADS WorldCat Munir A. , Wang J., Zhu Z. & Zhou H. S. ( 2011 ) Mathematical modeling and analysis of a magnetic nanoparticle-enhanced mixing in a microfluidic system using time-dependent magnetic field. IEEE Trans. Nanotechnol., 10 , 953 – 961 . Google Scholar Crossref Search ADS WorldCat Richardson J. , Hill A., Luxton R. & Hawkins P. ( 2001 ) A novel measuring system for the determination of paramagnetic particle labels for use in magneto-immunoassays. Biosens. Bioelectron., 16 , 1127 – 1132 . Google Scholar Crossref Search ADS WorldCat Saha B. , Evers T. H. & Prins M. W. J. ( 2011 ) How antibody surface coverage on nanoparticles determines the activity and kinetics of antigen capturing for biosensing. Anal. Chem., 86 , 8158 – 8166 . Google Scholar Crossref Search ADS WorldCat Sharif E. , Kiely J. & Luxton R. ( 2013 ) Novel immunoassay technique for rapid measurement of intracellular proteins using paramagnetic particles. J. Immunol. Meth., 388 , 78 – 85 . Google Scholar Crossref Search ADS WorldCat Shevkoplyas S. S. , Siegel A. C., Westervelt R. M., Prentiss M. G. & Whitesides G. M. ( 2007 ) The force acting on a superparamagnetic bead due to an applied magnetic field. Lab Chip ., 7 , 1294 – 1302 . Google Scholar Crossref Search ADS PubMed WorldCat Tate J. & Ward G. ( 2004 ) Interferences in immunoassay the force acting on a superparamagnetic bead due to an applied magnetic field. Clin. Biochem. Rev., 25 , 105 – 120 . OpenURL Placeholder Text WorldCat Xu H. , Jian R. L. & Williams D. E. ( 2006 ) Effect of surface packing density of interfacially adsorbed monoclonal antibody on the binding of hormonal antigen human chorionic gonadotrophin. The J. Phys. Chem. B ., 110 , 1907 – 1914 . Google Scholar Crossref Search ADS WorldCat © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.
TI - Mathematical modelling of a magnetic immunoassay
JF - IMA Journal of Applied Mathematics
DO - 10.1093/imamat/hxx034
DA - 2017-12-01
UR - https://www.deepdyve.com/lp/oxford-university-press/mathematical-modelling-of-a-magnetic-immunoassay-u3FxHwpdyh
SP - 1253
VL - 82
IS - 6
DP - DeepDyve
ER -