Modeling convection-diffusion-reaction systems for microfluidic molecular communications with surface-based receivers in Internet of Bio-Nano Things

We consider a microfluidic molecular communication (MC) system, where the concentration-encoded molecular messages are transported via fluid flow-induced convection and diffusion, and detected by a surface-based MC receiver with ligand receptors placed at the bottom of the microfluidic channel. The overall system is a convection-diffusion-reaction system that can only be solved by numerical methods, e.g., finite element analysis (FEA). However, analytical models are key for the information and communication technology (ICT), as they enable an optimisation framework to develop advanced communication techniques, such as optimum detection methods and reliable transmission schemes. In this direction, we develop an analytical model to approximate the expected time course of bound receptor concentration, i.e., the received signal used to decode the transmitted messages. The model obviates the need for computationally expensive numerical methods by capturing the nonlinearities caused by laminar flow resulting in parabolic velocity profile, and finite number of ligand receptors leading to receiver saturation. The model also captures the effects of reactive surface depletion layer resulting from the mass transport limitations and moving reaction boundary originated from the passage of finite-duration molecular concentration pulse over the receiver surface. Based on the proposed model, we derive closed form analytical expressions that approximate the received pulse width, pulse delay and pulse amplitude, which can be used to optimize the system from an ICT perspective. We evaluate the accuracy of the proposed model by comparing model-based analytical results to the numerical results obtained by solving the exact system model with COMSOL Multiphysics.

Introduction their sensing operations through molecular signals in blood flow, where convection and diffusion act simultaneously on the transport of molecules [1]. Furthermore, it can also find use in microfluidic networked lab-on-a-chip devices, which is an emerging technology to diversify the point-of-care medical applications and increase the efficiency of on-chip diagnostics [18]. Moreover, imitating the transport of molecules with convection and diffusion in confined geometries like vascular and neuro-synaptic channels, similar microfluidic configurations can find application in organs-on-chips and artificial synapses relying on molecular information and communication technologies [3,[19][20][21].
The overall process, which covers the release of ligands by the transmitter in the form of a finite-duration concentration pulse, molecular transport in laminar flow through microfluidic channel, and the molecular detection by the surface receiver equipped with finite number of ligand receptors is a highly nonlinear and time-varying process, and does not yield an analytical solution for the ligand and bound receptor concentration fields. Therefore, it necessitates the application of computationally-expensive numerical methods, like finite element analysis (FEA). In this study, we develop an end-to-end analytical model that can capture the expected time course of the received signal in terms of number of bound receptors. The model is based on quasi-steady state two-compartment model, which we tailor to the time-varying characteristics of the microfluidic MC system. The resulting model captures the effect of channel and receiver geometry, and system parameters regarding the fluid flow and ligand-receptor reaction. It accounts for the nonlinearities caused by laminar flow resulting in parabolic velocity profile and finite number of receptors resulting in saturation of the receiver. The effect of interplay between reaction and transport rates, which can lead to depletion layer over the receiver surface is also covered. Based on the developed model, we derive approximate analytical expressions for the received pulse delay, pulse amplitude and pulse width to help characterize and optimize the system from communication theoretical perspective. The analytical results are compared to numerical solutions obtained with the exact system model by using COMSOL Multiphysics, which is a finite element analysis simulation software.
The remainder of the paper is organized as follows. In the Methods section, we present the exact model of the communication system with nonlinear partial differential equations, and we introduce our approximate model and derive the communication theoretical metrics for the received signal. A comparative analysis of the analytical and numerical results is provided in the Results and Discussions section. Finally, the concluding remarks are given at the end of the paper.

Communication system model
In this section, we represent the end-to-end model of the microfluidic communication channel as a system of partial differential equations. We utilize a two-dimensional (2-d) model considering a microfluidic channel with rectangular cross section as shown in Fig 1A. 2-d models are proved effective in modeling the molecular transport, especially when there is an obvious interplay between convection, diffusion and surface reaction, as the uniformity of the molecular concentration along y-direction is disturbed above the reactive surface [22]. On the other hand, for cases, where there is no reactive surface, 1-d models can succesfully capture the effect of convection and diffusion [23].
Using a similar notation to that of [22], we define three orthogonal domains: (i) bulk domain O b , where the convection and diffusion of ligands occur, (ii) reacting surface domain O rx denoting the biorecognition layer of the receiver, where the ligand-receptor reaction occurs, and (iii) non-reacting surface domains O nr defining the walls of the microfluidic channel excluding the receiver surface, as demonstrated in Fig 1B. We also define two dependent variables c = c(x, y, t) denoting the ligand concentration in space and time domain, and R = R(x, t) denoting the bound receptor concentration at the receiver surface.
Hence, in a 2-d convection-diffusion system the propagation of ligands is governed by where r 2 ¼ @ 2 @x 2 þ @ 2 @y 2 is the two dimensional Laplace operator, u x (y) is the flow velocity as a function of distance to the surrounding walls. We assume fully developed laminar flow in the microfluidic channel yielding a parabolic flow velocity profile, i.e., where u is the maximum flow velocity. D is the effective diffusion coefficient taking into account the effect of Taylor-Aris type dispersion [14]. For a channel with rectangular crosssection, it is given by where the intrinsic diffusion coefficient is denoted by D 0 [14]. Here, h ch and w ch denote the height and width of the microfluidic channel, respectively. No-flux boundary condition is assigned to the non-reacting walls of the channel, i.e., @c @y ¼ 0; y 2 Ω nr : ð4Þ Flux condition at the reactive boundary, where ligand-receptor reaction occurs is given by where R is the bound receptor concentration, and R denotes the reactive flux. We define the inlet and outlet boundary conditions as where c in (t) is the transmitted signal. Assuming no surface diffusion for receptors at the receiver surface, bound receptor concentration as a function of time can be written as Given the finite receptor concentration at the receiver, ligand-receptor binding reaction can be described by the first-order Langmuir kinetics giving the reactive flux: Lastly, the initial conditions for the system are defined as cðx; y; 0Þ ¼ 0; ðx; yÞ 2 Ω b : ð11Þ The exact system model presented above is not analytically tractable and necessitates numerical methods to compute the ligand and bound receptor concentration.

Proposed model
We develop an analytical model that approximates the time course of the mean number of bound receptors on the receiver surface for a rectangular ligand concentration pulse transmitted at the channel inlet, as shown in Fig 1C. Prior to modeling, it is worth elaborating briefly on the impact of conditions resulting from the competition between ligand transport and ligand-receptor reaction. In convection-diffusion-reaction systems, if the convective/diffusive transport in the channel supply ligands much more quickly than the receptors can bind them, then the system becomes reaction-limited, implying that transport dynamics has negligible effect on the resulting waveform for the bound receptor concentration. In such cases, ligand concentration near the receiver surface can be assumed equal to the concentration of ligands supplied at the channel inlet, and the well-mixed condition can be assumed for modeling the ligand-receptor reaction [22].
In practical systems, however, the reaction-limited condition often does not hold, because the concentration of supplied ligands is not continuous, giving rise to a concentration gradient between the channel inlet and the reactive surface [24,25]. Furthermore, in microfluidic channels, the fluid flow is usually laminar leading to a parabolic flow velocity profile above the reactive surface, implying a concentration gradient from the center of the channel toward the reactive surface [22]. Moreover, when the system is utilized for communication purposes, reaction rates at the receiver surface should be kept large enough compared to the transport rate in order not to cause intersymbol interference (ISI). Because of these reasons, we assume that the system would be operating either in transport-limited regime, where the reaction rates are much larger than the ligand transport rate, or under partial mass transport limitations, caused by transport and reaction rates comparable to each other. In such cases, mass transport limitations could have substantial impact on the time course of bound receptor concentration, and well-mixed condition for the surface reactions often does not hold.
A compartmental approach is developed in [26] to model the ligand-receptor kinetics affected by mass transport limitation. It is based on dividing the space domain into two compartments, in each of which the ligand concentration can be assumed steady. This steady-state assumption proved to be effective in describing the isolated association and dissociation phases (two-stage process). The two-compartment model is widely employed in BIAcore analyses to determine the affinity of various ligand-receptor pairs [24].
We make use of the two-compartment model by tailoring it to the peculiarities of the system under investigation. We will demonstrate that this simple yet effective strategy can be used to reveal the characteristics of the microfluidic MC channel. Two-Compartment model. Here, we review the basics of two-compartment model that we modify to propose an approximate model for the microfluidic communication channel.
In the two-compartment model, ligand concentration of the first compartment is assumed to be equal to the concentration at the channel inlet c = c 0 [26]. The concentration in the second compartment covering the receptors c| y = 0 could be different than the bulk concentration due to the binding reactions occurring at the receiver surface. Concentration gradient between two compartments results in a flux of ligands expressed by where k T is the mass transport parameter, which is a function of the channel and receiver geometries, and the system parameters regarding diffusion and convection of ligands [25], i.e., where A rx = l rx w rx is the receiver surface area with the receiver width w rx = w ch . Here, C T is defined as where d rx is the minimum distance of the channel inlet to the receiver, l rx is the length of receiver along the x-axis, and F = h ch × w ch × u is the maximum flow rate. The steady-state solution of this two-compartment model for ligand concentration near the receiver surface is obtained by equating the ligand flux to the surface J given in Eq (12) to the reactive flux R given by Eq (9), i.e., where the position dependent surface receptor concentrations in Eq (9) are integrated over the receiver surface to obtain the number of receptors in units of mol, i.e., Substituting this expression into Eq (9) yields a nonlinear differential equation representing the evolution of the number of bound receptors, i.e., The solution of Eq (16) for the association phase is then obtained by setting the initial conditions as c(x, y, 0) = c 0 and R(0) = 0 [25,27], i.e., where The solution of Eq (16) for the dissociation phase can be found by setting the initial conditions as c(x, y, 0) = 0 and N R (0) = N R,0 , with N R,0 being the number of bound receptors at the start of the dissociation phase [25,27], i.e., Proposed model. We propose an approximate model built upon the two-compartment model. In two-compartment model, for each of the phases, e.g., association and dissociation phases, the bulk concentration in the first compartment is assumed constant and equal to its initial value. However, in our system, the finite duration input to the system results in pluglike passage of the transmitted ligands over the receiver surface as shown in Fig 1C, making it impractical to directly apply the two-compartment model.
Our strategy is to consider the whole process again as a two-phase process consisting of association and dissociation phases, and then, set the start time instants of association and dissociation phases to reflect the plug-like passage. We also define an effective plug length and effective plug concentration in order to elaborate the system within the framework of twocompartment model.
We start by modeling the propagation in the microfluidic convection-diffusion channel. Once we are done with the propagation, we incorporate the receiver reactions into the model with a modified two-compartment model. Distance between transmitter and receiver d is taken as the distance of the center of the receiver to the entrance of the channel where transmission occurs: When the transmitter sends an impulse signal in the form of surface concentration, i.e., N/ (A ch N A ) = 1(mol/m 2 ), neglecting the reaction at the receiver for the moment, the ligand concentration anywhere in the channel at any time can be written as [14] The concentration function is Gaussian with variance s 2 imp ¼ 2Dt, which is varying with time, making the subsequent calculations analytically untractable. Therefore, during the passage of the resulting ligand plug, we neglect the dispersion and assume that the variance is constant at the receiver location, making the plug a Gaussian function with variance s 2 R ¼ 2Dd=u, moving at the flow velocity. Hence, end-to-end impulse response of the communication channel from the transmitter to the receiver location becomes In this study, we assume a more practical signal, i.e., finite-length rectangular pulse, as the input. Accordingly, transmitter releases a total of N molecules at a constant rate μ T for a specified pulse duration T p uniformly through the entire surface of the channel inlet; thus, the transmitted signal can be written as where rect(t) = 1 for −0.5 < t < 0.5 is the rectangular function, and μ T = N/T p is the transmission rate in molecules/s. The convection-diffusion channel excluding the ligand-receptor reactions is a linear timeinvariant (LTI) system; thus, the response to a rectangular pulse can be found via convolution, i.e., where N A is the Avogadro's number.
We define the plug delay as the time of the passage of the peak ligand concentration through the center of the receiver, i.e., To utilize the two-compartment model, we approximate the Gaussian ligand concentration pulse passing over the receiver as a finite duration pulse with a pulse length set to cover approximately 95-96% of the ligands, as done in [28]. This effective plug length is computed as where s R ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2Dd=u p . Then the effective plug passage is given as Every point along the receiver surface is assumed to be exposed to a stationary concentration between a time window, whose end-points are marked by association and dissociation times. Taking the plug delay as the central point in time, we set the association time as and the dissociation time is set to Then, the effective plug concentration is calculated as the time-average of the concentration passing through the receiver surface between t a and t d , i.e., where Q(t) is given by Given the definitions regarding the effective plug concept, the time course of mean number of bound receptors for a system where transmission starts at t = 0 obtained by the two-compartment model can be given by where Θ[.] denotes the Heaviside step function, and the modified parameters are given as where the transport parameter is set to k Ã T ¼ k T Â k to be optimized for our system. Received pulse characteristics. In this section, based on the proposed model, we attempt to characterize the microfluidic MC channel by providing analytical expressions for three metrics: pulse delay, pulse amplitude, pulse width. A similar approach has been previously taken for diffusion-based MC channel in [29]. In the next section, we will compare our results to that obtained through finite element analysis numerical solutions in COMSOL Multiphysics for different system settings.
Pulse delay, different from the plug delay given in Eq (28), is defined as the time instant, at which the number of bound receptors reaches its peak value. In our model, the number of bound receptors given by Eq (35) is monotonically increasing for t t d , and monotonically decreasing for t > t d . Therefore, the pulse delay can be written as Then, the pulse amplitude, as the peak value of the received signal, is given by Another important metric is the pulse width. As done in [29] for MC, we define it as the time interval, at which the pulse magnitude is greater than half of its peak value. Pulse width has implications for the achievable bandwidth as it determines the extent of intersymbol interference (ISI). Based on our model, pulse width is calculated as where N Ã R;eq , α Ã , β Ã , and γ Ã are given in Eqs (36)-(39). The above analytical expressions for the received pulse delay, pulse amplitude and pulse width are of paramount importance for MC engineering, as they enable an optimization framework that can be utilized to optimize the overall system from ICT perspective, and help develop advanced communication schemes, such as optimum detection methods and modulation techniques.

Results and discussions
In this section, we first estimate the transport parameter k Ã T ¼ k Â k T , and then, evaluate the acuracy of the proposed analytical model by comparing the calculations under different conditions to the results obtained via finite element analyses in COMSOL Multiphysics. The default values for the system parameters used in the analyses are listed in Table 1.
To find the optimal value of the free parameter k, we perform nonlinear least square estimation in MATLAB using Levenberg-Marquardt optimization algorithm. The curve fitting is conducted on the numerical results obtained via COMSOL for 35 different scenarios, in each of which we vary only one parameter from its default setting at a time. The results of the optimization with the corresponding system parameters are presented in Fig 2. The mean and standard deviation of the obtained data are 0.5138 and 0.3616, respectively. Setting k to its mean, we rewrite the transport parameter optimized for the communication scenario as Using the optimized transport parameter, we compare the time course of normalized number of bound receptors obtained via Eq (35) under different scenarios to the numerical results of COMSOL experiments. The results are presented in Fig 3. Clearly, our analytical model approximates quite well the numerical solution, justifying the accuracy of the model.
In the second part, we analyze the capability of our model to reflect the characteristics of the received signal under various conditions, this time, using the metrics defined in the previous section. Accordingly, the pulse delay obtained via Eq (40) is compared to the numerically computed results in Fig 4. As is seen, the simple expression given in Eq (40) is quite accurate in approximating the exact results and following the trends with varying parameter values. It is worth noting that the received pulse delay is not affected by the molecular transmission rate, surface receptor concentration, and the binding/unbinding rates of the ligand-receptor pairs. As expected, it is most affected by the minimum TX-RX distance d rx and flow velocity u.
The analyses are repeated for the normalized pulse amplitude, which is the ratio of the peak number of bound receptors to the total number of receptors. Expected pulse amplitude is quite important for the communication system design, as it determines the received SNR. The results obtained via Eq (41) are compared to numerical results in Fig 5. The parameters range between values corresponding to sparsely occupied and saturated receiver. The results reveal that almost all parameters have substantial effect on the pulse amplitude, and the exact numerical results and the corresponding trends are well approximated by the proposed analytical model, no matter whether the receiver is saturated or sparsely occupied. One interesting result observed by both numerical and analytical calculations is obtained for varying flow velocity. As shown in Fig 4G, given the fixed molecular transmission rate, the pulse amplitude is decreasing when the flow velocity is higher or lower than a certain optimal value. This is because for low velocities, the degree of attenuation through dispersion until the plug arrives at the receiver location increases, and for high velocities, the time duration, at which the ligand plug and receiver are in contact decreases, leading to a lower number of bound receptors. The same trend was also observed and discussed in detail in [30][31][32].
The third set of analyses are conducted for the pulse width, which is an important parameter since it is directly linked to the extent of ISI and achievable bandwidth. Given a fixed symbol duration, ISI increases with pulse width. The analytical results compared to the numerical calculations are presented in Fig 6. The proposed model quite well reflects the characteristic trends observed under varying conditions. The results reveal that the unbinding rate of the employed ligand-receptor pair is the most critical parameter to adjust the received signal pulse width.
Our model, which is quite simple, practical, not requiring the computationally expensive numerical methods, is able to accurately capture the design tradeoffs, and could be used to design efficient and reliable microfluidic MC system before the final implementation.
Finally, we would like to discuss about an alternative communication scheme that can be targeted to improve the performance of the microfluidic MC with surface-based receivers. The performance could be substantially improved if multiple types of ligand-receptor pairs with different binding characteristics, are employed in the communication system [4]. This way, the transmitter can employ molecular division multiplexing, similar to the code division multiplexing in conventional EM communication systems, to boost the communication rate by simultaneously transfering multiple messages in the same channel without causing significant interference. This can be reliazed in two ways. (i) Employing multiple receiver antennas, e.g., surface-based biosensors, with each having different kind of receptors that bear affinity to different kinds of ligands. In this scheme, the receiver can process the output of each antenna separately with minimum interference. (ii) Employing multiple types of receptors, at a single receiver antenna, corresponding to different kinds of ligands. In this scheme, the output of different kinds of receptors are superposed at the receiver output; therefore, this scheme requires employing more complex signal processing to discriminate the contributions of different messages conveyed through different ligands. Fortunately, there are some studies proposing frequency domain detection techniques to exploit the diversity in the critical frequency of ligandreceptor binding noise.

Conclusion
In this paper, we develop an analytical model for microfluidic MC systems with surface-based receivers equipped with ligand receptors to approximate the time course of the number of bound receptors. The model is based on the two-compartment model of convection-diffusionreaction systems, which is tailored to the peculiarities of the communication system. The comparison of the analytical model-based results with the numerical results obtained by solving the exact model in COMSOL proves the accuracy of the developed model, which is quite successful in capturing the nonlinearities of the system. We also provide analytical expression for received pulse width, pulse amplitude and pulse delay to provide a framework that would help optimize the microfluidic systems from communication theoretical perspective.