Transfer Function Models for Cylindrical MC Channels with Diffusion and Laminar Flow

Contents:

1. Preliminaries
2. Videos of the propagation of the particles in the cylinder for a uniform release supplementary to Fig. 4 (see [1, Sec. V-D])
3. Videos of the propagation of the particles in the cylinder for a point release supplementary to Fig. 5 (see [1, Sec. V-E])
4. Videos of the propagation of the particles in the cylinder for analysis of accuracy (see [1, Sec. VI-B])
5. Remarks regarding the implementation (see [1, Sec. VI-A])
6. Performance analysis
7. References

1. Preliminaries

We note that the variables in the presented formulas are defined in [1].

1.1 Summary of the proposed model

The considered advection-diffusion problem with laminar flow is mathematically described by the PDE [1, Eq. (4)] $$\begin{array}{r}\frac{\partial }{\partial t}p(\mathbf{x},t)=D\text{div}\left(\text{grad}(p(\mathbf{x},t))\right)-\text{div}(p(\mathbf{x},t)v(r){\mathbf{e}}_z),\end{array}$$ with the radius dependent flow velocity [1, Eq. (5)] $$\begin{array}{r}v(r)=v_0(1-\frac{r^2}{R_0^2}).\end{array}$$

The solution of the PDE is derived by a transfer function approach as described in Section IV in [1]. Finally, the concentration $p$ in the volume can be represented in terms of a concentration Green’s function in the time domain [1, Eq. (63)] $$\begin{array}{rl}& p(\mathbf{x},t)={\int }_{V}g(t,\mathbf{x}|\xi )p_{\mathrm{init}}(\xi )\mathrm{d}\xi ,\end{array}$$ with Green’s function \begin{align*} g(t, \mathbf{x}\vert \mathbf{\xi}) = \mathbf{c}1^{\mathrm{\scriptsize T}}(\mathbf{x}) \mathrm{e}^{\mathbf{\mathcal{A}}\mathrm{c}(v_0)t}\tilde{\mathbf{c}}^_4(\mathbf{\xi}). \end{align}

For numerical evaluation and analysis of the proposed model a formulation in the discrete-time domain, i.e., in terms of a state-space description is preferred [1, Eqs. (65), (66)] \begin{align*} \bar{\mathbf{y}}^\mathrm{d}[k] &= \mathbf{\mathcal{A}}{\mathrm{c}}^\mathrm{d}(v_0) \bar{\mathbf{y}}^\mathrm{d}[k - 1] + \bar{\mathbf{f}}^\mathrm{d}\mathrm{e}[k] + \bar{\mathbf{y}}_\mathrm{init}\delta[k],\\[2pt] p^\mathrm{d}[\mathbf{x},k] &= \mathbf{c}_1^{\mathrm{\scriptsize T}}(\mathbf{x})\bar{\mathbf{y}}^\mathrm{d}[k]. \end{align*}

1.2 Parameters

The proposed model has been derived and is evaluated in terms of normalized physical quantities with respect to a reference length $\rho$ and a reference time $\tau$. For numerical evaluation, the parameters in Table 1 are considered.

Parameter	Value	Normalized value
Radius $R_0$	$100\mu \text{m}$	$1$
Length $Z_0$	$1\text{mm}$	$10$
Flow velocity $v_0$	$50\mu {\text{m}}^{-1}$	$50$
TX/RX distance $d$	$100\mu \text{m}$	$1$
Diffusion coefficient $D$	$2.5\cdot {10}^{-12},\dots ,5\cdot {10}^{-9}{\text{m}}^2{\text{s}}^{-1}$	$2.5\cdot {10}^{-2},\dots ,50$
Reference length $\rho$	$R_0$
Reference time $\tau$	$1\cdot {10}^2\text{s}$

Table 1: Physical parameters for numerical evaluation

1.3 Characterization into regimes

The proposed model is classified with respect to existing analytical models. Therefore, the dispersion factor $\alpha$, is introduced [2, Eq. (20)] $$\begin{array}{r}\alpha =\frac{Dd}{v_{\mathrm{eff}}{R}_0^2}=\frac{2Dd}{v_0{R}_0^2},\end{array}$$ to distinguish between different regimes. Dispersion factor $\alpha$ is used to classify the considered scenarios into different regimes where the different regimes are realized by changing the diffusion coefficient $D$. For numerical evaluation, the diffusion coefficients $D$ and resulting $\alpha$ values in Table 2 are considered.

$D$ in $\text{m}{\text{s}}^{-2}$	$2.5\cdot {10}^{-12}$	$2.5\cdot {10}^{-11}$	$2.5\cdot {10}^{-10}$	$7.5\cdot {10}^{-10}$	$1.25\cdot {10}^{-9}$	$2.5\cdot {10}^{-9}$	$5\cdot {10}^{-9}$
$\alpha$	$1\cdot {10}^{-3}$	$1\cdot {10}^{-2}$	$0.1$	$0.3$	$0.5$	$1$	$2$

Table 2: Considered diffusion coefficients $D$ and dispersion factors $\alpha$

1.4 Validation

For validation, the proposed model is compared to the results from particle based simulations (PBS) and the well-known models for flow-dominant and dispersive regime, see e.g. [3, Eqs. (11), (16)]. The parameters of the PBS are further described in [1, Sec. V-C].

2. Uniform release of particles

Supplementary to Fig. 3 in [1, Sec. V-D], we present videos of the particle concentration $p(\mathbf{x},t)$ in the $y$-$z$-plane ($y=r\sin (\frac{\pi }{2})$). Different values of $\alpha$ are considered.

We note that the value $\alpha$ has been derived based on the normalized TX/RX distance $d=1$ (see Table 1). Therefore, when the particles pass $z=2$ in the videos, the regime becomes more dispersive as $\alpha$ increases, i.e., $\alpha \sim d$.

The videos for a specific $\alpha$ are selected by clicking on the corresponding button

$\alpha$	Eigenvalues	Release type
varied	$Q=6000$, $\mathbf{q}=[0,30,200]$	uniform

$\alpha =1\cdot {10}^{-3}$ $\alpha =1\cdot {10}^{-2}$ $\alpha =0.1$ $\alpha =0.3$ $\alpha =0.5$ $\alpha =1$ $\alpha =2$

3. Point release of particles

Supplementary to Fig. 5 in [1, Sec. V-E], we present videos of the particle concentration $p(\mathbf{x},t)$ in the $y$-$z$-plane ($y=r\sin (\frac{\pi }{2})$) for a point release of particles. The point releases are centered on the radial axis at three positions $r_e=0.25,0.5$, and $0.75$.

We note that the values $\alpha$ have been derived based on the normalized TX/RX distance $d=1$ (see Table 1). Therefore, when the particles pass $z=1$ in the videos, the regime becomes more dispersive and $\alpha$ increases, i.e., $\alpha \sim d$.

The videos for a specific release position $r_e$ are selected by clicking on the corresponding button

$\alpha$	Eigenvalues	Release type
$\alpha =1\cdot {10}^{-3}$	$Q=12.6\cdot {10}^4$, $\mathbf{q}=[20,30,200]$	point

$r_e=0.75$ $r_e=0.5$ $r_e=0.25$

$\alpha$	Eigenvalues	Release type
$\alpha =1\cdot {10}^{-2}$	$Q=12.6\cdot {10}^4$, $\mathbf{q}=[20,30,200]$	point

$r_e=0.75$ $r_e=0.5$ $r_e=0.25$

$\alpha$	Eigenvalues	Release type
$\alpha =0.1$	$Q=12.6\cdot {10}^4$, $\mathbf{q}=[20,30,200]$	point

$r_e=0.75$ $r_e=0.5$ $r_e=0.25$

$\alpha$	Eigenvalues	Release type
$\alpha =0.5$	$Q=12.6\cdot {10}^4$, $\mathbf{q}=[20,30,200]$	point

$r_e=0.75$ $r_e=0.5$ $r_e=0.25$

$\alpha$	Eigenvalues	Release type
$\alpha =1$	$Q=12.6\cdot {10}^4$, $\mathbf{q}=[20,30,200]$	point

$r_e=0.75$ $r_e=0.5$ $r_e=0.25$

In addition to Fig. 5 in [1, Sec. V-E] we present the concentration at the RX over time for $\alpha =1$ (not shown in the paper):

$\alpha$	Eigenvalues	Release type
$\alpha =2$	$Q=12.6\cdot {10}^4$, $\mathbf{q}=[20,30,200]$	point

$r_e=0.75$ $r_e=0.5$ $r_e=0.25$

In addition to Fig. 5 in [1, Sec. V-E] we present the concentration at the RX over time for $\alpha =2$ (not shown in the paper):

4. Analysis of accuracy

Supplementary to [1, Sec. VI], we present videos of the particle concentration $p(\mathbf{x},t)$ in the $y$-$z$-plane ($y=r\sin (\frac{\pi }{2})$).

4.1 Numerical reflections due to absorbing boundary

As described in [1, Sec. VI-B], the numerical evaluation of the proposed model leads to undesired effects due to the restriction of the cylinder in $z$-direction.

$\alpha$	Eigenvalues	Release type
$\alpha =1\cdot {10}^{-3}$	$Q=6000$, $\mathbf{q}=[0,30,200]$	uniform

Instead of leaving the cylinder at $z=Z_0$, particles are reflected and re-enter the cylinder at $z=0$. The re-entering process is accompanied by undesired reflections.

The occurring effect is well investigate in the modeling of wave-propagation in the free space and it can be avoided by the design of a perfectly matched layer. In future work, we want to overcome this effect by adopting the techniques presented in [4].

4.2 Undesired oscillations for small values of $\alpha$

As described in [1, Sec. VI-B], the numerical evaluation of the proposed model is limited in the flow dominant regime for very small values of $\alpha$. The amplitude, duration and shape of the uniform release of particles is perfectly captured at the RX position (see [1, Fig. 8c]), while undesired oscillations occur afterwards. The video shows the emergence of these oscillations.

$\alpha$	Eigenvalues	Release type
$\alpha =1\cdot {10}^{-5}$	$Q=6000$, $\mathbf{q}=[0,30,200]$	uniform

In the video it can be seen, that the laminar shape of the flow profile is not captured perfectly. Furthermore, it is superimposed by undesired oscillations.

These oscillations originate from the well-known Gibb’s phenomenon which generally occurs when discontinuous and peaky functions are approximated by a finite series of continuous sine or cosine functions. Therefore, for the problem at hand, the phenomenon occurs for peaky particle releases that do not change their shape during transmission in the flow based regime. For larger values of $\alpha$, the oscillations do not occur because the initially sharp peak is blurred by diffusion. Although, Gibb’s phenomenon cannot be avoided completely for scenarios with very small $\alpha$ and a peaky initial release, it can be suppressed by increasing $Q$ (see [1, Fig. 9c]).

4.3 Dependence of accuracy on number of eigenvalues $Q$

Supplementary to Figs. 8a and b in [1, Sec. VI-B] we present videos of the particle concentration $p(\mathbf{x},t)$ in the $y$-$z$-plane ($y=r\sin (\frac{\pi }{2})$) for different values of $Q$. The videos emphasize the dependency of the numerical evaluation of the proposed model on the number of eigenvalues $Q$.

The number of eigenvalues is [1, Eq. (67)] $$\begin{array}{rlr}& Q=(2N+1)\cdot M\cdot L& \mathbf{q}=[N,M,L].\end{array}$$

4.3.1 Uniform release of particles

The following videos correspond to the scenario in Fig. 8a, i.e., a uniform release of particles for different values of $Q$.

$\alpha$	Eigenvalues	Release type
$\alpha =0.1$	varied	uniform

The videos for a specific number of eigenvalues $Q$, defined by vector $\mathbf{q}$ are selected by clicking on the corresponding button

$\mathbf{q}=[0,30,200]$ $\mathbf{q}=[0,20,150]$ $\mathbf{q}=[0,10,50]$ $\mathbf{q}=[0,2,50]$ $\mathbf{q}=[0,2,10]$

• The videos show, that the proposed model perfectly captures the influence of diffusion and flow on the propagation of the particles for $Q=6000$ and $Q=3000$.
• For $Q=500$ it can be seen that for $t=0$ the initial uniform distribution is interfered by reflections that originate from Gibb’s phenomenon, i.e., the number of eigenvalues is too low to represent the initial distribution correctly in the volume. Furthermore, the video shows that these reflections attenuate very fast so that the received concentration is not interfered, see green curve in [1, Fig. 8a].
• For $Q=100$ and $Q=20$ the reflections become stronger and are not attenuating, instead they are propagating in the cylinder and interfere the concentration at the receiver, see black and orange curves in [1, Fig. 8a].

4.3.2 Point release of particles

The following videos correspond to the scenario in Fig. 8b, i.e., a point release of particles for different values of $Q$.

$\alpha$	Eigenvalues	Release type
$\alpha =1\cdot {10}^{-2}$	varied	point

The videos for a specific number of eigenvalues $Q$, defined by vector $\mathbf{q}$ are selected by clicking on the corresponding button

$\mathbf{q}=[20,30,200]$ $\mathbf{q}=[20,30,100]$ $\mathbf{q}=[20,10,100]$ $\mathbf{q}=[5,10,100]$ $\mathbf{q}=[2,10,100]$

• The videos show, that the proposed model perfectly captures the influence of diffusion and flow on the propagation of the particles for $Q=12.6\cdot {10}^4$ and $Q=6.3\cdot {10}^4$.
• For $Q=2.1\cdot {10}^4$ it can be seen that for $t=0$ the initial point distribution is interfered by reflections that originate from Gibb’s phenomenon, i.e., the number of eigenvalues is too low to represent the initial distribution correctly in the volume. Furthermore, the video shows that these reflections attenuate very fast so that the received concentration is not interfered, see green curve in [1, Fig. 8b].
• For $Q=5\cdot {10}^3$ and $Q=2\cdot {10}^3$ an interesting effect occurs for the point release that originates from lowering the number $N$ of Bessel functions. A mirror source is created symmetrically to the considered point release.

5. Implementation

In the following we provide further details on the implementation of the proposed model in Matlab. The complete code can be downloaded via GitHub.

Symmetry of Bessel functions

As mentioned in [1, Secs. VI-A and V-E], the number of eigenvalues $Q$ can be reduced exploiting the symmetry of Bessel functions.

Starting with the particle concentration $P$ formulated in terms of [1, Eq. (26)] $$\begin{array}{r}P(\mathbf{x},s)=\sum _{\mu =0}^{Q-1}\frac{1}{N_{\mu }}\overline{Y}(\mu ,s)K_1(\mathbf{x},\mu ).\end{array}$$ Exploiting the purpose of $\mu$ as index tupel $(n,m,\nu )\to \mu$, the sum can be reformulated in terms of three sums $$\begin{array}{rl}P(\mathbf{x},s)& =\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}\sum _{m=0}^{M-1}\sum _{\nu =0}^{L-1}\frac{1}{N_{n,m,\nu }}\overline{Y}(n,m,\nu ,s)K_1(\mathbf{x},n,m,\nu ).\end{array}$$ Furthermore, the sum over $n$ can be divided into the positive ($n>0$) and negative ($n<0$) values and a term for $n=0$ yielding $$\begin{array}{rl}P(\mathbf{x},s)& =\frac{1}{N_{0,m,\nu }}\overline{Y}(0,m,\nu ,s)K_1(\mathbf{x},n,m,\nu )\\ & +\sum _{n=-\mathrm{\infty }}^{-1}\sum _{m=0}^{M-1}\sum _{\nu =0}^{L-1}\frac{1}{N_{n,m,\nu }}\overline{Y}(n,m,\nu ,s)K_1(\mathbf{x},n,m,\nu )\\ & +\sum _{n=1}^{\mathrm{\infty }}\sum _{m=0}^{M-1}\sum _{\nu =0}^{L-1}\frac{1}{N_{n,m,\nu }}\overline{Y}(n,m,\nu ,s)K_1(\mathbf{x},n,m,\nu )\\ & =P^0(\mathbf{x},s)+P^{-}(\mathbf{x},s)+P^{+}(\mathbf{x},s),\end{array}$$ where the sum in $P^{-}(\mathbf{x},s)$ can be expressed in terms of the positive $n$ $$\begin{array}{rl}{P}^{-}(\mathbf{x},s)& =\sum _{n=-\mathrm{\infty }}^{-1}\sum _{m=0}^{M-1}\sum _{\nu =0}^{L-1}\frac{1}{N_{n,m,\nu }}\overline{Y}(n,m,\nu ,s)K_1(\mathbf{x},n,m,\nu )\\ & =\sum _{n=1}^{\mathrm{\infty }}\sum _{m=0}^{M-1}\sum _{\nu =0}^{L-1}\frac{1}{N_{-n,m,\nu }}\overline{Y}(-n,m,\nu ,s)K_1(\mathbf{x},-n,m,\nu )\end{array}$$ Exploiting the symmetry of Bessel functions, i.e., $J_{-n}(x)=(-1{)}^n{J}_n(x)$, the individual terms for $n<0$ in $P^{-}(\mathbf{x},s)$ can be related to the terms for $n>0$ in $P^{+}(\mathbf{x},s)$ by the relations $$ \begin{align*} N_{-n,m,\nu} &= N_{n,m,\nu} \\[8pt] \bar{Y}(-n,m,\nu,s) &= (-1)^n\bar{Y}^(n,m,\nu,s)\\[8pt] K_1(\mathbf{x},-n,m,\nu) &= (-1)^n K_1^(\mathbf{x},n,m,\nu) \end{align*} $$Finally,\$P(\mathbf{x},s)\$canbeexpressedas$$ \begin{align*} P(\mathbf{x},s) &= P^0(\mathbf{x},s) + P^{-}(\mathbf{x},s) + P^{+}(\mathbf{x},s) \\ &= P^0(\mathbf{x},s) + P^{+}(\mathbf{x},s) + P^{+}(\mathbf{x},s) \\ &= P^0(\mathbf{x},s) + 2\mathrm{Re}\{P^{+}(\mathbf{x},s)\} \end{align} $$ By this reformulation it is sufficient to evaluate the proposed model for the values $n\ge 0$ and the negative ones follow from the positive ones. This modifies the derivation of $Q$ in [1, Eq. (67)] $$\begin{array}{rlrlr}& Q=(2N+1)\cdot M\cdot L& \to & & Q=(N+1)\cdot M\cdot L\end{array}$$

Decomposition

As mentioned in [1, Sec. VI-A] the block-diagonal structure of the state matrix ${\mathcal{A}}_{\mathrm{c}}^{\mathrm{d}}$ can be exploited to speed up the numerical evaluation. The block structuring is shown in Fig. 2. The matrix consists of $N+1$ blocks of size $(M\cdot L)\times (M\cdot L)$.

Fig. 2: Structuring of feedback matrix ${\mathcal{A}}_{\mathrm{c}}^{\mathrm{d}}$

The block structuring in Fig. 2 depends strongly on how the eigenvalues $s_{\mu }=s_{n,m,\nu }$ are counted. To illustrate the counting of eigenvalues, a small example is shown in Table 3 for $M=2$ and $L=2$.

$\mu$	$n$	$m$	$\nu$	$s_{\mu }=s_{n,m,\nu }=-D,(k_{n,m}^2+{\lambda }_{\nu }^2)$
0	0	0	0	$s_0=s_{0,0,0}=-D,(k_{0,0}^2+{\lambda }_0^2)$
1	0	0	1	$s_1=s_{0,0,1}=-D,(k_{0,0}^2+{\lambda }_1^2)$
2	0	1	0	$s_2=s_{0,1,0}=-D,(k_{0,1}^2+{\lambda }_0^2)$
3	0	1	1	$s_3=s_{0,1,1}=-D,(k_{0,1}^2+{\lambda }_1^2)$
4	1	0	0	$s_4=s_{1,0,0}=-D,(k_{1,0}^2+{\lambda }_0^2)$
5	1	0	1	$s_5=s_{1,0,1}=-D,(k_{1,0}^2+{\lambda }_1^2)$
$⋮$

Table 3: Counting of eigenvalues for $M=2$, $L=2$

Motivated by the block-structuring, the vectors and matrices in the state equation [1, Eq. (65)] can be rearranged as follows $$ \begin{align*} &\bar{\mathbf{y}}^\mathrm{d}[k] = \begin{bmatrix} \bar{\mathbf{y}}0^\mathrm{d}[k]\\ \vdots \\ \bar{\mathbf{y}}n^\mathrm{d}[k]\\ \vdots \\ \bar{\mathbf{y}}N^\mathrm{d}[k]\end{bmatrix}, &\bar{\mathbf{f}}\mathrm{e}^\mathrm{d}[k] = \begin{bmatrix} \bar{\mathbf{f}}{\mathrm{e},0}^\mathrm{d}[k]\\ \vdots \\ \bar{\mathbf{f}}{\mathrm{e},n}^\mathrm{d}[k]\\ \vdots \\ \bar{\mathbf{f}}{\mathrm{e},N}^\mathrm{d}[k]\end{bmatrix}, &&\bar{\mathbf{y}}\mathrm{init}^\mathrm{d} = \begin{bmatrix} \bar{\mathbf{y}}{\mathrm{init},0}^\mathrm{d}\\ \vdots \\ \bar{\mathbf{y}}{\mathrm{init},n}^\mathrm{d}\\ \vdots \\ \bar{\mathbf{y}}{\mathrm{init},N}^\mathrm{d}\end{bmatrix}, \end{align*} $$$$ \begin{align*} \mathcal{A}\mathrm{c}^\mathrm{d}(v_0)&= \begin{bmatrix}\mathcal{A}{\mathrm{c},0}^\mathrm{d}(v_0) & & 0 & & 0\\ 0 & \ddots & & & \\ \vdots & & \mathcal{A}{\mathrm{c},n}^\mathrm{d}(v_0) & & \\ & & &\ddots & \\ 0 & & & & \mathcal{A}_{\mathrm{c},N}^\mathrm{d}(v_0)\end{bmatrix}, \end{align*} $$ where the involved block matrices ${\mathcal{A}}_{\mathrm{c},n}^{\mathrm{d}}$ are of size $(M\cdot L)\times (M\cdot L)$ and the vectors of size $(M\cdot L)\times 1$. Codeblock 1 shows how this block structuring affects the implementation. Instead of a direct evaluation of state equation [1, Eq. (65)], the individual blocks are evaluated separately.

Codeblock 1: Block-wise evaluation of the state equation.

Each block ${\mathcal{A}}_{\mathrm{c},n}^{\mathrm{d}}$ of ${\mathcal{A}}_{\mathrm{c}}^{\mathrm{d}}$ may still be fully occupied. Therefore, the state matrix ${\mathcal{A}}_{\mathrm{c},n}$ is decomposed into a diagonal matrix ${\mathrm{\Lambda }}_n$ of eigenvalues and a matrix $\mathbf{U}n$ofeigenvectors$$ \begin{align*} \mathcal{A}{\mathrm{c},n}(v_0) &= \mathbf{U}_n\mathbf{\Lambda}n\mathbf{U}n^{-1}\\ \mathcal{A}{\mathrm{c},n}^{\mathrm{d}}(v_0) &= \mathrm{e}^{\mathcal{A}{\mathrm{c},n} (v_0)T } = \mathbf{U}_n\mathrm{e}^{\mathbf{\Lambda}_n T}\mathbf{U}_n^{-1} \end{align*} $$ Inserting this relation into the state equation [1, Eq. (65)] and multiplying by ${\mathbf{U}}^{-1}$ from the left hand side leads to a modified state equation $$ \begin{align*} \mathbf{U}_n^{-1}\bar{\mathbf{y}}_n^\mathrm{d}[k] &= \mathrm{e}^{\mathbf{\Lambda}_n T}\mathbf{U}_n^{-1}\bar{\mathbf{y}}_n^\mathrm{d}[k - 1] + \mathbf{U}n^{-1}\bar{\mathbf{f}}^\mathrm{d}{\mathrm{e},n}[k] + \mathbf{U}n^{-1}\bar{\mathbf{y}}{\mathrm{init},n}\delta[k], \\ \bar{\mathbf{u}}_n^\mathrm{d}[k] &= \mathrm{e}^{\mathbf{\Lambda}_n T}\bar{\mathbf{u}}_n^\mathrm{d}[k - 1] + \mathbf{U}n^{-1}\bar{\mathbf{f}}^\mathrm{d}{\mathrm{e},n}[k] + \mathbf{U}n^{-1}\bar{\mathbf{y}}{\mathrm{init},n}\delta[k], \end{align*} $$whereauxiliarysystemstatesaredefinedas$$ $$\begin{array}{r}{\overline{\mathbf{u}}}_n^{\mathrm{d}}={\mathbf{U}}_n^{-1}{\overline{\mathbf{y}}}_n^{\mathrm{d}}.\end{array}$$ $$

The modified state equation can be evaluated more efficiently than the unmodified version [1, Eq. (65)] as the multiplication by a fully occupied matrix is replaced by the multiplication by a diagonal matrix. Therefore, it can be evaluated very efficiently by first-order filters (see Codeblock 2).

Codeblock 2: Efficient evaluation of the modified state equation in terms of first-order filters

6. Analysis of runtime and complexity

In the following we present a short analysis of the runtime and complexity of the Matlab implementation of the proposed model.

6.1 Complexity

As described in [1, Sec. VI-A], many calculations can be done in advance. Therefore, only the complexity of the numerical evaluation of the state equation in Codeblock 2 is presented.

The filter operation

which evaluates one line of the modified state equation performs one addition and one multiplication. The complexity depends on the temporal length of the input signal or rather on the observation time of the system. Denoting this time by $T_{\mathrm{dur}}$, the complexity is $\mathcal{O}(T_{\mathrm{dur}})$. This operation has to be performed $(M\cdot L)$ times for one function call (see Codeblock 2)

Transforming the auxiliary system states $\overline{\mathbf{u}}$ into the system states $\overline{\mathbf{y}}$

exhibits a complexity of $\mathcal{O}((M\cdot L{)}^2\cdot T_{\mathrm{dur}})$

6.2 Runtime

The analysis of runtime strongly depends on the efficiency of the implementation, the programming language and the machine on that it is evaluated. In the following, the runtime of the Matlab implementation of the proposed model is presented. Particularly, the runtime of the evaluation of the state equation in Codeblock 1 is shown in Table 3 for different values of $Q$.

The numerical evaluation has been performed with the provided Matlab implementation on a Macbook-Pro with a $2.8$ GHz Quad-Core Intel Core i7 processor and $16$GB RAM.

$Q$	$\mathbf{q}$	state equation
$500$	$[0,10,50]$	$(N+1)\times 0.05$s
$3000$	$[0,20,150]$	$(N+1)\times 0.8$s
$6000$	$[0,30,200]$	$(N+1)\times 3.5$s
$2.1\cdot {10}^4$	$[20,10,100]$	$(N+1)\times 0.1$s
$6.3\cdot {10}^4$	$[20,30,100]$	$(N+1)\times 1$s
$12.5\cdot {10}^4$	$[20,30,200]$	$(N+1)\times 5$s

Table 4: Runtimes for different values of $Q$

7. References

[1] M. Schäfer, W. Wicke, L. Brand, R. Rabenstein and R. Schober, “Transfer Function Models for Diffusion and Laminar Flow in Cylindrical Systems”, submitted to IEEE Trans. Mol. Biol. Multi-Scale Commun., 2020, [online]: https://arxiv.org/abs/2007.01799

[2] V. Jamali, A. Ahmadzadeh, W. Wicke, A. Noel, and R. Schober, “Channel Modeling for Diffusive Molecular Communication—A Tutorial Review,” Proceedings of the IEEE, vol. 107, no. 7, pp. 1256–1301, 2019.

[3] W. Wicke, T. Schwering, A. Ahmadzadeh, V. Jamali, A. Noel, and R. Schober, “Modeling Duct Flow for Molecular Communication,” in Proc. IEEE Global Commun. Conf. (GLOBECOM), 2018, pp. 206–212.

[4] 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