Brownian-motion based simulation of stochastic reaction–diffusion systems for affinity based sensors

In this work, we develop a 2D algorithm for stochastic reaction–diffusion systems describing the binding and unbinding of target molecules at the surfaces of affinity-based sensors. In particular, we simulate the detection of DNA oligomers using silicon-nanowire field-effect biosensors. Since these devices are uniform along the nanowire, two dimensions are sufficient to capture the kinetic effects features. The model combines a stochastic ordinary differential equation for the binding and unbinding of target molecules as well as a diffusion equation for their transport in the liquid. A Brownian-motion based algorithm simulates the diffusion process, which is linked to a stochastic-simulation algorithm for association at and dissociation from the surface. The simulation data show that the shape of the cross section of the sensor yields areas with significantly different target-molecule coverage. Different initial conditions are investigated as well in order to aid rational sensor design. A comparison of the association/hybridization behavior for different receptor densities allows optimization of the functionalization setup depending on the target-molecule density.


Introduction
The label-free, fast, and reliable detection of biomolecules in liquids has seen great experimental progress in recent years. Label-free detection of minuscule amounts of DNA oligomers and proteins has been demonstrated using nanowire field-effect sensors [1][2][3][4][5][6]. However, to enable the optimal sensor design, it is crucial to provide fundamental understanding of the characteristics of the sensing processes. The first step here is the determination of the rates at which surface interactions occur. This has been investigated by several groups employing different approaches, ranging from models based on the Poisson-Boltzmann equation [7], on the master-equation formalism [8] or on generalized Langmuir forms for the isotherms [9], but also in one of our previous works [10]. The most important feature that is covered by these approaches is the fact that the surface hybridization rate decreases with increasing receptor density at the surface. With this knowledge, topics of recent interest on the whole sensing process include the reduction of the noise levels [11][12][13][14][15][16], which enables lower detection limits, or the determination of waiting times until the biomolecules have been detected. Mathematical models for the average sensor response and their simulation including studies of the effect of device parameters have previously been reported in [17][18][19][20][21][22] for biomolecule detection and also for gas detection [23].
In order to obtain a more realistic model of the sensing process, it is not sufficient to consider the hybridization/ dissociation dynamics at the surface. Rather, it is crucial to take the limited transport of the DNA molecules through the liquid into account, as especially in the case of low targetmolecule concentrations, the hybridization characteristics are changed significantly by this process. A thorough investigation and a quantification of the resulting effects are the objective of this study. Recently, we presented a onedimensional model describing the surface interactions as stochastic processes and including a transport model for the biomolecules in the liquid [10,24]. The simulations of the surface processes for different parameters yielded some important, yet unexpected results for the optimal system design. This investigation is now continued with a twodimensional model, which allows the consideration of more complicated and hence more realistic geometries. Using a Brownian-motion approach for the diffusion of target molecules results in much faster simulations compared to the previous algorithm. The Brownian motion is coupled to a stochastic ordinary differential equation for binding and unbinding at the sensor surface. Partitioning the surface allows to distinguish different types of regions, where the hybridization/dissociation statistics differ quantitatively as well as qualitatively. In particular, this feature enables the analysis of the effects of off-sensor hybridization. These results yield important insights regarding the optimal design of the devices, in particular regarding favorable sensor geometries and functionalization densities. This paper is organized as follows. In section 2, the coupled model is described. The algorithm as implemented is presented in section 3. Numerical results can be found in section 4, where we investigate the binding and unbinding behavior at different areas of the device, for different initial conditions and also for different receptor densities. The findings are discussed in section 5.

The coupled model
The principle of biomolecule detection with nanowire sensors in field-effect transistor configuration is simple. The surface of the device, i.e.the nanowire and parts of the substrate, are functionalized with the receptor molecules matching the target molecules to be detected, which guarantees the selectivity of the device. The liquid solution surrounding the device is then enriched with the target molecules. A sketch of the crosssection of the device is shown in figure 1.
The hybridization of the target with its matching receptor at the device surface yields a change in the surface-charge density, which in turn influences the conductivity of the whole nanowire. Measuring the conductance of the device therefore enables detection of the biomolecules.
However, there are many effects arising in the whole system that are not yet well understood. Some of them lead to fluctuations in the signal, where the part stemming from biological noise shall be quantified by the simulations described in the following.

Association and dissociation
The dynamics at the surface are based on the association (or hybridization) of a free near-surface target molecule with a matching receptor molecule fixed at the surface, forming a hybridized complex, and on the inverse process of dissociation, where the complex breaks and the target molecule is reinserted into the liquid. We are mainly interested in very low amounts of target molecules, because the noise increases in this regime which is relevant for the detection limit. In this case, the mass-action approach is not valid, and we describe the dynamics in terms of probabilities for the respective reactions to take place.
The probability that an association event occurs for a specific target molecule within a time intervalτ is given by where r a is the hybridization constant, P is the total number of receptor molecules, and PT is the number of probe target complexes. Accordingly, the probability that a dissociation event occurs for a specific PT-complex during an intervalτ is given by where r d is the dissociation constant. This description is the discretization of a chemical Langevin-type equation [ Here, the dB i are independent Wiener processes and T denotes a free near-surface target molecule. In particular, the second and the third term in equation (3a) are responsible for the noise in the surface processes. Ignoring them leads to the corresponding deterministic rate equation. The second term represents the noise arising in the hybridization process, while the last term describes the one due to the desorption process. The prescribed initial conditions here indicate that there are no PT-complexes at the surface at the beginning of the simulation.
It is important to mention that not only is the surface of the nanowire functionalized with the receptor molecules, but so is usually the whole surface of the substrate. Therefore, the whole device surface needs to be considered for the assocation/dissociation dynamics, since off-sensor association influences the target-molecule concentration in the liquid and therefore the adsorption at the sensor.

Coupling to diffusion
In order to properly simulate the surface dynamics, we also need to take into account that the surface reactions are linked to the transport of the molecules in the liquid. Here we use the diffusion equation where D denotes the diffusion constant and ν the outward pointing normal vector to the domainΩ. In this work, the value of the diffusion constantD is derived by assuming a stiff rod-like geometry for the DNA oligomer. Then, we use the equation to findD, where the constants are explained in table 1. There are also other approaches to model the diffusion constant of DNA oligomers, but their use only leads to a small quantitative change in the final simulation results (see [10] for a discussion).
The right-hand side of the Neumann boundary conditions (4b) is defined such that the flux at the part of the boundary representing the device surface is exactly the net change in PT-complexes and is zero everywhere else. The net flux is described by equation (3), which couples the surface reactions to diffusion.
Note that the diffusion equation considers target-molecule concentrations rather than single particles. The transition to single particles is achieved via the random-walk based description of the diffusion.
Different initial conditions T 0 r can be used in (4c) and are discussed in the results section of this work.

The algorithm
In contrast to [10], we consider a two-dimensional domain here (see figure 1). This allows the implementation of more complicated, but also more realistic geometries in order to identify regions with special properties within the system. In fact, since the device is uniform along the nanowire, the investigation of the processes in two-dimensional cross sections is sufficient to cover all the effects taking place, so that three-dimensional simulations are not needed for elongated devices. The Brownian-motion based approach used here has several advantages compared to the box-based approach used in [10]. On the one hand, modeling the diffusion of the particles by a random-walk is much easier to implement than considering transitions between boxes in two or even three dimensions, since keeping track between admissible transitions is very awkward. On the other hand, all particles are moved at once, which allows an increase of the time step and hence yields less computational cost. However, the memory requirements are larger in this case, which is due to the necessity of keeping track of the position of every particle in the domain.
The time step, which is chosen to be in the range of a few microseconds, is small enough such that the particles only move a few nanometers in every time interval. Therefore, the use of the Euler-Maruyama scheme is sufficient for the simulation of the movement of the target molecules.
The algorithm, which is an extension of the standard random walk [30], is presented in figure 2. The definition of molecules located near the surface in Step3 is indicated by the dashed line in figure 1, which is of course not true to scale. We chose a value of 25 nm for the distance to the surface, since at this distance, the two ends of the DNA oligomers can contact and initiate the hybridization process.
To investigate the surface dynamics in different parts of the device, we furthermore partitioned the surface parts of 25 nm length and calculated the dynamics separately for each part (see the dotted lines in figure 1 and also figure 3 below) as well as for the whole nanowire. Running the algorithm once yields one realization of the stochastic process. To obtain the statistical quantities we are interested in, many simulations must be performed. Fortunately, since every run is independent of all the others, the work can be parallelized on many computing cores in a straightforward manner.
The algorithm presented in this section was implemented within the Julia environment [31,32].
The signal-to-noise ratio (SNR) is an important sensor characteristic and is shown in the results below. It is defined as the ratio

Results
The simulation results and their implications for sensor design are discussed in the following. In this work, we consider two topics. First, we investigate the hybridization behavior at different locations of the device surface, where we use a fixed receptor density and a fixed number of molecules, but vary the width of the simulation domain. Second, we investigate the total sensor behavior for different receptor densities, but fixed geometries. The numerical values for the reaction constants are taken from [10], while all the other parameters except those for the geometry of the system are taken from [28], where DNA hybridization was investigated experimentally. We chose a square-shaped cross-section of the nanowire for the simulations, because on the one hand it is a good approximation of the actual sensor shape, and on the other hand it allows a simple implementation of the partitioned surface speeding up simulations. However, other shapes like trapezoids or semicircular cross-sections can also be simulated using the algorithm presented here. The other parameters used are listed and explained in tables 1 and 2.

Partitioning the surface
Due to the geometry, the surface coverage of PT-complexes varies depending on location. We investigate the hybridization-dissociation processes at different locations by partitioning the surface: It turns out that we can distinguish three different regions, which will be called the edge, middle, and corner regions. Their locations are shown in figure 3. As will be seen throughout the remainder of the results section, the behavior varies significantly at these locations, but is similar for the same type of location for symmetric initial conditions.
In the following, we will consider the expected value, the variance, and the SNR both in equilibrium and as functions of time. To obtain the statistics for the respective quantities, we sampled from 2000 realizations for each specific setup.
In this set of simulations, we keep the number of target molecules fixed at 100 molecules and the geometry of the nanowire fixed at a square with a side length of100 nm. We only vary the width of the whole domain between 150 nm and 5 μm, yielding target-molecule concentrations between 2.7 and 92 μM. The receptor density used here is 3×10 12 cm −2 . Regarding the initial condition, we distributed the particles uniformly in the upper half of the simulation domain.

Signal
The equilibrium values of the signal PT ( )  are shown in figure 4(a). The values increase for increasing target-molecule concentrations, but are always highest in the edge regions and lowest in the corner regions. This is due to the fact that the volume from which molecules can bind to the surface is larger for the edge regions compared to the other parts, while the situation is vice versa for the corner regions.
Note that the average surface-charge density is very close to the surface-charge density in the middle region as is indicated by the black line in figure 4(a).
Next, we consider the time evolution of the hybridization at the highest target-molecule concentration investigated, which was 92 μM in this set of simulations. This particular value for the concentration stem from prescribing the total number of target molecules and the width of the domain. This setup yields a remarkable result, which is shown in figure 4(b). The number of probe-target complexes in the edge regions achieves its maximum after approximately 20 s; afterwards the number slightly decreases to reach its equilibrium value. However, this effect is not reflected in the overall signal, as is indicated by the black curve.

Noise
The results for the noise show a more complicated behavior than the results for the signal. The equilibrium values for the noise in figure 5(a) show the behavior already encountered in the simulations of flat surfaces [10], where the maximum is seen at medium target concentrations. The interesting difference here is that the position of the maximum depends on the type of surface region.
Another remarkable result is the fact the surface location (edge, middle or corner) with the highest noise in equilibrium also depends on the target concentration. For low densities, Table 2. Reaction parameters with respect to the receptor density at the surface. Data from [10].  the highest noise is encountered at the edges, while it is exactly the opposite when shifting to higher target densities. This inversion effect takes place around a target-molecule concentration of 30 M m and can be seen in figure 5(a). Finally, unlike for the signal, the overall noise is not similar to the middle corner behavior. In fact, for higher target concentrations, the decrease in the equilibrium noise is much more pronounced than for each of the single regions, while the maximum is reached at target-molecule densities similar to the edge region. Note that in this image, the overall noise is shown and not an average per area as for the signal.
Considering the time evolution again for 92 μM target concentration (see figure 5(b)), there is a maximum after a few seconds at all three locations, which is sharp in the edge and middle regions. At the edges, there is even a second oscillation resulting in a minimum after a few more seconds. The significant maximum is also observable for the overall noise (black line). This means that for higher target concentrations, there is a waiting time until the noise reaches its optimal value.

Signal-to-noise ratio
Despite the inversion effect in the noise, there is no intersection of the lines for the three regions in the SNR shown in figure 6. The SNR appears to increase according to a power law at all three locations.

Asymmetrical initial conditions
In order to investigate how the sensor behavior depends on different initial conditions, we considered asymmetrical initial conditions, which model realistic devices where the target molecules approach the nanowire transducer from one side.
Here, we initially put all particles in a small portion of the domain at the very right as indicated in figure 7(a) and simulated the time evolution of the system. As can be seen from figure 7(b), there is only a small difference between the left and the right half of the sensor, which even vanishes after tenths of a second. This clearly shows that the influence of the asymmetry in the initial conditions is negligible under the investigated conditions.

Different receptor densities
Now, we investigate the hybridization behavior for different receptor densities. In contrast to the previous experiments, we fixed the domain size here at 150×2000 nm, and varied the number of target molecules between1 and 150. The size of the nanowire was left unchanged. Again, we considered the signal, the noise, and the SNR.  Considering the binding efficiency in figure 8(a), the fact that it approaches 100% with an increasing number of target molecules is not surprising. As already encountered in [10], different numerical values for different receptor densities are observed. However, the absolute density of complexes in equilibrium is independent of the receptor density for low concentrations (see figure 8(b)), and a difference in the equilibrium value is only obtained at higher concentrations.
The standard deviation, which represents the noise in the system, is shown in figure 9(a). Again, the noise is almost equal for all investigated receptor densities at low target concentrations. The characteristic behavior that has already been discussed above is observed at concentrations higher than 10 μM. It is also remarkable that the position of the maximum in the noise also depends on the target-molecule density, which was not the case in the one-dimensional simulations.
The SNR shows a more complicated pattern in this set of simulations. For higher receptor densities, the SNR is similar to the already investigated 2D case for different domain widths and for the 1D simulations. However, for lower densities, the increase with increasing target-molecule concentration becomes larger between 10 and 100 μM. This behavior is shown in figure 9(b).
It is important to note here that the best SNR values are obtained for the lowest receptor density. This is somewhat counter-intuitive, and also contradicts the naive notion that a higher receptor density yields a better response. The simulations show that this is not necessarily the case, since in this case the low noise is responsible for a better distinctness of the signal from the noise.

Discussion
The numerical experiments show that the biomolecule detection depends on the geometry of the nanowire sensor, especially on the number and size of edge and corner regions. This dependence is not only seen during the evolution of the signal, but also in equilibrium. In particular, the investigation of a partitioned surface shows that in regions with a larger amount of liquid nearby (as at the edges) adsorption is much more likely and therefore yields equilibria similar to higher target-molecule concentrations. In corner regions, this behavior is just the other way round. This qualitative behavior can also be expected for other shapes like trapezoids or with semicircular cross-sections. In order to put these findings to use in experimental setups, it is straightforward to aim for sensor geometries where the accessible liquid volume is increased for most of the nanowire surface. Of course, a semicircular geometry seems to be a very favorable candidate in this sense, and trapezoid cross sections seem to be a good compromise between the optimal shape and manufacturability.
While the occurence of some hindrance effects in surface adsorption for very large receptor densities leading to lower signals has already been known [28], this study shows that at higher target-molecule densities the noise levels strongly increase with increasing receptor densities. However, since the quality of the signal decreases with increasing noise levels, the findings indicate that one should actually aim at having low receptor densities at the nanowire surface in order to obtain a high SNR whenever the target-molecule concentration is high.
Altogether, the numerical results indicate that it is not sufficient to only consider the average probe-target density at the surface. Rather, it is important to also take into account the fluctuations of this quantity as well as steric effects. Maximization of the SNR, and not only the signal, is important since a large noise level may prohibit the proper detection of the signal. The optimal number of receptor molecules then depends on the target-molecule concentration.

Conclusions
In this work, we developed a 2D algorithm and simulator for a coupled reaction-diffusion system representing the transport and association-dissociation of target molecules at nanowire biosensors. A random-walk based approach was used for the diffusion of the target molecules and coupled with a stochastic-simulation algorithm for the association and dissociation of the target molecules with their receptors at the surface. This approach makes it possible to investigate noise and fluctuations.
Because of the geometry of the nanowire transducer, the behavior of the association-dissociation process is not uniform across the sensor surface. In fact, the data at nanowire edges shows higher PT-complex densities, while corner areas show lower densities. This is due to the size of the liquid area in which hybridization is possible.
The consideration of asymmetrical initial conditions, which model the initial phase of detection when the target molecules approach the nanowire from one side, revealed that there is only a very little effect on the system behavior. This is due to the fast diffusion at these length scales compared to the association probability.
The dependence of the hybridization behavior on the receptor density is also a crucial factor in optimization of sensors of this kind. The investigation of this topic showed that there is almost no difference in a low target-molecule regime, while at higher concentrations, it appears to be better to employ lower receptor densities at the surface.