The application of communication engineering principles to biotechnological problems has prompted a large interest and has motivated the investigation of the transport of molecules as carriers of information, see  for a review of recent literature.
The transport of molecules, or more generally particles, differs fundamentally from the propagation of electromagnetic waves and therefore requires new theoretical models for system design and analysis . Though negligible at macroscale, the transport by diffusion becomes relevant at nanoscale. However, diffusion can be slow and thus relying only on diffusion might be limited to small distances such as in intracellular communication. In fact, in biotechnology, active transport mechanisms induced by an external force, fluid flow, or a combination thereof are well established . In particular, magnetic forces are attractive because the magnetization of biological particles is often negligible and special purpose magnetic nanoparticles can be engineered and tailored to an application by adapting their size and composition . In this way, magnetic nanoparticles can be selectively detected and externally steered in a preferred direction by a magnetic field. Furthermore, magnetic nanoparticles can be made biocompatible by suitable coating to avoid reactions with undesired reactants and allow for specific binding to others .
Motivated by the concept of using magnetic nanoparticles for molecular communication proposed in [4, 5] and the practicality of studying blood vessels [6, 7, 8], here we consider particle transport subject to horizontal fluid flow, diffusion, and a vertical magnetic force in a cylinder model of a blood vessel. For simplicity, we assume a uniform axial flow, which is a typical model known as plug flow . Although blood vessels may have a complicated shape, modeling them as cylinders is well-accepted . To study the distribution of particles along a segment of the vessel wall, a line receiver model is considered.
When neglecting diffusion, the particle trajectory in a cylinder due to fluid flow and magnetic force can be accurately modeled even for more complicated flow patterns . However, modeling diffusion, flow, and particle drift by magnetic force in a cylinder is analytically demanding due to the non-symmetry of the problem. Therefore, usually finite element numerical simulation is employed to study this type of system, e.g., in drug-targeting [11, 12]. However, numerical solutions lead to algorithms with typically long run-times and low flexibility due to changing channel parameters.
A different modeling technique is the application of functional transformations, which are based on the modal expansion of an initial boundary value problem [13, 14]. This technique has been already applied successfully, e.g. in the field of sound synthesis and circuit modeling, see [15, 16] for further references. This kind of modeling technique offers several benefits: It leads to semi-analytical models in terms of a state-space description which provide insights into the diffusion process. Furthermore, these models allow the adaptation of their boundary behavior using feedback techniques [15, 17].
This paper proposes a semi-analytical model for a diffusive cylindrical molecular communication channel in the presence of a magnetic force and flow. The model leads to a low run-time simulation algorithm, which provides several benefits compared to purely numerical models. The derivation of the model follows the mathematically more detailed description in , where particle diffusion under the influence of a magnetic force in a circular disk is simulated. Our results might be useful for studying molecular communication in, e.g. blood vessels, flow reactors or microfluidic channels. It is also applicable to the study of diffusion with horizontal slurry flow in tubes where particles sediment due to gravity [19, Chapter 13].
The paper is structured as follows: In Section II, a physical description of the scenario under study is presented. In Section III, a semi-analytical model for the diffusive cylindrical channel is derived and a simulation algorithm is proposed. The validity of the model is shown in Section IV by simulation, and the algorithm is analyzed and its benefits are discussed. Finally, Section V concludes the paper and proposes topics for future work.
Ii Problem Description
Fig. 1 shows a simplified cylindrical model of a blood vessel. Information carrying nano particles are emitted into the channel of radius by a point source at . They are diffusing with a constant diffusion coefficient and are dragged by a uniform flow in -direction with velocity . The particle movement is controlled by a vertical magnetic force causing a constant uniform drift velocity , see  for details. For modeling the reception, the particle concentration is desired at a point on the vessel wall where a receiver nanomachine is mounted. As an example generalization, for probing the dynamics of the particle concentration over an extended axial distance, a line receiver of length is considered. The particle concentration over time in the receiver is denoted by (the -index indicates the position in
-direction), the particle flux is denoted by vector.
Before the physical description of the scenario in Fig. 1
in terms of partial differential equations (PDEs) is established, the variables in the system are normalized. Introducing a reference timeand a reference length in terms of the radius of the cylinder and the diffusion coefficient leads to the following normalization
In the following, concentrations and fluxes in the cylindrical geometry are formulated in polar coordinates, using the mapping and , particularly . With this normalization and substitution a normalized dimensionless physical model is considered in the following sections.
Ii-B Physical Model
Due to the geometry of the problem, the average concentration in a line receiver, , following a point release, can be expressed as
Here, is the accumulated concentration in the line segment for a one-dimensional drift-diffusion process in -direction. It is obtained by integration [1, Eq. (23)]
Furthermore, is the concentration of a diffusion process in a circular cross-section of a cylinder. Its calculation is the main task in the following sections. We note that a similar separation as in (3) also holds for the flux .
The particle concentration in the circular disk , under the influence of the vertical drift velocity is given by the following PDE
with the unit vector in -direction , the delta impulse and gradient grad and divergence div in cylindrical coordinates. The emission of particles is described by the initial condition (5d) and the zero flux at by the boundary condition (5c).
Ii-C Transformation of Variables
The presence of the drift term in (5b) highly scales the complexity of a direct solution of the system (5a), (5b). Therefore, the problem is transformed into an auxiliary problem with the auxiliary particle concentration by application of the substitution , see also 
The main task in the following sections is to derive a semi-analytical model for the auxiliary concentration with the non-trivial boundary condition (6c). In the end, this result is used to obtain a semi-analytical simulation model for the desired concentration (3) at the receiver.
Ii-D Vector Formulation
Before the transfer function model for the considered diffusion process is constructed in Section III, the auxiliary PDE in (6a), (6b), the boundary conditions (6c), and the initial conditions (6d) are rearranged into a unified vector form according to . The reformulation of the PDE in (6a), (6b) in terms of a matrix valued spatial differential operator defined on the circular disk with and leads to
with a mass matrix , a matrix of damping parameters
, and the identity matrix. The involved block matrices and operators are defined as
Vector contains the independent variables of the system: the auxiliary particle concentration and the flux in - and -direction
where denotes the transposed. According to the dimensions of and in (9), the matrices and vectors in (8) have sizes and , respectively. The initial conditions for the auxiliary particle concentration in (6d) are rearranged similar to the vector of variables in (9)
In addition to the desired boundary behavior (6c), a second boundary condition is defined
with a general boundary excitation . The second set of conditions (11) is used to design a general transfer function model for the diffusion process in Section III-A. The derived model serves as a blueprint for a diffusion process, where different kinds of boundary conditions can be incorporated. In Section III-B, the desired boundary behavior (6c) is realized by a modification of the blueprint model.
Ii-E Laplace Transformation
where is the frequency domain equivalent of the vector of variables in (9). The complex frequency is denoted by . The frequency domain representation of the PDE in (12) and the boundary conditions (6c), (11) are the basis for the subsequent transformations.
Iii Semi-Analytical Model
In this section, a transfer function model for the circular diffusion process is derived. First, a model for the process described by the PDE in (12) with the second set of boundary conditions in (11) is established using a modal expansion of the spatial differential operator
into eigenfunctions[13, 14]. Then, this model is formulated as a state-space description and is forced to fulfill the boundary conditions in (6c) by the design of a feedback loop .
Iii-a Transfer Function Model
To derive a transfer function for the given 2-D diffusion process (12), the spatial differential operator is expanded into eigenfunctions (see Section III-A3) based on the Sturm-Liouville theory .
Iii-A1 Sturm-Liouville Transformation
For the transformation of space variables, a Sturm-Liouville transformation is defined in terms of the integral transform 
with the eigenfunctions . The inverse transformation is a generalized Fourier series with the eigenfunctions and the scaling factor
where is a scalar representation of the vector of variables in (9) in the temporal and spatial transform domain. Index is the index of a spatial frequency variable for which (7) has nontrivial solutions [14, 15].
Iii-A2 Transform Domain Representation
with the transform domain representation of the initial conditions from (10) and the transformed boundary term given by
with the -th order Bessel functions of the first kind and is the frequency domain equivalent of the boundary excitation in (11). The values represent the spatial eigenfrequencies of the diffusion process. Solving (15) for the transformed vector of variables leads to a representation in terms of a multidimensional transfer function
Iii-A3 Eigenfunctions and Eigenfrequencies
The eigenfunctions for the forward and inverse transformation are derived by solving a dedicated eigenvalue problem . Its derivation is skipped here for brevity but it is described for similar problems in [21, 15]. The eigenfunctions in (13), (14) are obtained as
where . The eigenfrequencies of the diffusion process are calculated by the evaluation of boundary conditions similar to (11), i.e., from the real-valued zeros of .
Due to the bi-orthogonality of the eigenfunctions, the scaling factor is calculated by evaluation of the integral 
Iii-A4 State-Space Description
Grouping all elements of the matrices and vectors in (15) in the range of into vectors gives a state equation in the continuous frequency domain
with the matrices and vectors
The diagonal matrix contains the eigenvalues of the system on its main diagonal. The reformulation of the inverse transformation in (14) leads to an output equation for the state-space description,
Here, is the -th element of eigenfunction in (19). The complete state-space description based on state equation (21) and output equation (24) is shown in Fig. 2 (switch open). The matrices and vectors in (21), (24) are theoretically of infinite size (see range of in (14)), and therefore they have to be interpreted in an operator sense. For numerical evaluation, the sums in (14) have to be truncated so that the matrices become finite and computable.
This state-space model was obtained for the boundary conditions in (11) and includes the yet undetermined general boundary term . Setting this term to zero corresponds to the scenario of particle diffusion in a cylinder with zero flux at the boundary [22, p. 378]. In the next section, the boundary terms are used to realize the desired boundary behavior (6c).
Iii-B Feedback Loop
Iii-B1 Feedback Matrix
The state-space model in (21), (24) is designed with the general boundary term . Now, the excitations are replaced by the desired ones (6c). The goal is an expression for in terms of the state vector of the state-space description. A feedback matrix for the system is constructed according to the concepts described in . The starting point is the general transformed boundary term (17). Inserting the desired boundary conditions (6c) into the general term in (11) (in the frequency domain) leads to
with the coefficients and a function defined in terms of Kronecker delta functions and the imaginary unit
Sorting the sums in (28) into vectors leads to an expression for the transformed boundary term realizing the realistic boundary conditions in terms of the state vector
with the feedback matrix
The system shown in Fig. 2 with open switch realizes the diffusion process (6a), (6b) with boundary conditions (11). By incorporation of (30) and closing the switch, the system fulfills the desired boundary conditions (6c).
Iii-B2 Modified State-Space Description
The modification of the state matrix shifts the eigenvalues such that the realistic boundary conditions (6c) are fulfilled.
Iii-C Simulation Algorithm
To obtain a simulation algorithm for the desired concentrations and fluxes in (9), the modified state-space description has to be transformed into the discrete-time domain.
Iii-C1 Impulse Invariant Transformation
Sampling both equations in (34) with the sampling interval as and applying an impulse invariant transform (IIT) results in a synthesis algorithm in the discrete-time domain for the vector of variables
Iii-C2 Transformation of Variables
To obtain a simulation algorithm for the particle concentration in the receiver in Fig. 1, three steps are applied. First, the output equation of the state-space model is restricted to get the auxiliary concentration
Third, the concentration in the cylindrical slice has to be recombined with the discrete-time equivalent of the concentration in a line segment (4) to obtain a solution for the complete cylinder
The individual steps, which are performed in the implementation and the underlying equations, are shown in Table I. Steps 1-3 can be performed before the discrete-time simulation. In step 1 the sums in (14) are truncated to and . This leads to the total number of terms . Step 3 constructs a state-space model of the diffusive channel with the generalized boundary term . In step 4 the desired boundary behavior is chosen - here (6c) - and the feedback matrix is constructed in step 5. With that matrix, the modified state matrix is computed in step 6. In step 7, the auxiliary concentration is computed with (35), (36) with a for-loop over time. Finally, the auxiliary concentration is transformed into the receiver concentration in steps 8,9.
The algorithm assumes fixed parameters (, , etc.), but variable parameters can be easily incorporated and not all steps have to be repeated. Changing e.g. the velocity , the algorithm stays the same for steps 1-3, a change of the horizontal flow affects only step 9. Also, to incorporate time-varying parameters, e.g. for a time-varying velocity , steps 5,6 have to be moved inside the for-loop in step 7. This can affect the runtime of the algorithm.
|1. Choose index range|
|2. Compute zeros of|
|3. Compute state and output matrices||(19), (20), (23), (25)|
|4. Choose desired boundary behavior||(6c)|
|5. Compute feedback matrix||(30), (31), (32)|
|6. Compute modified state matrix||(33)|
|7. for k = 1:simulationDuration compute auxiliary concentration||(35), (36)|
|8. Transformation of variables||(37)|
|9. Compute receiver concentration||(38)|
Iv Verification and Analysis
In this section, simulation results for the semi-analytical model described in Section III are presented. The model outcome is compared to particle-based simulations . Furthermore, the proposed semi-analytical model is analyzed in view of runtime and accuracy. Also, the benefits of a semi-analytical model over purely numerical models are discussed.
Iv-a Methods and Parameters
The radius of the cylinder is . The diffusion is characterized by diffusion coefficient and horizontal uniform flow with (see ). All parameters are normalized according to Section II-A with reference length and time . The injection of particles into the system is performed at , , .
The simulations use a Matlab-Implementation of the semi-analytical model with a total number of terms (, ). The presented results were produced by the algorithm in Table I with a normalized sampling time .
For validation, we employ a particle-based simulation of the advection-diffusion transport where the positions of particles are updated and tracked in discrete normalized time steps of length following diffusion and drift, see e.g. [1, Eq. (1)]. Within each time step, if a particle crosses the channel boundary, it is reflected back into the channel. For measuring the concentration within a line receiver, we consider a normalized cuboid of size approximating the line receiver in Fig. 1
. An estimate of the accumulated concentration is obtained by counting the number of observed particles within the cuboid and dividing by the surface areaand the number of released particles . For , results are averaged times.
The parameters, which are varied in the simulations, are velocity and the receiver position with a normalized length of . A total of four receiver positions are used with , , each placed at .
Iv-B Particle-Based Verification
In Fig. 3, the average concentration along the receiver (see Sec. IV-A) placed near the boundary of the cylinder is presented. The average concentration following an instantaneous point release of magnetic nano-particles represents the impulse response in a molecular communication system. For such a scenario, the increased signal strength induced by the drift velocity caused by a magnetic force is evaluated and can be directly related to the reliability of transmission . The figure shows the result of the semi-analytical model according to (38) (solid and dashed lines) and the outcome of the particle-based simulation (circle markers) for different drift velocities . The semi-analytical model and the particle-based simulations are in excellent agreement. The normalized width of the received signals in Fig. 3 can be well approximated by which accounts for the length of the line receiver. Hence, the diffusion along the -axis does not significantly contribute to the width of the received signal. In contrast, the height of the received signal critically depends on , , and the receiver position and is governed by the combination of cross-sectional diffusion and the vertical drift .
Iv-C Algorithm Analysis
Fig. 4 shows two plots of the concentration in the receiver simulated with the proposed model (38) at the bottom of the cylinder. The left plot follows the blue line in Fig. 3 and shows the particle concentration at with a drift velocity of . For the numerical analysis of the proposed model, the range of is fixed to and the range of is varied as . The circle markers are the result of the particle-based simulation. The right plot in Fig. 4 shows the receiver concentration at time as a function of for decreasing number of terms. As already shown in Section IV-B for (), the particle-based simulation and the results of the proposed semi-analytical are in excellent agreement.
The runtime of the algorithm is to simulate a normalized duration of starting at . For the semi-analytical model still agrees very well with the results of the particle-based model. The runtime decreases to . A further reduction of the number of terms leads to visible deviations, but the runtime drops to , . For comparison, the runtime of the particle-based simulation is approx. to simulate a normalized duration . The results in Fig. 4 reveal a trade-off between runtime and accuracy. This shows the flexibility of the proposed model: If the exact concentrations are needed, the number of terms, , can be chosen large at the cost of runtime. On the other hand, to study the rough behavior of a channel for different parameter sets, the number of terms can be reduced to run many simulations in a short time. The plot on the right hand side of Fig. 4 shows that the accuracy also depends on the velocity and, for a given , decreases for increasing .
Iv-D Benefits of the Semi-Analytical Model
The proposed semi-analytical state-space model (see Section III) offers many advantages over purely numerical models.
In particular, the semi-analytical description of the diffusion process allows to compute the channel impulse response in an analytical manner. Hence, solutions in the form of (35) can be related to parameters relevant for communication system design (e.g. inter-symbol interference) more directly than solutions for purely numerical models. Additionally, the complexity of the simulation algorithm is flexible (varying ), and can be adjusted to the simulation purpose (see Section IV-C). Knowing the eigenvalues of the diffusion process provides insight into the behavior, e.g. damping or asymptotic behavior. As the output variable in (35) also contains the fluxes , they can be obtained without any additional effort.
The formulation of the solution in terms of a state-space representation in Section III-A4 with the general boundary term can be extended to different kinds of boundary conditions. In this paper, the techniques from Section III-B are used to realize the boundary conditions (6c). By the same technique also other types of boundary behavior can be realized, e.g. semi-permeable membranes or non-reflective cell walls . The incorporation of different boundary conditions can be seen as an extension of the general state-space description, so the computation of the general model (Steps 1-3 in Table I) remains the same. This allows the simulation of complex diffusive systems/channels in a semi-analytical manner by separating the general model from the desired boundary behavior.
Also, the interconnection of several semi-analytical models for larger biochemical structures is possible by the techniques shown in Section III-B . E.g. several cylindrical models can be used to model a cascade of blood vessels, or several spherical models can realize complex interacting cell structures.
V Conclusion and Future Work
In this paper, we presented a semi-analytical model for the diffusion of particles influenced by a vertical force and horizontal uniform drift in a cylindrical shape modeling e.g. the transport of information carrying magnetic nano-particles within a blood vessel under the influence of a magnetic field. The semi-analytical model is based on a modal expansion of spatial differential operators, and is useful e.g., for studying the signal strength of the channel impulse response. After the design of a general model for the diffusion process, the boundary behavior of the system is adjusted by a feedback loop. The validity of the model is verified by comparisons to particle-based simulation results. An analysis of the algorithm reveals several benefits (e.g. in runtime and flexibility) compared to purely numerical models.
Future work might investigate extensions to non-uniform laminar flow, more complex receiver models, and different geometries. Based on the time-variant source signals in , the excitation of the model by realistic time- and space-variant source signals will be considered.
-  N. Farsad, H. B. Yilmaz, A. Eckford, C.-B. Chae, and W. Guo, “A comprehensive survey of recent advancements in molecular communication,” IEEE Commun. Surv. Tutorials, vol. 18, no. 3, pp. 1887–1919, Feb. 2016.
-  R. L. Fournier, Basic Transport Phenomena in Biomedical Engineering. CRC press, 2017.
-  Q. A. Pankhurst, J. Connolly, S. Jones, and J. Dobson, “Applications of magnetic nanoparticles in biomedicine,” J. Phys. D: Appl. Phys., vol. 36, no. 13, pp. 167–181, 2003.
-  W. Wicke, A. Ahmadzadeh, V. Jamali Kooshkghazi, H. Unterweger, C. Alexiou, and R. Schober, “Molecular communication using magnetic nanoparticles,” in Proceedings of IEEE Wireless Communications and Networking Conference (WCNC), 2018, pp. 1–6.
-  H. Unterweger, J. Kirchner, W. Wicke, A. Ahmadzadeh, D. Ahmed, V. Jamali, C. Alexiou, G. Fischer, and R. Schober, “Experimental molecular communication testbed based on magnetic nanoparticles in duct flow,” presented at IEEE SPAWC 2018. [Online]. Available: https://arxiv.org/abs/1803.06990
-  L. Felicetti, M. Femminella, G. Reali, P. Gresele, M. Malvestiti, and J. N. Daigle, “Modeling CD40-based molecular communications in blood vessels,” IEEE Trans. Nanobiosci., vol. 13, no. 3, pp. 230–243, Jul. 2014.
-  Y. Chahibi, M. Pierobon, S. O. Song, and I. F. Akyildiz, “A molecular communication system model for particulate drug delivery systems,” IEEE Trans. Biomed. Eng., vol. 60, no. 12, pp. 3468–3483, Jun. 2013.
-  P. He, Y. Mao, Q. Liu, P. Liò, and K. Yang, “Channel modelling of molecular communications across blood vessels and nerves,” in Proc. IEEE ICC, 2016, pp. 1–6.
-  O. Levenspiel, Chemical Reaction Engineering, 3rd ed. John Wiley & Sons, 1999.
-  E. Furlani and K. Ng, “Analytical model of magnetic nanoparticle transport and capture in the microvasculature,” Phys. Rev. E, vol. 73, no. 6, p. 061919, Jun. 2006.
-  A. Nacev, C. Beni, O. Bruno, and B. Shapiro, “The behaviors of ferromagnetic nano-particles in and around blood vessels under applied magnetic fields,” J. Magn. Magn. Mater., vol. 323, no. 6, pp. 651–668, 2011.
-  S. Afkhami and Y. Renardy, “Ferrofluids and magnetically guided superparamagnetic particles in flows: a review of simulations and modeling,” J. Eng. Math., vol. 107, no. 1, pp. 231–251, 2017.
-  R. F. Curtain and H. Zwart, An Introduction to Infinite-Dimensional Linear Systems Theory. New York, NY, USA: Springer-Verlag New York, Inc., 1995.
-  R. V. Churchill, Operational Mathematics, 3rd ed. Boston, Massachusetts: Mc Graw Hill, 1972.
-  R. Rabenstein, M. Schäfer, and C. Strobl, “Transfer function models for distributed-parameter systems with impedance boundary conditions,” International Journal of Control, 2017.
-  G. Antonini, “Spectral models of lossy nonuniform multiconductor transmission lines,” IEEE Transactions on Electromagnetic Compatibility, vol. 54, no. 2, pp. 474–481, April 2012.
-  J. Deutscher and C. Harkort, “Parametric state feedback design of linear distributed-parameter systems,” International Journal of Control, vol. 82, no. 6, pp. 1060–1069, May 2009.
-  M. Schäfer, W. Wicke, R. Rabenstein, and R. Schober, “An nD model for a cylindrical diffusion-advection problem with an orthogonal force component,” in 23rd International Conference on Digital Signal Processing (DSP 2018), 2018.
-  C. E. Brennen, Fundamentals of Multiphase flow. Cambridge University Press, 2005.
-  J. Pérez Guerrero, L. G. C. Pimentel, T. Skaggs, and M. T. Van Genuchten, “Analytical solution of the advection–diffusion transport equation using a change-of-variable and integral transform technique,” Int. J. Heat Mass Transfer, vol. 52, no. 13-14, pp. 3297–3304, 2009.
-  M. Schäfer and R. Rabenstein, “A transfer function model for 1D diffusion processes,” in 10th International Workshop on Multidimensional (nD) Systems (nDS 2017), 2017.
-  H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids, 3rd ed. New York: Oxford University Press, 1946.
-  J. C. Willems, “The behavioral approach to open and interconnected systems,” IEEE Control Systems, vol. 27, no. 6, pp. 46–99, Dec 2007.