Transfer Function Models for Cylindrical MC Channels with Diffusion and Laminar Flow
Maximilian Schäfer, Wayan Wicke, Lukas Brand, Rudolf Rabenstein, Robert Schober
11 Transfer Function Models for Cylindrical MCChannels with Diffusion and Laminar Flow
Maximilian Schäfer, Wayan Wicke, Lukas Brand, Rudolf Rabenstein,and Robert Schober
Abstract
The analysis and design of advection-diffusion based molecular communication (MC) systems incylindrical environments is of particular interest for applications such as micro-fluidics and targeteddrug delivery in blood vessels. Therefore, the accurate modeling of the corresponding MC channelis of high importance. The propagation of particles in these systems is caused by a combination ofdiffusion and flow with a parabolic velocity profile, i.e., laminar flow. The propagation characteristicsof the particles can be categorized into three different regimes: The flow dominant regime where theinfluence of diffusion on the particle transport is negligible, the dispersive regime where diffusion hasa much stronger impact than flow, and the mixed regime where both effects are important. For thelimiting regimes, i.e., the flow dominant and dispersive regimes, there are well-known solutions andapproximations for particle transport. In contrast, there is no general analytical solution for the mixedregime, and instead, approximations, numerical techniques, and particle based simulations have beenemployed. In this paper, we develop a general model for the advection-diffusion problem in cylindricalenvironments which provides an analytical solution applicable in all regimes. The modeling procedureis based on a transfer function approach and the main focus lies on the incorporation of laminar flowinto the analytical model. The properties of the proposed model are analyzed by numerical evaluationfor different scenarios including the uniform and point release of particles. We provide a comparisonwith particle based simulations and the well-known solutions for the limiting regimes to demonstratethe validity of the proposed analytical model.
I. I
NTRODUCTION
Recently, the application of communication engineering principles to biomedical problemshas spawned the emerging interdisciplinary research field of molecular communication (MC).Comprehensive descriptions of MC can be found in [1], [2], while a tutorial review of theoreticalconcepts and modeling techniques is provided in [3]. MC is ubiquitous in natural biologicalsystems and has a high potential for bio-medical applications such as targeted drug delivery,health monitoring [2], [4], [5], and micro-fluidic channel design [6]. Besides medical applications,MC may be applied in industrial settings, e.g., for monitoring of chemical reactors and pipelines a r X i v : . [ c s . ET ] J u l [7]. The main difference between MC and classical communications is the means of transportof information from the transmitter (TX) to the receiver (RX). While classical communicationsystems rely on transport by electro-magnetic or acoustic wave propagation, mostly in free space,motivated by biological systems several different transport mechanisms have been considered forMC. These mechanisms include diffusion, gap-junction, and molecular motor based transport [3].In fluid environments, diffusion often occurs together with advection, which is prevalent, e.g., inblood vessels or pipelines. In this case, the particles are diffusing randomly and are additionallyaffected by a background flow. The flow in blood vessels and pipelines is characterized asPoiseuille flow, which exhibits a specific laminar flow profile with a radial dependence of theflow velocity [3].The accurate modeling of MC channels is crucial for the analysis of naturally occurring MCsystems and the design of artificial MC systems. As the analysis and design of advection-diffusionbased MC systems in cylindrical environments is of particular interest, e.g., for micro-fluidicapplications and targeted drug delivery systems, corresponding models have been extensivelystudied [6], [8]–[18]. Hereby, the most challenging aspect of the modeling is the correct in-corporation of the parabolic flow profile which introduces a coupling between the axial andcross-sectional particle distributions. Therefore, many existing models resort to the common plug-flow simplification, which assumes a uniform axial flow in the cylinder [19]. Based on Green’sfunctions, the authors of [8] present an MC channel model for a cylindrical environment withplug flow, a first-order degradation reaction, and partially absorbing boundaries. Advection anddiffusion of magnetic nano-particles in a cylinder is considered in [9], [10], where the particlesare also affected by an external magnetic force.There are only a few analytical models in the MC literature which consider Poiseuille flow.In [11], the impulse response of a three-dimensional (3D) advection-diffusion channel is derivedby approximating the laminar flow profile by a piecewise function for the axial distribution ofparticles. A Markovian-based channel model is presented in [12], where the cross section of thecylinder is divided into rings and the laminar flow profile is approximated by the mean of theflow velocity in each ring. A heuristic parametric model is proposed in [13] for micro-fluidicMC channels with surface-based receivers. For modeling the influence of Poiseuille flow on thepropagation of particles, it is convenient to categorize the transport process into three differentregimes, namely the flow dominant, dispersive, and mixed regimes [3, Sec. II-B], [17]. In [15],an analytical model for the flow dominant regime is presented for both uniform and point release of particles. The effect of diffusion is neglected in this regime. Dispersion is used in [16] tomodel the interplay of diffusion and laminar flow, which is also known as Taylor dispersionwhere an effective diffusion coefficient is utilized together with a plug flow approximation [17].The resulting model is applicable in the dispersive regime, where the interaction of diffusionand laminar flow yields a uniform distribution of particles in the cross section [15]. Thisapproximation for particle propagation is applied for the modeling of MC channels in, e.g., [6],[14], [15], [18]. In the mixed regime, the particles are affected by both diffusion and laminar flowand neither is negligible. Therefore, the solutions for the flow dominant and dispersive regimesare not applicable and either numerical techniques or the simplified models in [11]–[13] havebeen employed. To the best of the authors’ knowledge, a general analytical model for cylindricalMC channels with diffusion and laminar flow, which is applicable in all three regimes, has notbeen reported, yet.In this paper, we establish a general analytical model for the transport of particles by diffusionand laminar flow in cylindrical MC channels, see Fig. 1. The starting point for the modeling isthe well-known advection-diffusion equation, a partial differential equation (PDE). Subsequently,a transfer function model (TFM) is established. The TFM approach is based on the modalexpansion of a PDE into a set of eigenfunctions and eigenvalues, and provides a representationof the problem in a spatio-temporal transform domain [20], [21]. Finally, the solution of the PDEis represented as the output of a state-space description (SSD) and in terms of a concentrationGreen’s function (CGF) [22]. The TFM approach has been applied for the modeling of cylindricaland spherical MC systems [9], [23], where it has been used to realize complex boundaryconditions. However, laminar flow was not considered in [9], [23]. Therefore, in this paper,the TFM approach is extended to incorporate the influence of laminar flow. To this end, thePDE is first reduced to a simple diffusion equation which is solved in terms of an SSD of anopen loop system. Then, the influence of laminar flow is incorporated via a feedback system thatis attached to the open loop SSD to form a closed loop SSD. The design of feedback systemsis well known in control theory, see e.g., [24]. Here, this approach is adopted to incorporate theinfluence of laminar flow. The main contributions of this paper can be summarized as follows: • We derive a general analytical model for the transport of particles by diffusion and laminarflow in cylindrical MC channels, which is applicable in all three particle propagationregimes. • The proposed model can be formulated either in terms of a CGF for analytical analysis or v ( r ) e z zZ y xr R ϕ Fig. 1:
Diffusing particles (blue circles) in a cylinder of volume V subject to horizontal laminar flow v ( r ) . Thecylinder has radius R , length Z , and the boundaries are fully reflective. an SSD for efficient numerical evaluation. • Uniform and point release of the particles are considered for analysis. Comparisons withresults from particle based simulations (PBS) and known solutions for the flow dominant[15] and dispersive regimes [16] verify the validity of the proposed model.The remainder of this paper is structured as follows: Section II presents the considered advection-diffusion problem and introduces its mathematical description. Section III establishes a TFM ofthe 3D diffusion process, i.e., the open loop SSD. In Section IV, the influence of laminar flowis incorporated via a feedback system that is attached to the open loop SSD to form a closedloop SSD. The validity of the derived model is verified in Section V via numerical evaluation.Section VI further analyses the proposed model and discusses its practical implementation.Finally, Section VII concludes the paper and presents several topics for further research.II. S
YSTEM M ODEL AND M ATHEMATICAL P RELIMINARIES
A. System Model
The cylindrical volume V shown in Fig. 1 can be characterized by vector x = [ r , ϕ, z ] incylindrical coordinates and its radial boundary ∂ V r and axial boundary ∂ V z as follows V = { x = [ r , ϕ, z ] T (cid:12)(cid:12) ≤ r ≤ R , − π ≤ ϕ ≤ π, ≤ z ≤ Z } , (1) ∂ V r = { x = [ r , ϕ, z ] T (cid:12)(cid:12) r = R , − π ≤ ϕ ≤ π, ≤ z ≤ Z } , (2) ∂ V z = { x = [ r , ϕ, z ] T (cid:12)(cid:12) ≤ r ≤ R , − π ≤ ϕ ≤ π, z = , Z } . (3)The diffusion and flow of particles in V are described by an initial-boundary value problem(IBVP) in terms of the particle concentration p ( x , t ) in mol / m and the vector of particle flux i ( x , t ) in mol m − s − . The IBVP consists of a set of PDEs defined on (1), a set of boundaryconditions (BCs) defined on (2), (3), and a set of initial conditions (ICs) defined on (1) for t = . The PDE that describes the particle concentration p ( x , t ) in volume V under the influence ofdiffusion and flow is the advection-diffusion equation [25, Eq. (5.22)] ∂∂ t p ( x , t ) = D div ( grad ( p ( x , t ))) − div ( p ( x , t ) v ( x )) , (4)where ∂ / ∂ t denotes the partial derivative with respect to time and the operators div (·) and grad (·) denote the divergence and gradient operators in cylindrical coordinates, respectively. Constant D is the diffusion coefficient in m s − and v ( x ) is the velocity vector. For the case that theconsidered scenario in Fig. 1 represents a straight channel with no-slip boundary conditions, thevelocity profile is referred to as Poiseuille flow. Assuming that the channel in Fig. 1 containsa Newtonian fluid with viscosity η , the velocity vector simplifies to a radius dependent laminarflow velocity [25, Ch. 3] v ( x ) = v ( r ) e z , v ( r ) = v (cid:32) − r R (cid:33) , (5)where e z is the unit vector in z -direction and v = v eff , while v eff is the mean velocityin the channel. Decomposing the PDE in (4) into a continuity equation and a concentrationgradient equation and exploiting (5) yields a set of two PDEs describing the dynamics of particleconcentration p ( x , t ) and flux i ( x , t ) in volume (1) as follows ∂∂ t p ( x , t ) + div ( i ( x , t )) = f s ( x , t ) , x ∈ V , < t ≤ ∞ , (6) i ( x , t ) + D grad ( p ( x , t )) − v ( r ) p ( x , t ) e z = , x ∈ V , < t ≤ ∞ , (7)where the vector of fluxes i ∈ R × contains the components of the three coordinate directions i ( x , t ) = (cid:104) i r ( x , t ) i ϕ ( x , t ) i z ( x , t ) (cid:105) T (8)with (·) T denoting transposition. Function f s in (6) denotes a space and time-dependent sourcefunction that can be used to model particle injection into the channel. In addition to PDEs (6)and (7), a set of boundary and initial conditions is defined p ( x , t ) (cid:12)(cid:12) z = , Z = , x ∈ ∂ V z , < t ≤ ∞ , (9) i r ( x , t ) (cid:12)(cid:12) r = R = , x ∈ ∂ V r , < t ≤ ∞ , (10) p ( x , t ) (cid:12)(cid:12) t = = p init ( x ) , x ∈ V , t = . (11)The boundary conditions of the cylinder in z -direction (9) imply a cylinder with absorbingboundary, i.e., particles can leave the cylinder at z = , Z . We note that mostly channels ofinfinite length have been considered in the literature, e.g., [8], [9]. However, due to the applied modeling approach (see Section III) a bounded domain has to be chosen. The radial boundariesof the cylinder are fully reflective (10), and therefore the particle flux i r in radial direction iszero on ∂ V r . The initial distribution of the particles p init in V is defined by the IC (11). B. Vector Formulation
For the application of the proposed modeling approach in Section III, PDEs (6), (7) and initialconditions (11) are reformulated in a unified vector formulation as follows [22] (cid:20) ∂∂ t D − L (cid:21) y ( x , t ) = f e ( x , t ) + v flow ( x , t ) , x ∈ V , < t ≤ ∞ , (12) L = A + ∇ B , (13) y ( x , t ) (cid:12)(cid:12) t = = y init ( x ) , x ∈ V , t = , (14)where the vector of variables y ∈ R × contains the physical quantities of the PDEs (6), (7) y ( x , t ) = (cid:104) p ( x , t ) i T ( x , t ) (cid:105) T . (15)The temporal derivatives are captured by a temporal differential operator ∂∂ t D including capaci-tance matrix D ∈ R × , and the spatial differential operator L ∈ C × is composed of parametermatrix A ∈ R × and operator ∇ B containing spatial derivatives. Matrices D , A , and operator ∇ B are given by D = , A = − I , ∇ B = − D grad − div , (16)with identity matrix I ∈ R × . In accordance with the vector of variables y in (15), the sourcefunction f s in (6) is arranged into the vector f e ∈ C × in (12) and initial condition (11) intovector y i in (14) as follows f e ( x , t ) = (cid:104) f s ( x , t ) (cid:105) T , y init ( x ) = (cid:104) p init ( x ) (cid:105) T . (17)The vector valued flow term in (7) is included in vector v flow ∈ C × and is moved to the righthand side of (12) v flow ( x , t ) = (cid:104) v ( r ) p ( x , t ) (cid:105) T . (18) C. Laplace Transformation
Before the proposed TFM can be derived, the mathematical time domain description (12)is transformed into the continuous frequency domain. Application of the one-sided Laplacetransform
L{·} to (12) - (14) yields an equivalent description in the frequency domain [ s D − L ] Y ( x , s ) = F e ( x , s ) + V flow ( x , s ) + D y init ( x ) , x ∈ V , s ∈ C , (19)where variables in the continuous frequency domain are denoted by upper-case letters and dependon the complex frequency variable s , i.e., L{ y ( x , t )} = Y ( x , s ) . The temporal derivatives ∂∂ t in(12) are replaced by multiplications with complex frequency variable s ∈ C . The term D y init onthe right hand side of (19) arises from the ICs (11).III. O
PEN L OOP T RANSFER F UNCTION M ODEL
In this section, the vector V flow containing the term responsible for laminar flow is omitted,which simplifies (19) to [ s D − L ] Y ( x , s ) = F e ( x , s ) + D y init ( x ) , x ∈ V , s ∈ C . (20)Omitting V flow reduces (19) to a 3D diffusion problem. In the following, after some initialremarks regarding the modeling approach, its individual components are introduced. Then, theproposed approach is applied to (20) with BCs (9), (10) yielding a model for 3D diffusion inthe cylinder. The derived model is formulated in terms of an SSD constituting the open loopsystem that forms the basis for the incorporation of V flow in Section IV. A. Initial Remarks
The applied modeling approach is based on the modal expansion of the vector PDE (20)into an infinite set of bi-orthogonal eigenfunctions K ( x , µ ) ∈ C × and ˜ K ( x , µ ) ∈ C × wherethe functions K are the primal eigenfunctions and ˜ K are their adjoints [22]. Furthermore, eacheigenfunction K , ˜ K is associated with an eigenvalue s µ , where the infinitely many eigenvaluesdefine the discrete spectrum of the spatial differential operator L [20], [21]. Although the exactform of the eigenvalues and eigenfunctions will be derived later in Section III-D, index µ ∈ Z is already introduced here to count the eigenvalues and eigenfunctions.The infinite number of eigenvalues and eigenfunctions is necessary to ensure convergenceof the analytical solution. Nevertheless, for numerical evaluation in Section V only a finitenumber of eigenvalues can be considered. Therefore, the number of eigenvalues is truncated to µ = , . . . , Q − in the subsequent sections [26, Chap. 4.8]. B. Forward and Inverse Transformation
The solution Y of (20) in terms of an SSD is established by the application of a pair oftransformations that are introduced in the following. Their main purpose is to transform thespatial derivatives in operator L , where the transformation should have a similar effect as theLaplace transform does for the temporal derivatives. Based on well-known concepts from operatortheory and functional analysis, the proposed spatial transformations constitute a forward and aninverse Sturm-Liouville transformation (SLT) [20].
1) Forward Transformation:
Forward transformation
T {·} performs an expansion of Y ( x , s ) into a set of Q adjoint eigenfunctions ˜ K ( x , µ ) , where an individual expansion coefficient ¯ Y ( µ, s ) can be defined in terms of a scalar product or an integral in V ¯ Y ( µ, s ) = (cid:10) DY ( x , s ) , ˜ K ( x , µ ) (cid:11) = ∫ V ˜ K H ( x , µ ) DY ( x , s ) d x , (21)where (·) H denotes the conjugate-transpose. Arranging the adjoint eigenfunctions ˜ K ( x , µ ) into amatrix ˜ C ( x ) ∈ C × Q and the expansion coefficients ¯ Y ( µ, s ) into a vector ¯ Y ( s ) ∈ C Q × , ¯ Y ( s ) = (cid:2) ¯ Y ( , s ) , . . . , ¯ Y ( Q − , s ) (cid:3) T , ˜ C ( x ) = (cid:2) ˜ K ( x , ) , . . . , ˜ K ( x , Q − ) (cid:3) , (22)the forward transformation T {·} is defined in terms of a vector-valued scalar product
T { Y ( x , s )} = ¯ Y ( s ) = (cid:10) DY ( x , s ) , ˜ C ( x ) (cid:11) = (cid:10) DY ( x , s ) , ˜ K ( x , ) (cid:11) ... (cid:10) DY ( x , s ) , ˜ K ( x , Q − ) (cid:11) , (23)wherein matrix ˜ C acts as transformation kernel. In the following, variables that are transformedwith (23) are denoted by an overbar and by the variable’s dependence on µ in the scalar case.
2) Differentiation Theorem:
The most important part of the forward transformation is the def-inition of a suitable differentiation theorem enabling the replacement of the spatial derivatives bythe multiplication with a frequency domain variable. Therefore, to fit the forward transformation(23), a differentiation theorem for operator L in (20) is defined as follows [22] T { L Y ( x , s )} = (cid:10) L Y ( x , s ) , ˜ C ( x ) (cid:11) = A ¯ Y ( s ) + ¯ Φ ( s ) . (24)Diagonal matrix A = diag (cid:0) s , . . . , s Q − (cid:1) ∈ C Q × Q contains all Q eigenvalues s µ , which actas spatial frequency variables for the transformation. Vector ¯ Φ ∈ C Q × arises from the BCs –analogous to the ICs in the Laplace transform (19) – and contains the transformed boundaryvalues [22]. For the problem at hand, due to the homogeneous BCs (9), (10), ¯ Φ vanishes, i.e., ¯ Φ = . Therefore, it is omitted in the subsequent derivations.
3) Application to the PDE:
Applying forward transformation (23) to PDE (20) and exploitingdifferentiation theorem (24) yields a representation of (20) in a spatio-temporal transform domain [ s D − L ] Y ( x , s ) = F e ( x , s ) + D y init (cid:12)(cid:12) (cid:104) · , ˜ C ( x )(cid:105) s (cid:104) DY ( x , s ) , ˜ C ( x )(cid:105) − (cid:104) L Y ( x , s ) , ˜ C ( x )(cid:105) = (cid:104) F e ( x , s ) , ˜ C ( x )(cid:105) + (cid:104) D y init ( x ) , ˜ C ( x )(cid:105) s ¯ Y ( s ) − A ¯ Y ( s ) = ¯ F e ( s ) + ¯ y init . (25)Here, vectors ¯ F e ∈ C Q × and ¯ y init ∈ C Q × are the transform domain representations of the externalsources F e and the initial conditions y init ( x ) .
4) Inverse Transformation:
For forward transformation (23), an inverse transformation T − {·} is defined. The inverse transformation exploits the discrete nature of the spectrum of operator L , which allows its formulation in terms of a generalized Fourier series [20], [21] T − (cid:8) ¯ Y ( s ) (cid:9) = Y ( x , s ) = Q − (cid:213) µ = N µ ¯ Y ( µ, s ) K ( x , µ ) . (26)To fit the formulation of forward transformation (23), the sum in (26) is reformulated in terms ofa matrix-valued transformation kernel C ∈ C × Q containing primal eigenfunctions K and scalingfactor N µ T − (cid:8) ¯ Y ( s ) (cid:9) = Y ( x , s ) = C ( x ) ¯ Y ( s ) , C ( x ) = (cid:20) N K ( x , ) , . . . , N Q − K ( x , Q − ) (cid:21) . (27)Scaling factor N µ originates from the bi-orthogonality of eigenfunctions K and ˜ K and is requiredfor the formulation of an inverse transformation [26, Chap. 4.7] N µ = (cid:104) DK ( x , µ ) , ˜ K ( x , µ )(cid:105) . (28) C. Eigenfunctions and Eigenvalues
To obtain an analytical solution, eigenfunctions K , ˜ K , eigenvalues s µ , and scaling factor N µ have to be derived. These variables are derived based on the underlying physical systemexploiting specific properties of Sturm-Liouville theory. These properties are not presented indetail here but can be found in the relevant literature, e.g., [20], [21], [27].
1) Eigenfunctions:
Eigenfunctions K and ˜ K are derived by evaluation of the correspondingeigenvalue problems, which are well established for PDEs in the form of (20). The eigenvalueproblem for the primal eigenfunctions K is [22] L K ( x , µ ) = s µ DK ( x , µ ) , (29) where the eigenfunctions have to fulfill homogeneous BCs on ∂ V r , ∂ V z that are closely relatedto BCs (9), (10). The evaluation of (29) strongly depends on the exact form of L which, inthe considered case, consists of gradient and divergence operators (16). In this particular case,classical separation of variables can be applied and the resulting solution is well known in thecontext of heat transfer [28], and has been recently used in the context of MC in [8]. Similar to(29), an eigenvalue problem for the adjoint eigenfunctions ˜ K can be established, but is omittedhere for brevity [22]. The resulting eigenfunctions can be organized in vector form as follows K ( x , µ ) = J n ( k n , m r ) e j n ϕ sin ( λ ν z )− D k n , m J (cid:48) n ( k n , m r ) e j n ϕ sin ( λ ν z )− D j nr J n ( k n , m r ) e j n ϕ sin ( λ ν z )− D λ ν J n ( k n , m r ) e j n ϕ cos ( λ ν z ) , ˜ K ( x , µ ) = k n , m J (cid:48) n ( k n , m r ) e j n ϕ sin ( λ ν z ) j nr J n ( k n , m r ) e j n ϕ sin ( λ ν z ) λ ν J n ( k n , m r ) e j n ϕ cos ( λ ν z ) J n ( k n , m r ) e j n ϕ sin ( λ ν z ) , (30)where J n ( x ) denotes the n th order Bessel function of the first kind, J (cid:48) ( x ) = ∂∂ x J ( x ) is thecorresponding derivative, and j = √− . The values k n , m and λ ν with indices m ∈ Z and ν ∈ Z are related to the eigenvalues s µ and will be provided in the following.
2) Eigenvalues:
Eigenvalues s µ are derived from homogeneous BCs that have to be fulfilledby K and ˜ K . The corresponding derivation is omitted for brevity, but a detailed description canbe found in, e.g., [22]. The BC (9) on ∂ V z is a condition for the first entry of y in (15), andtherefore, the first entry K of K in (30) has to fulfill a homogeneous BC on ∂ V z yielding thecondition K ( x , µ ) (cid:12)(cid:12) z = Z ! = → sin ( λ ν Z ) ! = , (31)where relation (31) is fulfilled by wave-numbers λ ν of the form λ ν = ν π Z . The BC (10) on ∂ V r is a condition for the second entry of y in (15), and therefore, the second entry K of K in (30)has to fulfill homogeneous BCs on ∂ V r yielding the condition K ( x , µ ) (cid:12)(cid:12) r = R ! = → k n , m J (cid:48) n ( k n , m R ) ! = , (32)where k n , m is the m -th real-valued root of J (cid:48) n ( k n , m R ) . Finally, the eigenvalues s µ of the systemare defined in terms of the roots k n , m and wave-numbers λ ν s µ = s n , m ,ν = − D (cid:16) k n , m + λ ν (cid:17) . (33)This equation reveals the relevance of index µ , i.e., its purpose to count the individual eigenvalues s µ . s µ depends on the k n , m , i.e., on order n and root index m in (32), and on index ν of wave-numbers λ ν in (31). Particularly, µ represents an index tupel, i.e., ( n , m , ν ) → µ .
3) Scaling Factor:
Scaling factor N µ is defined in (28) and due to the geometrical separabilityof the cylindrical system, it is separated into three individual components N µ = (cid:104) DK ( x , µ ) , ˜ K ( x , µ )(cid:105) = N r ( µ ) N ϕ ( µ ) N z ( µ ) . (34)Evaluating the scalar product in (34) in terms of an integral, see (21), the individual componentsof N µ can be obtained as follows N r ( µ ) = ∫ R J n ( k n , m r ) r d r = R k n , m = R k n , m (cid:16) k n , m − n R (cid:17) J n ( k n , m R ) k n , m (cid:44) (35) N ϕ ( µ ) = ∫ π − π e j ( n − n ) ϕ d ϕ = π (36) N z ( µ ) = ∫ Z sin λ ν z d z = λ ν = Z λ ν (cid:44) (37) D. Open Loop Transfer Function Model
Based on the preceding sections, the open loop TFM as the solution of the 3D diffusion processin (20) can be obtained. The transform domain representation of PDE (20) in (25) serves as a state equation with state vector ¯ Y . The inverse transformation (27) acts as output equation .Together, state equation (25) and output equation (27) constitute an SSD in the s -domain s ¯ Y ( s ) = A ¯ Y ( s ) + ¯ F e ( s ) + ¯ y init , (38) Y ( x , s ) = C ( x ) ¯ Y ( s ) . (39)The vector valued solution of PDE (20), i.e., vector Y in (39), contains different physical quanti-ties. For the subsequent derivations and analysis in Sections V and VI, the particle concentrationis of interest. Therefore, output equation (39) is reduced to deliver a solution for the particleconcentration in the cylinder by restricting C in (27) to its first row, i.e., vector c T ∈ C × Q P ( x , s ) = c T ( x ) ¯ Y ( s ) , (40) c T ( x ) = (cid:20) N K ( x , ) , . . . , N Q − K ( x , Q − ) (cid:21) , (41)where K ( x , µ ) = J n ( k n , m r ) e j n ϕ sin ( λ ν z ) is the first entry of K in (30). We note that the represen-tation of the solution by (38) and (40) has to be transformed into the continuous or discrete-timedomain for analysis and numerical evaluation. IV. C
LOSED L OOP T RANSFER F UNCTION M ODEL
The previously derived open loop model (38), (40) constitutes a solution for 3D diffusion inthe cylinder. In this section, the previously excluded flow term V flow is reincorporated into themodel to obtain a solution for the considered advection-diffusion process with laminar flow. Toprovide a clear starting point, (20) is rewritten with the flow term included [ s D − L ] Y ( x , s ) = F e ( x , s ) + D y init ( x ) + V flow ( x , s ) . (42)Applying forward transformation (23) to (42) leads to a representation in the spatio-temporaltransform domain s ¯ Y ( s ) = A ¯ Y ( s ) + ¯ F e ( s ) + ¯ V flow ( s ) + ¯ y init , (43)where vector ¯ V flow ∈ C Q × denotes the transform domain representation of V flow , i.e., ¯ V flow ( s ) = (cid:104) V flow ( x , s ) , ˜ C ( x )(cid:105) = ∫ V ˜ C H ( x ) V flow ( x , s ) d x . (44)A direct evaluation of (44) is not possible in closed form as V flow itself depends on the particleconcentration P , see (18). In this paper, our objective is to establish an analytical solution forthe considered advection-diffusion problem. Therefore, we express ¯ V flow in (44) in terms of thesystem states ¯ Y and suitable feedback matrices. Finally, the derived expressions will introducea feedback system which extends the open loop state equation (38) to account for laminar flow. A. Decomposition of the Flow Profile
Due to the structure of the flow profile v ( r ) in (5), it can be decomposed into a uniform flowterm and a parabolic term. Starting with (44) and exploiting the structure of ˜ C in (22) and V flow in (18) leads to a representation with separate uniform flow and parabolic terms ¯ V flow ( s ) = ∫ V ˜ c ∗ ( x ) v ( r ) P ( x , s ) d x = v ∫ V ˜ c ∗ ( x ) P ( x , s ) d x (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) = ¯ V uni ( s ) − v R ∫ V ˜ c ∗ ( x ) P ( x , s ) r d x (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) = ¯ V par ( s ) . (45)Here, vector ˜ c T ∈ C × Q is the third row of ˜ C in (22) ˜ c T ( x ) = (cid:2) ˜ K ( x , ) , . . . , ˜ K ( x , Q − ) (cid:3) , (46)where ˜ K ( x , µ ) = λ ν J n ( k n , m r ) e j n ϕ cos ( λ ν z ) denotes the third entry of ˜ K in (30). In the following,the terms ¯ V uni and ¯ V par are considered separately. B. Uniform Flow
In the following, the uniform flow term ¯ V uni is reformulated in terms of a feedback matrixand the open loop system states ¯ Y . The starting point is the representation of ¯ V uni in terms ofan integral in (45) ¯ V uni ( s ) = v (cid:104) P ( x , s ) , ˜ c T ( x )(cid:105) = v ∫ V ˜ c ∗ ( x ) P ( x , s ) d x . (47)Now, particle concentration P is expressed in (40) in terms of the system states ¯ Y and vector c T . Inserting this representation into (47) leads to ¯ V uni ( s ) = v ∫ V ˜ c ∗ ( x ) c T ( x ) d x ¯ Y ( s ) = v K uni ¯ Y ( s ) , (48)where feedback matrix K uni ∈ C Q × Q is defined as K uni = (cid:10) c T ( x ) , ˜ c H ( x ) (cid:11) = ∫ V ˜ c ∗ ( x ) c T ( x ) d x . (49)Exploiting the similar structure of c and ˜ c , (48) can be simplified exploiting, e.g., integral andorthogonality theorems for the involved Bessel functions, and a closed-form expression for (49)can be obtained. Furthermore, we note that matrix K uni is independent of the flow velocity v ,but depends on the geometry of the cylinder. Thus, in practice, matrix K uni has to be calculatedonly once for a given cylinder geometry and can subsequently be scaled depending on the flowvelocity v . C. Parabolic Flow Profile
Analogous to the uniform flow term in (47), (48), the parabolic flow term ¯ V par is reformulated.Starting point is its representation in terms of an integral in (45) ¯ V par ( s ) = − v R (cid:104) P ( x , s ) r , ˜ c T ( x )(cid:105) = − v R ∫ V ˜ c ∗ ( x ) P ( x , s ) r d x . (50)Similar to (48), the particle concentration is expressed by (40), which is inserted into (50). Thisleads to a representation of ¯ V par in terms of matrix K par ∈ C Q × Q and system states ¯ Y , ¯ V par ( s ) = − v R ∫ V ˜ c ∗ ( x ) c T ( x ) r d x ¯ Y ( s ) = − v R K par ¯ Y ( s ) , (51)where K par = (cid:10) c T ( x ) r , ˜ c H ( x ) (cid:11) = ∫ V ˜ c ∗ ( x ) c T ( x ) r d x . (52)Similar to K uni , matrix K par can also be pre-calculated and depends only on the geometry of thecylinder. However, in contrast to K uni , matrix K par can not be obtained in closed form exceptfor Bessel functions of order n = . s − A − R K par K uni v C ( x ) ¯ F e ( s ) + ¯ y init Y ( x , s ) ¯ V flow ( s ) ¯ Y ( s ) Fig. 2:
Closed loop state space description of the advection-diffusion system in (6), (7) based on the derived TFMwith state equation (53) and output equation (39).
D. Closed Loop Transfer Function Model
Inserting ¯ V uni and ¯ V par into (45) and subsequently ¯ V flow into (43) leads to the closed loop stateequation with modified velocity dependent state matrix A c s ¯ Y ( s ) = A c ( v ) ¯ Y ( s ) + ¯ F e ( s ) + ¯ y init , A c ( v ) = A + v (cid:32) K uni − R K par (cid:33) , (53)while output equation (39) remains unchanged. Together, state equation (53) and output equation(39) constitute the closed loop TFM which is a solution to the considered advection-diffusionproblem in the s -domain. In Fig. 2, the complete SSD is illustrated. The figure clearly showsthat matrices K uni and K par constitute feedback matrices in the proposed SSD model.Structures as shown in Fig. 2 are well known in control theory for the design of parametricfeedback control systems [24]. In the considered scenario, the feedback structure serves adifferent purpose but its principle influence on the overall system behavior is similar. In particular,the open loop system is characterized by matrix A containing the eigenvalues of the 3D diffusionprocess in the absence of flow (see Section III). The flow term is incorporated via the twofeedback matrices which act on matrix A . Particularly, feedback matrices K uni and K par shiftthe eigenvalues in A which results in the new matrix A c . Therefore, the impact of flow hasbeen reduced to a modification of the eigenvalues in A resulting in the new state matrix A c which fully captures the behavior of the advection-diffusion system including laminar flow. E. Initial Conditions and External Sources
The closed loop state equation (53) of the advection-diffusion process contains – yet unspec-ified – source terms, i.e., functions ¯ F e and ¯ y init of the external sources and initial conditions,respectively. Both functions and their properties are discussed in the following.
1) Initial Conditions:
Via initial conditions y init in (14), an initial distribution of particles p init can be defined in the cylinder volume V . To obtain the transform domain representation ¯ y init in(53), the initial conditions y init in the space domain have to be transformed as in (25) ¯ y init = (cid:104) D y init ( x ) , ˜ C ( x )(cid:105) = ∫ V ˜ C H ( x ) D y init ( x ) d x . (54)The integral can be simplified exploiting the structure of D in (16) and y init in (17) as follows ¯ y init = ∫ V ˜ c ∗ ( x ) p init ( x ) d x = (cid:104) p init ( x ) , ˜ c T ( x )(cid:105) , (55) ˜ c T ( x ) = (cid:2) ˜ K ( x , ) , . . . , ˜ K ( x , Q − ) (cid:3) , (56)where ˜ c T ∈ C × Q is the fourth row of matrix ˜ C in (22) and ˜ K ( x , µ ) = J n ( k n , m r ) e j n ϕ sin ( λ ν z ) isthe fourth entry of ˜ K in (30).
2) External Sources:
Via function F e in (19), i.e., via function f s in (6), it is possible to modelthe spatial and temporal distributions of the injected particles. We assume that the correspondingfunction f s in (6) in the continuous-time domain is separable, i.e., f s ( x , t ) = f t ( t ) · f x ( x ) , (57)where f t models the temporal pulse shaping of an injection and f x models the spatial distribution.Similar to the initial conditions (54), the transform domain representation ¯ F e is obtained by thetransformation in (25) ¯ F e ( s ) = (cid:104) F e ( x , s ) , ˜ C ( x )(cid:105) = ∫ V ˜ C H ( x ) F e ( x , s ) d x , (58)which can be simplified by exploiting the structure of (17) and the separability assumed in (57) ¯ F e ( s ) = ∫ V ˜ c ∗ ( x ) F s ( x , s ) d x = F t ( s ) ∫ V ˜ c ∗ ( x ) F x ( x ) d x . (59) F. Relation to Green’s Function
As mentioned in Section III, the proposed modeling approach is based on the modal expansionof an IBVP. Modeling MC channels via modal expansions is well established for regular shapessuch as cylinders, see, e.g., [8], and often the channel is finally modeled in terms of a CGF.Although the presented approach differs from classical ones, especially due to its ability toincorporate complex flow profiles, the obtained solution can be related to a representation in terms of a CGF. To this end, first, the continuous-time equivalents of state equation (53) andoutput equation (39) are derived by application of an inverse Laplace transform s ¯ Y ( s ) = A c ( v ) ¯ Y ( s ) + ¯ y init , P ( x , s ) = c T ( x ) ¯ Y ( s ) , (60) (cid:2) L − {·} (cid:2) L − {·} ¯ y ( t ) = e A c ( v ) t ¯ y init , p ( x , t ) = c T ( x ) ¯ y ( t ) . (61)Function ¯ F e is omitted for the following considerations as the CGF is usually only derived withinitial conditions. Inserting ¯ y ( t ) and ¯ y init from (54) into output equation (61) leads to p ( x , t ) = c T ( x ) e A c ( v ) t (cid:104) p init ( x ) , ˜ c T ( x )(cid:105) . (62)Exploiting the integral formulation of (54) and rearranging (62) yields a representation of theconcentration in terms of a Green’s function, i.e., the CGF of the advection-diffusion problem p ( x , t ) = ∫ V g ( t , x | ξ ) p init ( ξ ) d ξ , g ( t , x | ξ ) = c T ( x ) e A c ( v ) t ˜ c ∗ ( ξ ) , (63)with spatial integration variable ξ . The Green’s function, g , in (63) is expressed in terms ofthe eigenfunctions c in (41) and ˜ c in (56), and modified state matrix A c , which includes theimpact of laminar flow. The influence of source function ¯ f e can be incorporated by convolutionwith (63). G. Interpretation in Terms of Transfer Functions
The proposed model is based on transfer functions. To make this fact more explicit, therepresentation in terms of an SSD is reformulated by inserting state equation (53) into outputequation (40) and solving for the concentration P ( x , s ) = c T ( x ) ¯ H ( s , D , v ) (cid:2) ¯ F e ( s ) + ¯ y init (cid:3) , ¯ H ( s , D , v ) = ( s I − A c ( v , D )) − , (64)where ¯ H denotes the transfer function. For clarity, the dependence of A c on diffusion coefficient D and flow velocity v in (53) is highlighted in the transfer function. In (64), the particleconcentration P is expressed in terms of transfer function ¯ H which is excited in the transformdomain by an input signal, i.e., by external sources ¯ F e and initial conditions ¯ y init . Hereby, transferfunction ¯ H models the influence of the cylindrical channel on the injected particles, i.e., theirpropagation based on advection and diffusion. The representation in (64) is a compact descriptionof the advection-diffusion process in the frequency domain and also allows an analysis of theprocess in terms of its spectrum. Transfer functions of the form of (64) are well known in linear operator and control theory, where they are referred to as resolvent operators that are used tostudy the spectral properties of linear operators [24], [26]. H. Discrete-time Transfer Function Model
While the previously derived representations in terms of the CGF (63) and the transfer function(64) provide compact descriptions in the continuous-time domain and the frequency domain,respectively, we also derive a representation in the discrete-time domain for numerical evaluation.To this end, an impulse-invariant transformation [29] is applied to state equation (53) and outputequation (40) to obtain a representation in the discrete-time domain ¯ y d [ k ] = A dc ( v ) ¯ y d [ k − ] + ¯ f de [ k ] + ¯ y init δ [ k ] , A dc ( v ) = e A c ( v ) T , (65) p d [ x , k ] = c T ( x ) ¯ y d [ k ] , (66)where t = kT , T is the sampling interval, and discrete-time state matrix A dc is defined in termsof a matrix exponential. Variables in the discrete-time domain are indicated by superscript (·) d and δ [ k ] denotes a delta impulse in the discrete-time domain.V. N UMERICAL E VALUATION
In this section, the proposed analytical model is numerically evaluated, i.e., its accuracyis verified by PBS, and the results are compared to existing solutions for the flow dominantand dispersive regimes [3], [15], [16]. Supplementary material including videos and figures isprovided in [30].
A. Simulation Parameters
The proposed model has been derived and is evaluated in terms of normalized physicalquantities with respect to a reference length ρ and a reference time τ , and therefore, the modelcan be applied to problems at different scales, i.e., nano, micro, or macro scale. For numericalevaluation, the parameters in Table I have been used, which may model, e.g., micro-fluidic ducts[15], [25], but can also be scaled to model small capillaries [17]. In the following, all parameters,except the diffusion coefficient D , are kept constant. For all numerical evaluations, the discrete-time SSD (65), (66) with sampling interval T = · − s was employed, and the number ofeigenvalues Q was chosen as Q = ( N + ) · M · L , q = [ N , M , L ] , (67) TABLE I: Physical parameters for numerical evaluation
Parameter Value Normalized valueRadius R µ m 1 Length Z Flow velocity v µ m s − TX/RX distance d µ m 1 Diffusion coefficient D . · − m s − . . . · − m s − . · − . . . Reference length ρ R Reference time τ · s where N denotes the maximum orders of Bessel functions J n used in (30), (32), i.e., n = − N , . . . , N , M is the number of roots k n , m in (32) for each order n , i.e., m = , . . . M − , and L denotes the number of wave-numbers λ ν in (31), i.e., ν = , . . . , L − . The selection of the valuesof N , M , and L in (68) directly affects the accuracy of the proposed model and is discussed indetail in Section VI-B. B. Initial Conditions
For the analysis and numerical evaluation of the proposed model, we consider two differentinitial distributions of the particles in the cylinder, i.e., a uniform distribution and a pointdistribution. For the initial distributions, the following raised cosine function is defined f i ( χ, χ , χ e ) = (cid:16) + cos ( πχ ( χ − χ e )) (cid:17) χ e − χ ≤ χ ≤ χ e + χ , (68)with a spatial width χ and center position χ e .
1) Uniform Distribution:
In the considered uniform distribution, the particles are uniformlydistributed in the r - ϕ -plane of the cylinder. In z -direction, the initial distribution of particles iscentered at z = z e and spread over z as defined by (68). The 3D uniform distribution is definedby specifying initial conditions p init in (11) as follows p init ( x ) = p uniform ( x ) (cid:66) f i ( z , z , z e ) , (69)with normalized width z = . and normalized center position z e = . The considered uniformdistribution is shown in the plots on the left hand side of Fig. 3. We note that instead of the raised cosine function (68), any other smooth function can be used to model the initial distributionof particles.
2) Point Distribution:
Furthermore, a point distribution centered at x TX = [ r e , ϕ e , z e ] is con-sidered. The particles are distributed as defined in (68) in all spatial directions. The 3D pointdistribution is defined by specifying initial conditions p init in (11) as follows p init ( x ) = p point ( x ) (cid:66) f i ( r , r , r e ) f i ( ϕ, ϕ , ϕ e ) f i ( z , z , z e ) , (70)with normalized widths r = . , ϕ = π , z = . . The center positions in ϕ - and z -direction are ϕ e = π and z e = , respectively, while the radial center position is varied as r e = . , . , . .The considered point distribution is shown in the plots on the left hand side of Fig. 5.Mostly, MC channel models are derived and analyzed by assuming a point release of particlesthat is defined in terms of δ -impulses, see e.g. [8, Eq. (11)], [9, Eq. (5d)] for cylindricalenvironments. However, such impulsive releases are unrealistic, as in practical systems particlescannot be released from an infinitesimal point, but the assumption simplifies the derivation ofthe channel model and the channel impulse response. To go one step towards more realisticchannel models, we consider a point release of particles over a non-zero volume (70). Anotherbenefit of the adopted release profile is the suppression of Gibbs phenomenon which otherwisemay occur for all modeling techniques based on modal expansions or CGFs [29]. C. Validation Parameters
For validation, we employ PBS of the considered advection-diffusion process. For PBS, theconcentration is estimated by counting the number of observed particles in a cuboid volume V cube = . × . × . centered at receiver position x RX . For the uniform release scenario inSection V-D, we released N TX = · particles and their positions are updated in discrete timesteps ∆ t = − s . The PBS results are averaged over realizations of the process. For thepoint release in Section V-E, N TX = · particles are released and their positions are updatedwith ∆ t = · − s and the results are averaged over · realizations.The proposed model is evaluated at point x RX , which is the center of the cuboid used for PBS.Using the uniform concentration assumption in the cuboid [31], output equation (66) becomes p dcube [ k ] = V cube · c T ( x RX ) ¯ y d [ k ] . (71)Furthermore, the proposed model is classified with respect to existing analytical models. There-fore, the dispersion factor α , is introduced [3, Eq. (20)] α = D d v eff R = D d v R (72) TABLE II:
Considered diffusion coefficients and resulting dispersion factors D in m s − . · − . · − . · − . · − . · − . · − · − α · − · − . . . to distinguish between different regimes. As summarized in [3, Sec. D-2] the transport of particlesin the considered scenarios can be categorized into three regimes, i.e., the flow dominant regime( α (cid:28) ), the dispersive regime ( α (cid:29) ), and the mixed regime. Furthermore, there are well-known solutions and approximations for the flow dominant regime, see [15, Eq. (16)], and thedispersive regime, see [3], [15, Eq. (11)], [16].In the following, dispersion factor α is used to classify the considered scenarios into differentregimes where the different regimes are realized by changing the diffusion coefficient D , seeTable II. We note that other parameters may also be varied to evaluate the model in differentregimes, see (72). Varying the diffusion coefficient D allows, e.g., to analyze the behavior ofparticles of different sizes in a given channel. D. Uniform Release
In this section, the proposed model is evaluated for a uniform release of particles modeledby initial condition (69). For numerical evaluation of the proposed model, a total number of Q = eigenvalues is used, i.e., q = [ , , ] . Here, N = is due to the initial conditionin (69), i.e., as the initial distribution is radially symmetrical, only Bessel functions of order n = contribute to the solution. The receiver is placed at x RX = [ , π / , ] , where the concentration iscalculated based on (71).In Fig. 3, the concentration of the particles after a uniform release is presented for the flowdominant ( α = · − , top row), mixed ( α = . , center row), and dispersive ( α = , bottomrow) regime at times t = , , and of the process. The figure illustrates the differencesbetween the defined regimes. In the flow dominant regime (see top row of Fig. 3), the influenceof flow is dominant and diffusion has little impact. The characteristic parabolic profile of thelaminar flow, see (5), becomes obvious over time, with maximum velocity v ( ) = v in the centerof the cylinder and zero velocity v ( R ) = at the boundaries. The spatial spreading of the initialdistribution is preserved at the considered RX position (see also Fig. 4a). In the mixed regime(see center row of Fig. 3), the flow profile is blurred by diffusion. In fact, both flow and diffusioninfluence the propagating particles. With increasing distance from the TX position, the initial − R − R R R × RX t = 0 s α = 10 − y ( n o r m a li ze d ) × RX t = 2 s α = 10 − × RX t = 4 s α = 10 − − R − R R R × RX t = 0 s α = 0 . y ( n o r m a li ze d ) × RX t = 2 s α = 0 . × RX t = 4 s α = 0 . − R − R R R × RX t = 0 s α = 2 y ( n o r m a li ze d ) × RX t = 2 s α = 2 × RX t = 4 s α = 2 Fig. 3:
2D concentration p ( x , t ) in the y - z -plane ( y = r sin ( π ) ) of the cylinder at times t = , , for a uniformrelease. Different values α = · − , . , are employed (top to bottom). uniform distribution is spread over space. This effect becomes even stronger in the dispersiveregime (see bottom row of Fig. 3), where the impact of diffusion dominates the impact of flowon the propagating particles. The initial distribution of the particles is noticeably spread overspace already after t = (see second plot in bottom row of Fig. 3).Fig. 4 shows the particle concentration at RX position x RX for different values of α . Thefigure shows results for the numerical evaluation of the proposed model (red color) and PBSas a ground truth (gray color). Furthermore, the existing solutions for the flow dominant (bluecolor) and dispersive (green color) regimes are shown. The most important observation fromFig. 4 is that the proposed model perfectly matches the PBS results for all considered regimes,which underlines the ability of the model to provide a solution valid for all regimes. In bothlimiting cases, the existing solutions for the flow dominant regime (Fig. 4a) and the dispersiveregime (Fig. 4f) also provide a good estimate for the received concentration. The plots in Fig. 4highlight the influence of the different regimes on the propagation of the particles, and reinforcethe observations made in Fig. 3. The peakiness of the uniform particle release profile is evidentin the flow dominant regime (Fig. 4a), and the temporal width t peak of the released concentrationcan be related to the spatial width of the uniform release, i.e., t peak ≈ z v = . . For α = · − and . (Figs 4b, c) the peak is spread by diffusion, but still recognizable. In both figures, themismatch between the known solutions for the limiting regimes and the results from PBS areobvious. The effect of diffusion starts to become dominant for α = . and . in Figs. 4d, ! . . . . . · ! C o n ce n tr a t i o n PBSFlow regimeDispersive regimeTFM model (a) α = · − ! . . . . . · ! PBSFlow regimeDispersive regimeTFM model (b) α = · − · ! PBSFlow regimeDispersive regimeTFM model (c) α = . · − Time in s C o n ce n tr a t i o n PBSFlow regimeDispersion regimeTFM model (d) α = . · ! Time in s PBSFlow regimeDispersive regimeTFM model (e) α = . · ! Time in s PBSFlow regimeDispersive regimeTFM model (f) α = Fig. 4:
Concentration p dcube at x RX over time for a uniform release of particles (69). Different values of α areconsidered. e. In this case, the tail of the received concentration increases, which is directly related to thespatial spreading of the uniform release, see center row of Fig. 3. Furthermore, both figures showthat for larger values of α the known solution for the dispersive regime (green color) starts toprovide a better estimate for the received concentration. For α = (dispersive regime) in Fig. 4f,diffusion clearly dominates. The peak is completely spread and, compared to the other scenariosin Fig. 4, even the temporal location of the maximum received concentration occurs earlier dueto the high diffusion. In this scenario, the known solution for the dispersive regime provides agood estimate for the concentration. E. Point Release
In this section, a point release of the form in (70) is considered. For numerical evaluation Q = . · is used, i.e., q = [ , , ] . Because of the point release, a large value of N is necessary to correctly represent the propagation of the particles in the 3D volume, seeSection VI-B. Due to the symmetry of Bessel functions, i.e., J − n ( x ) = (− ) n J n ( x ) , the negativevalues of n = − N , . . . , − are not evaluated separately, which reduces (68) to Q = ( N + ) · M · L . − R − R R R × RX t = 0 s α = 10 − y ( n o r m a li ze d ) × RX t = 2 s α = 10 − × RX t = 4 s α = 10 − − R − R R R × RX t = 0 s α = 0 . y ( n o r m a li ze d ) × RX t = 2 s α = 0 . × RX t = 4 s α = 0 . − R − R R R × RX t = 0 s α = 1 y ( n o r m a li ze d ) × RX t = 2 s α = 1 × RX t = 4 s α = 1 Fig. 5:
2D concentration p ( x , t ) in the y - z -plane ( y = r sin ( π ) ) of the cylinder at times t = , , for a pointrelease. Different values α = · − , . , are employed (top to bottom). Point releases at three different release positions with r e = . , . , and . are considered,while the receiver position is fixed at x RX = [ . , π / , ] .Fig. 5 shows the concentration for a point release at r e = . in the cylinder for differentvalues of α . Although the general propagation behavior for a point release is similar to that for auniform release (see Fig. 3), Fig. 5 highlights some differences and reveals significant effects ofpractical relevance. The top row of Fig. 5 shows the propagation for α = · − . As previouslydiscussed, flow dominates diffusion in this case and the initial point is not spread spatially, butits initial shape is still distorted over time due to the parabolic flow profile, see top row of Fig. 3.For α = . the point starts to spread by diffusion, see middle row of Fig. 5. Due to the radialrelease position at r e = . , and the the zero flow v ( R ) = at the boundary, a certain percentageof particles accumulate at the cylinder wall for t = where they are only affected by diffusion.Furthermore, particles start to propagate into the lower part of the cylinder. This effect is evenmore pronounced for α = , see bottom row of Fig. 5. In this case, after t = , the particlesare distributed over the entire radial plane of the cylinder. The effects that arise for increasing α and t raise the question for which α the initial TX position is completely forgotten for a givenRX position.This question is further investigated in Fig. 6, which shows the concentration (71) at position x RX for different values of α after a point release on the radial axis for r e = . , . and . .The figure shows the numerical evaluation of the proposed model (red, blue, green colors) and for (a) α = · − (b) α = . (c) α = . Fig. 6:
Concentration p cube at x RX over time for a point release of particles at radial release positions r e = . , . , . . Different values of α are considered. comparison results from PBS (gray color). For α = · − (flow dominant regime, Fig. 6a), thedifferences in the received concentration caused by different release positions are clearly visible.For α = . (mixed regime, Fig. 6b) the received concentration starts to increase simultaneouslyfor r e = . and r e = . , while the received concentration starts to increase later for r e = . ,which is due to the zero flow at the boundary. For α = . (mixed regime, Fig. 6c) the receivedconcentration increases simultaneously for all considered release positions r e . In this scenario,the previously mentioned effect becomes evident, i.e., it is not possible to determine the releaseposition r e based on the received concentrations. This effect becomes even more pronounced forlarger α , see supplementary material in [30, Sec. 3]. We note that the numerical results obtainedwith the proposed TFM perfectly match the PBS results for all scenarios considered in Fig. 6.VI. I MPLEMENTATION AND A NALYSIS
The previous section has shown that the evaluation of the proposed model perfectly matchesthe results obtained with PBS for all considered scenarios. The proposed TFM is an analyticalsolution for the advection-diffusion problem and provides a compact description valid for allregimes. Therefore, the model closes the gap between existing solutions for the flow dominantand dispersive regimes.In this section, we provide a short overview of the implementation of the proposed discrete-time SSD (65), (66). Furthermore, the accuracy of the model is analyzed and its limitations arediscussed. Finally, the benefits of the proposed model are shortly summarized. Cylinder: R , Z Diffusion: D Advection: v k n , m , λ ν s µ K uni , K lam N µ c T ( x ) A A c ( v ) (32), (31)(34) (33)(40) (24) (53)(49), (52) (53) Fig. 7:
Calculation of vectors and matrices in the SSD (65), (66) and their dependencies on physical parameters D , v , and geometrical parameters R , Z . A. Remarks on the Implementation
To analyze the dynamics of the particles in the cylinder, the derived discrete-time model in(65), (66) has to be numerically evaluated. The dependence of variables c and A c in (65), (66)on the physical and geometrical parameters is illustrated in Fig. 7. In the figure, green boxesindicate a dependence on geometrical parameters R , Z , red boxes a dependence on the diffusioncoefficient D , and blue boxes a dependence on the flow velocity v . As can be observed, thevalues k n , m and λ ν only depend on the geometry of the cylinder and can be computed independentfrom diffusion coefficient D and flow velocity v . The same applies for both feedback matrices K . By exploiting these limited dependencies, many calculations needed for the evaluation ofthe SSD (65), (66) can be performed once in advance and do not have to be repeated if theparameters change.Depending on the number of eigenvalues Q , see (68), a straightforward implementation of theSSD (65), (66) may lead to high computational costs. Particularly, the calculation and subsequentmultiplication of the potentially fully occupied matrix A dc and system states ¯ y d in (65) istime consuming. Therefore, state equation (65) should be modified to speed up the requiredmultiplications. To this end, the block diagonal structure of matrix A dc can be exploited, i.e.,the matrix consists of ( N + ) blocks of size ( M · L ) × ( M · L ) . This block structure allowsthe block-wise calculation of (65). Each of the resulting blocks can be further simplified by aneigendecomposition which allows for fast evaluation by a parallel structure of filters. Furthermore,due to the symmetry of the Bessel functions, the number of blocks can be reduced to ( N + ) ,see Section V-E. These modifications can be applied to enable fast evaluation of the proposedmodel in, e.g., MATLAB. They are further described in the supplementary material providedtogether with the MATLAB code in [30]. · − Time in s C o n ce n tr a t i o n PBS Q = 6000 : q = [0 , , Q = 3000 : q = [0 , , Q = 500 : q = [0 , , Q = 100 : q = [0 , , Q = 20 : q = [0 , , (a) α = . , uniform release . . . . . . . · − Time in s C o n ce n tr a t i o n PBS Q = 12 . · , q = [20 , , Q = 6 . · , q = [20 , , Q = 2 . · , q = [20 , , Q = 5 · , q = [5 , , Q = 2 · , q = [2 , , (b) α = · − , point release · − Time in s C o n ce n tr a t i o n PBSFlow Regime Q = 6000 : q = [0 , , Q = 15 · : q = [0 , , · (c) α = · − , uniform release Fig. 8:
Accuracy analysis for (a) uniform release with α = . , (b) point release width α = · − , and (c)uniform release with α = · − . B. Analysis of Accuracy
Although the proposed model involves an infinite sum, see (26), it represents an analyticalsolution for the considered advection-diffusion process. Mathematically, the solution only con-verges if the number of eigenvalues Q → ∞ . For numerical evaluation and analysis, this numberhas to be restricted to a finite value. This implies a trade-off between the complexity and accuracyof the model, see also [9, Section IV-C]. In practice, the number of eigenvalues Q has to bechosen such that the accuracy requirements of the desired scenario are met. For example, ifonly a rough impression of the system behavior is desired, a low number Q would be sufficient.For an accurate evaluation of the concentration in the complete volume, a higher value of Q isneeded.In Figs. 8a and b, the proposed model is evaluated for different values of Q and PBS resultsare provided as ground truth. The scenarios considered in Figs. 8a and b are identical to thosein Figs. 4c and 6a ( r e = . ). The accuracy for uniform release and point release are analyzedseparately because for the uniform release only Bessel functions of order n = contribute to thesolution. For uniform release, the proposed model converges to the PBS result for Q = .Reducing the number of eigenvalues to Q = or even Q = , the proposed model is stillin good agreement with the PBS results. Only for Q = and Q = a significant differencecan be observed. In the point release case, a higher number of eigenvalues is necessary forconvergence. Fig. 8b shows that the proposed model converges to the PBS results for Q = . · , . . . , . · eigenvalues. For smaller values Q = · and Q = · , the proposedmodel is not converging to the PBS results. Compared to a uniform release, larger values of Q are necessary for convergence for a point release because also Bessel functions of order n (cid:44) contribute to the solution. Particularly, large values of N are necessary to correctly represent theparticle propagation in the volume (see videos for small values of N in [30, Sec. 4.3]).Fig. 8c illustrates a limitation of the proposed model, i.e., of its numerical evaluation, arisingfor very low values of α , e.g., α = · − . The figure shows the numerical evaluation of theproposed model for different values of Q while PBS results and the known solution for the flowdominant regime are provided as ground truth. From t = up to t ≈ . , it can be observedthat the evaluation of the proposed model still perfectly matches the PBS results and the knownsolution for the flow dominant regime. Particularly, the amplitude, the duration, and the shapeof the received concentration from a uniform release are captured by the proposed model, seethe zoomed excerpts in Fig. 8c, while undesired oscillations occur after the peak. As explainedin Section IV, in the proposed model, the influence of flow is incorporated by a shift of theeigenvalues, the accuracy of which depends on the number of eigenvalues Q . However, for verysmall values of α , the influence of flow dominates and therefore, the shift of the eigenvaluesis very large and the proposed model is not converging for the considered values of Q . Byincreasing the number of eigenvalues significantly to Q = · (green curve in Fig. 8c), theamplitude of the oscillations can be reduced but the effect cannot be suppressed completely.The analytical form of the proposed model, e.g., in terms of a CGF in (63), is not necessarilyrestricted to a finite number of eigenvalues Q . Therefore, the limitation for very small values of α does not affect the validity of the proposed analytical model, but only its numerical evaluationwhere Q has to be finite. Hence, for very small values of α , using the known solution for the flowdominant regime may be preferable, see, e.g., [15, Eq. (16)]. As a rule of thumb, the number ofeigenvalues Q required for an accurate representation decreases for increasing values of α andfor increasing spatial spreading of the initial particle distribution. This is because fewer spatialeigenfunctions are needed to approximate smooth functions compared to peaky ones.In Section II, the cylinder is bounded in z -direction by BCs (3) and (9). Particularly, BC (9)corresponds to an absorbing boundary, and therefore all particles leave the cylinder for t → ∞ .This is true for the analytical formulation of the proposed model, but for its numerical evaluationsome undesired effects occur due to the numerical restriction of the z -direction to Z , see videoin [30, Sec. 4.1]. Instead of leaving the cylinder at z = Z , particles are reflected and re-enterthe cylinder at z = where the re-entering is accompanied by undesired reflections in thecylinder. In future work, we plan to overcome this effect by adopting the techniques proposed in[32]. For numerical evaluation of the model in the proposed form, these effects can be avoided, TABLE III:
Different formulations of the proposed model.
PDE ∂∂ t p ( x , t ) = D div ( grad ( p ( x , t ))) − v ( r ) div ( p ( x , t ) e z ( x )) Eq. (4),(5) continuous-time domainCGF p ( x , t ) = ∫ V g ( t , x | ξ ) p i ( ξ ) d ξ Eq. (63) continuous-time domainTFM P ( x , s ) = c T ( x ) ¯ H ( s , D , v ) (cid:2) ¯ F e ( s ) + ¯ y i (cid:3) Eq. (64) frequency domainSSD ¯ y d [ k ] = A dc ( v ) ¯ y d [ k − ] + ¯ f de [ k ] + ¯ y i δ [ k ] p d [ x , k ] = c T ( x ) ¯ y d [ k ] Eq. (65)Eq. (66) discrete-time domain e.g., by ensuring that no particles leave the cylinder during the observation time of the system.Therefore, we restricted the observation time to t obs ≤ Z − z e v =
18 s for the numerical evaluationin the considered scenarios.
C. Benefits of the Proposed Model
In Section V, the proposed model matches the PBS results for all considered scenarios.Despite the previously mentioned limitations for the numerical evaluation of the proposed model,it provides a general analytical description of the advection-diffusion process with laminarflow. Section IV introduced different equivalent formulations of the proposed model which aresummarized in Table III. The CGF in (63) and the representation in terms of transfer functionsin (64) provide an analytical description of the MC channel and allow a representation of thechannel response in analytical form for given TX models. The discrete-time SSD in (65), (66)provides a suitable model for numerical evaluation. The dependence of the convergence on thenumber of eigenvalues Q (see Section VI-B) also implies flexibility. Particularly, the value of Q can be adjusted depending on the objective of the investigation. To get a rough impression of thechannel behavior for a given input signal, a small Q is sufficient to perform many simulations in ashort time, while a large value of Q can be chosen for accurate simulation of particle propagation.The initial formulation of the model in terms of an SSD (see (40), (53)) can also be extendedas has been discussed for similar models in [9, Sec. IV-D]. For example, one possibility isthe incorporation of more complex or time-varying boundary conditions, which is described indetail in [22], [26], and has been applied to cylindrical and spherical MC systems in [9], [23].Furthermore, the SSD allows the interconnection of multiple systems (see [23, Sec. V]), which isbeneficial for the modeling of interconnected tube systems or cascades of blood vessels. Finally,SSD models can be exploited for derivation of parameter estimation algorithms based on Kalmanfilters [33] to determine relevant system parameters from measurements. VII. C
ONCLUSION
In this paper, we have proposed an analytical model for advection-diffusion processes incylindrical environments affected by laminar flow. The proposed model has been derived based ona transfer function approach, which provides a general, flexible, and extendable description of theMC channel, and can be formulated in terms of a CGF or an SSD depending on the application.The validity of the proposed model has been verified by numerical evaluation. Particularly, ithas been shown that the proposed solution matches the PBS in all considered regimes andcorresponds to the known solutions for the flow dominant and dispersive regimes. In contrastto all known models, the proposed solution is also applicable in the mixed regime where bothdiffusion and flow have a similar impact on particle propagation.The discussion of the benefits and limitations of the proposed model in Section VI suggeststhe following topics for future work: As discussed in Section VI, a large number of eigenvalues Q may be necessary to fully capture all effects of the propagation of the particles in the cylinder.The required value of Q can be reduced by a model reduction to the most dominant eigenvalues.Furthermore, to make the model even more comprehensive, the extension of the SSD to morecomplex boundary conditions, e.g., semi-permeable walls, and the inclusion of particle reactionsis an interesting direction for future work. Also, the incorporation of time-variant flows to model,e.g., the pumping of blood, and the analysis of more complex RX and TX models are of interest.R EFERENCES [1] T. Nakano, A. W. Eckford, and T. Haraguchi,
Molecular Communication . Cambridge: Cambridge University Press,2013. [Online]. Available: https://doi.org/10.1017/CBO9781139149693[2] N. Farsad, H. B. Yilmaz, A. Eckford, C. B. Chae, and W. Guo, “A Comprehensive Survey of Recent Advancements inMolecular Communication,”
IEEE Communications Surveys Tutorials , vol. 18, no. 3, pp. 1887–1919, thirdquarter 2016.[3] V. Jamali, A. Ahmadzadeh, W. Wicke, A. Noel, and R. Schober, “Channel Modeling for Diffusive Molecular Communi-cation - A Tutorial Review,”
Proceedings of the IEEE , vol. 107, no. 7, pp. 1256–1301, Jul. 2019.[4] L. Felicetti, M. Femminella, G. Reali, P. Gresele, M. Malvestiti, and J. N. Daigle, “Modeling CD40-based MolecularCommunications in Blood Vessels,”
IEEE Trans. Nanobiosci. , vol. 13, no. 3, pp. 230–243, Jul. 2014.[5] O. C. Farokhzad and R. Langer, “Impact of Nanotechnology on Drug Delivery,”
ACS Nano , vol. 3, no. 1, pp. 16–20,2009, pMID: 19206243. [Online]. Available: https://doi.org/10.1021/nn900002m[6] A. O. Bicen and I. F. Akyildiz, “System-theoretic Analysis and Least-Squares Design of Microfluidic Channels for Flow-induced Molecular Communication,”
IEEE Trans. Signal Process. , vol. 61, no. 20, pp. 5000–5013, Jul. 2013.[7] L. P. Gine and I. F. Akyildiz, “Molecular Communication Options for Long Range Nanonetworks,”
Computer Networks ,vol. 53, no. 16, pp. 2753 – 2766, 2009.[8] M. Zoofaghari and H. Arjmandi, “Diffusive Molecular Communication in Biological Cylindrical Environment,”
IEEETrans. NanoBioscience , vol. 18, no. 1, pp. 74–83, Jan 2019. [9] M. Schäfer, W. Wicke, R. Rabenstein, and R. Schober, “Analytical Models for Particle Diffusion and Flow in a HorizontalCylinder with a Vertical Force,” in Proc. IEEE Int. Conf. Commun. (ICC 2019) , Shanghai, China, May 2019, pp. 1–7.[10] 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,” in
Proc. IEEESPAWC , 2018. [Online]. Available: https://arxiv.org/abs/1803.06990[11] F. Dinc, B. C. Akdeniz, A. E. Pusane, and T. Tugcu, “A General Analytical Approximation to Impulse Response of 3-DMicrofluidic Channels in Molecular Communication,”
IEEE Trans. NanoBiosci. , vol. 18, no. 3, pp. 396–403, Mar. 2019.[12] Y. Lo, C. Lee, P. Chou, and P. Yeh, “Modeling Molecular Communications in Tubes with Poiseuille Flow and RobinBoundary Condition,”
IEEE Commun. Letters , vol. 23, no. 8, pp. 1314 – 1318, Jun. 2019.[13] M. Kuscu and O. B. Akan, “Modeling Convection-Diffusion-Reaction Systems for Microfluidic Molecular Communicationswith Surface-Based Receivers in Internet of Bio-Nano Things,”
PLOS ONE , vol. 13, no. 2, p. e0192202, Jul. 2018.[14] P. He, Y. Mao, Q. Liu, P. Liò, and K. Yang, “Channel Modelling of Molecular Communications across Blood Vessels andNerves,” in
Proc. IEEE Int. Conf. Commun. (ICC) , 2016, pp. 1–6.[15] W. Wicke, T. Schwering, A. Ahmadzadeh, V. Jamali, A. Noel, and R. Schober, “Modeling Duct Flow for MolecularCommunication,” in
Proc. IEEE Global Commun. Conf. (GLOBECOM) , 2018, pp. 206–212.[16] R. Aris and G. I. Taylor, “On the Dispersion of a Solute in a Fluid Flowing through a Tube,”
Proceedings of the RoyalSociety of London. Series A. Mathematical and Physical Sciences , vol. 235, no. 1200, pp. 67–77, Apr. 1956.[17] R. F. Probstein,
Physicochemical Hydrodynamics: An Introduction . John Wiley & Sons, 2005.[18] Y. Chahibi, M. Pierobon, and I. F. Akyildiz, “Pharmacokinetic Modeling and Biodistribution Estimation Through theMolecular Communication Paradigm,”
IEEE Trans. Biomed. Eng. , vol. 62, no. 10, pp. 2410–2420, May 2015.[19] O. Levenspiel,
Chemical Reaction Engineering , 3rd ed. John Wiley & Sons, 1999.[20] R. V. Churchill,
Operational Mathematics . Boston, Massachusetts: Mc Graw Hill, 1972.[21] R. Curtain and H. Zwart,
An Introduction to Infinite-Dimensional Systems Theory . New York: Springer-Verlag, 1995.[22] R. Rabenstein, M. Schäfer, and C. Strobl, “Transfer Function Models for Distributed-Parameter Systems with ImpedanceBoundary Conditions,”
Int. J. of Control , vol. 91, no. 12, pp. 2726–2742, Nov. 2017.[23] M. Schäfer, W. Wicke, W. Haselmayr, R. Rabenstein, and R. Schober, “Spherical Diffusion Model with Semi-PermeableBoundary: A Transfer Function Approach,” in
Proc. IEEE Int. Conf. Commun. (ICC 2020) , Dublin, Ireland, June 2020.[24] J. Deutscher,
Zustandsregelung verteilt-parametrischer Systeme . Heidelberg, Germany: Springer, 2012.[25] H. Bruus,
Theoretical Microfluidics , 1st ed. Oxford, UK: Oxford University Press, 2007.[26] M. Schäfer, “Simulation of Distributed Parameter Systems by Transfer Function Models,” doctoral thesis, Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU), 2020.[27] C. Eringen, “The Finite Sturm-Liouville-Transform,”
J. Mathematics, Oxford 2nd Series , vol. 5, pp. 120–129, 1954.[28] H. S. Carslaw and J. C. Jaeger,
Conduction of heat in solids , 3rd ed. New York: Oxford University Press, 1946.[29] B. Girod, R. Rabenstein, and A. Stenger,
Signals and Systems . West Sussex, UK: John Wiley & Sons Ltd, 1997.[30] M. Schäfer, W. Wicke, L. Brand, and R. Schober. (2020) Transfer Function Models for Cylindrical MC Channels withDiffusion and Laminar Flow. [Online]. Available: https://maximilianschaefer.org/publication/tfm-laminar[31] A. Noel, K. C. Cheung, and R. Schober, “Using Dimensional Analysis to Assess Scalability and Accuracy in MolecularCommunication,” in
IEEE Int. Conf. Commun. Workshops (ICC) , 2013, pp. 818–823.[32] J. Grant and M. Wilkinson, “Advection–Diffusion Equation with Absorbing Boundary,”
Journal of Statistical Physics ,vol. 160, no. 3, pp. 622–635, 2015. [Online]. Available: https://doi.org/10.1007/s10955-015-1257-2[33] M. Schäfer, A. Ruderer, and R. Rabenstein, “An Eigenfunction Approach to Parameter Estimation for 1D DiffusionProblems,” in