A new 3D mass diffusion–reaction model in the neuromuscular junction

Abstract

A three-dimensional model of the reaction-diffusion processes of a neurotransmitter and its ligand receptor in a disk shaped volume is proposed which represents the transmission process of acetylcholine in the synaptic cleft in the neuromuscular junction. The behavior of the reaction-diffusion system is described by a three-dimensional diffusion equation with nonlinear reaction terms due to the rate processes of acetylcholine with the receptor. A new stable and accurate numerical method is used to solve the equations with Neumann boundaries in cylindrical coordinates. The simulation analysis agrees with experimental measurements of end-plate current, and agrees well with the results of the conformational state of the acetylcholine receptor as a function of time and acetylcholine concentration of earlier investigations with a smaller error compared to experiments. Asymmetric emission of acetylcholine in the synaptic cleft and the subsequent effects on open receptor population is simulated. Sensitivity of the open receptor dynamics to the changes in the diffusion parameters and neuromuscular junction volume is investigated. The effects of anisotropic diffusion and non-symmetric emission of transmitter at the presynaptic membrane is simulated.

Appendix

1.1 Finite difference scheme for Neumann boundary condition

To obtain a new finite difference scheme without discretizing the Neumann boundary condition, Eq. (9a), we first design a mesh, where the distance between the actual left boundary and z ₁ is assumed to be θ ₁Δz, and the distance between the actual right boundary and z _K is θ ₂Δz, as shown in Fig. 4(c). For simplicity, we denote (A)(r, Φ, z, t) and Δz as A(z, t) and h, respectively, and then express the finite difference approximation of $\frac{\partial ^{2}A(z,t)}{ \partial z^{2}}$ at z ₁, which is the grid point next to the left boundary, as follows:

$$ \begin{array}{rll} b\frac{\partial ^{2}A(z_{1},t)}{\partial z^{2}}&=&\frac{a}{h^{2} }\big[A(z_{2},t)-A(z_{1},t)\big]\\ &&-\,\frac{1}{h}\frac{\partial A\left( z_{1}-\theta _{1}h,t\right) }{\partial z}, \label{eqA1} \end{array}$$

(18)

where a, b, θ ₁ are constants to be determined. If Eq. (18) is rewritten as follows:

$$\begin{array}{rl} \label{eqA1'} &b\frac{\partial ^{2}A(z_{1},t)}{\partial z^{2}}+\frac{1}{h} \frac{\partial A\left( z_{1}-\theta _{1}h,t\right) }{\partial z}\\ &=\frac{a}{ h^{2}}\left[A(z_{2},t)-A(z_{1},t)\right], \end{array} $$

(18')

one may see that the above equation is an improvement of the combined compact finite difference method (where the first and second-order derivatives are included (Chu and Fan 1998; Zhao et al. 2007, 2008) by introducing the parameter θ ₁ in order to raise the order of accuracy. The first-order derivative is kept in Eq. (18) so that the Neumann boundary condition can be applied directly without discretizing. Expanding each term of Eq. (18) into Taylor series at z ₁, we obtain the right-hand-side (RHS) result of Eq. (18) as follows:

$$ \begin{array}{rll} \text{RHS} &=&\frac{a}{h^{2}}\left[hA_{z}(z_{1},t)+\frac{h^{2}}{2}A_{zz}(z_{1},t)+ \frac{h^{3}}{6}A_{z^{3}}(z_{1},t)\right] \\ &&-\frac{1}{h}\!\left[\!A_{z}(z_{1},t){\kern-1pt} -{\kern-1pt} \theta _{1}hA_{zz}(z_{1},t){\kern-1pt} +{\kern-2pt} \,\frac{\theta _{1}^{2}h^{2}}{2}A_{z^{3}}(z_{1},t)\!\right] \\ &&+\,O\left(h^{2}\right) \\ &=&\frac{1}{h}\left[ a-1\right] A_{z}(z_{1},t) \\ &&+\,\left[\frac{a}{2}+\theta _{1}\right]A_{zz}(z_{1},t)+\frac{h}{2}\left[\frac{a}{3}-\theta _{1}^{2} \right]A_{z^{3}}(z_{1},t) \\ &&+O\left(h^{2}\right). \label{eqA2} \end{array}$$

(19)

Matching both sides gives

$$ a=1,\text{}b=\frac{1}{2}+\frac{\sqrt{3}}{3},\text{}\theta _{1}=\frac{\sqrt{3}}{3}. \label{eqA3} $$

(20)

Thus, substituting the values of a, b, θ ₁ in Eq. (20) into Eq. (18) and dropping the truncation error O(h ²), we obtain a second-order finite difference approximation for $\frac{\partial ^{2}A}{\partial z^{2}}$ at z ₁ as

$$ \begin{array}{rll} \frac{\partial ^{2}A(z_{1},t)}{\partial z^{2}}&\approx& \frac{a }{bh^{2}}\left[A(z_{2},t)-A(z_{1},t)\right]\\ &&-\,\frac{1}{bh}\frac{\partial A\left( z_{1}-\theta _{1}h,t\right) }{\partial z}. \label{eqA4} \end{array}$$

(21)

Symmetrically, we can express the finite difference approximation of $\frac{ \partial ^{2}A(z,t)}{\partial z^{2}}$ at z _K , which is the grid point next to the right boundary, as

$$ \begin{array}{rll} b^{*}\frac{\partial ^{2}A(z_{K},t)}{\partial z^{2}}&=&\frac{1}{h}\frac{ \partial A\left( z_{K}+\theta _{2}h,t\right) }{\partial z}\\ &&-\,\frac{a^{*}}{h^{2} }\left[A(z_{K},t)-A(z_{K-1},t)\right], \label{eqA5} \end{array}$$

(22)

where $a^{*},b^{*},\theta _{2}$ are constants to be determined. Again, matching both sides in Taylor series gives

$$ a^{*}=1,\text{}b^{*}=\frac{1}{2}+\frac{\sqrt{3}}{3},\text{ }\theta _{2}=\frac{\sqrt{3}}{3}, \label{eqA6} $$

(23)

and hence a second-order finite difference approximation at z _K for the right boundary can be obtained as

$$ \begin{array}{rll} \frac{\partial ^{2}A(z_{K},t)}{\partial z^{2}}&\approx& \frac{1}{b^{*}h}\frac{ \partial A\left( z_{K}+\theta _{2}h,t\right) }{\partial z}\\ &&-\,\frac{a^{*}}{ b^{*}h^{2}}\left[A(z_{K},t)-A(z_{K-1},t)\right]. \label{eqA7} \end{array}$$

(24)

If the number of interior grid points K is given, then the grid size and the coordinates of the grid points can be determined as follows:

$$ \begin{array}{rll} h&=&\frac{L}{K+\theta _{1}+\theta _{2}-1},\\ z_{k}&=&(k-1+\theta _{1})h,\text{}k=1,\cdots ,K. \label{eqA8} \end{array}$$

(25)

Based on the Neumann boundary condition, Eqs. (9a), (21) and (24) can be simplified to

$$ \frac{\partial ^{2}A(z_{1},t)}{\partial z^{2}}\approx \frac{a }{bh^{2}}\left[A(z_{2},t)-A(z_{1},t)\right], \label{eqA9a} $$

(26a)

$$ \frac{\partial ^{2}A(z_{K},t)}{\partial z^{2}}\approx -\frac{a^{*}}{ b^{*}h^{2}}\left[A(z_{K},t)-A(z_{K-1},t)\right]. \label{eqA9b} $$

(26b)

1.2 Stability of the finite difference scheme

To analyze the stability of the finite difference scheme, Eqs. (13a)–(13c), with the initial and boundary conditions, Eqs. (14a)–(16b), we first define following finite difference operators:

$$ \begin{array}{rll} P_{r}\left[A_{i,j,k}^{n}\right]&\equiv& r_{i+\frac{1}{2}}\frac{ A_{i+1,j,k}^{n}{\kern-1pt} -{\kern-1pt} A_{i,j,k}^{n}}{\left( \Delta r\right) ^{2}}{\kern-.5pt} \!-\!r_{i-\frac{1}{2}} \frac{A_{i,j,k}^{n}-A_{i-1,j,k}^{n}}{\left( \Delta r\right) ^{2}},\\[4pt] P_{\phi}\left[A_{i,j,k}^{n}\right]&\equiv& \frac{ A_{i,j+1,k}^{n}-2A_{i,j,k}^{n}+A_{i,j-1,k}^{n}}{\left( \Delta \phi \right) ^{2}},\\[4pt] P_{z}\left[A_{i,j,k}^{n}\right]&\equiv& \frac{ A_{i,j,k+1}^{n}-2A_{i,j,k}^{n}+A_{i,j,k-1}^{n}}{\left( \Delta z\right) ^{2}},\\[4pt] \nabla _{\bar{r}}A_{i,j,k}^{n}&\equiv& \frac{A_{i,j,k}^{n}-A_{i-1,j,k}^{n}}{ \Delta r},\\[4pt] \nabla_{\bar{\phi}}A_{i,j,k}^{n}&\equiv& \frac{ A_{i,j,k}^{n}-A_{i,j-1,k}^{n}}{\Delta \phi },\\[4pt] \nabla _{\bar{z}}A_{i,j,k}^{n}&\equiv& \frac{A_{i,j,k}^{n}-A_{i,j,k-1}^{n}}{ \Delta z},\\[4pt] W_{t}\left[A_{i,j,k}^{n}\right]&\equiv&\frac{ A_{i,j,k}^{n+1}+A_{i,j,k}^{n}}{2}. \end{array} $$

We assume that the values of (R), (AR), (A ₂ R), $(A_{2}R^{\rm open})$ at time step n + 1 have already been known since they are calculated ahead based on the Runge–Kutta method. As such, Eqs. (13a)–(13c) can be simplified to, at interior points, k = 2, ⋯ ,K − 1,

$$ \begin{array}{rll} \frac{A_{i,j,k}^{n+1}\!-\! A_{i,j,k}^{n}}{\Delta t}\! &=& \! \frac{D_{r}}{r_{i}} P_{r}\left\{\! W_{t}\big[A_{i,j,k}^{n}\big]\!\right\} \!+\! \frac{D_{\phi }}{r_{i}^{2}} P_{\phi }\left\{ W_{t}\big[A_{i,j,k}^{n}\big]\right\} \\ [4pt] &&+\,D_{z}P_{z}\left\{ W_{t}\big[A_{i,j,k}^{n}\big]\right\} \\ [4pt] &&-\,C_{i,j,k}^{\,n+\frac{1}{2}}\left(\frac{A_{i,j,k}^{n+1}+A_{i,j,k}^{n}}{2} \right){\kern-3pt} +{\kern-3pt} f_{i,j,k}^{\,n+\frac{1}{2}}, \label{eqA10a} \end{array}$$

(27a)

where 1 ≤ i ≤ I − 1, 0 ≤ j ≤ J − 1; at the location z ₁,

$$\begin{array}{rll} \frac{A_{i,j,1}^{n+1}\!-\!A_{i,j,1}^{n}}{\Delta t} \!&=&\! \frac{D_{r}}{r_{i}} P_{r}\left\{ W_{t}\big[A_{i,j,1}^{n}\big]\right\} \!+\! \frac{D_{\phi }}{r_{i}^{2}} P_{\phi }\left\{ W_{t}\big[A_{i,j,1}^{n}\big]\right\} \\ &&+\,\frac{a}{b}\frac{D_{z}}{ \Delta z}\nabla _{\bar{z}}\left\{ W_{t}\big[A_{i,j,1}^{n}\big]\right\} \\ &&-\,C_{i,j,1}^{n+\frac{1}{2}}\left(\frac{A_{i,j,1}^{n+1}+A_{i,j,1}^{n}}{2} \right)+f_{i,j,1}^{n+\frac{1}{2}}, \label{eqA10b} \end{array}$$

(27b)

where $a=1,b=\frac{1}{2}+\frac{\sqrt{3}}{3}$ and 1 ≤ i ≤ I − 1, 0 ≤ j ≤ J − 1; and at the location z _K ,

$$ \begin{array}{rll} \frac{A_{i,j,K}^{n+1}\!-\!A_{i,j,K}^{n}}{\Delta t} \!&=&\! \frac{D_{r}}{r_{i}} P_{r}\!\left\{\! W_{t}\big[\!A_{i,j,K}^{n}\big]\!\right\} \!+\! \frac{D_{\phi }}{r_{i}^{2}} P_{\phi }\!\left\{\! W_{t}\big[A_{i,j,K}^{n}\big]\!\right\} \\ &&-\,\frac{a^{*}}{b^{*}}\frac{D_{z} }{\Delta z}\nabla _{\bar{z}}\left\{ W_{t}\big[A_{i,j,K}^{n}\big]\right\} \\ &&-C_{i,j,K}^{n+\frac{1}{2}}\left(\frac{A_{i,j,K}^{n+1}+A_{i,j,K}^{n}}{2} \right)+f_{i,j,K}^{n+\frac{1}{2}} \\ \label{eqA10c} \end{array}$$

(27c)

where $a^{*}=1,b^{*}=\frac{1}{2}+\frac{\sqrt{3}}{3}$ and 1 ≤ i ≤ I − 1, 0 ≤ j ≤ J − 1. Here, $C_{i,j,k}^{n+\frac{1}{2}}$ and $f_{i,j,k}^{n+ \frac{1}{2}},$ which are from the values of species (R, AR, A ₂ R, $A_{2}R^{\rm open})$, are considered to be a positive coefficient and a source term, respectively.

Theorem 1

The finite difference scheme, Eqs. (27a)–(27c), with the initial and boundary conditions, Eqs. (14a)–(14b), is unconditionally stable with respect to the initial condition and source term.

Proof

Assume that A ₁ and A ₂ are two solutions obtained based on Eqs. (27a)–(27c) with same boundary conditions but different initial conditions and source terms f ₁ and f ₂, respectively. Letting A = A ₁ − A ₂ and f = f ₁ − f ₂, then A and f satisfy Eqs. (27a)–(27c) with the boundary condition

$$ \begin{array}{rll} A_{i,J,k}^{n}&=&A_{i,0,k}^{n},A_{i,-1,k}^{n}=A_{i,J-1,k}^{n},\\ A_{I,j,k}^{n}&=&0,A_{0,j,k}^{n}=\frac{\Delta \phi }{2\pi }\displaystyle\sum\limits_{j=0}^{J-1}A_{1,j,k}^{n}, \label{eqA11} \end{array}$$

(28)

for 0 ≤ j ≤ J − 1, 1 ≤ k ≤ K.

We then multiply Eq. (27a) by 2r _i ΔzΔt $[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}]$ for interior points, k = 2, ⋯ ,K − 1; multiply Eq. (27b) by $2\frac{b}{a}r_{i}\Delta z\Delta t$ $[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}]$; multiply Eq. (27c) by $2\frac{b^{*}}{a^{*}}$ r _i ΔzΔt $[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}]$; and then sum them over k, where 1 ≤ k ≤ K. This gives

$$ \begin{array}{rll} &&\Delta zr_{i}\left\{\!\frac{b}{a} \!\left(\!\left[\!A_{i,j,1}^{n+1}\right]^{2}\!\!-\!\left[\!A_{i,j,1}^{n}\right]^{2}\right)\!+\!\sum \limits_{k=2}^{K-1}\!\left(\!\left[\!A_{i,j,k}^{n+1}\right]^{2}\!\!-\!\left[\!A_{i,j,k}^{n}\right]^{2}\right)\right.\\ &&{\kern24pt} \left.+\,\frac{b^{*}}{ a^{*}}\left(\left[\!A_{i,j,K}^{n+1}\right]^{2}\!-\!\left[\!A_{i,j,K}^{n}\right]^{2}\right)\right\} \\ &&=2\Delta z\Delta tD_{r}\left\{\frac{b}{a}P_{r}\left\{ W_{t}\big[A_{i,j,1}^{n}\big]\right\} \left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right]\right. \\ &&{\kern6pt} +\sum\limits_{k=2}^{K-1}P_{r}\left\{ W_{t}\big[A_{i,j,k}^{n}\big]\right\} \left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right] \\ &&{\kern6pt} +\left.\frac{b^{*}}{a^{*}}P_{r}\left\{ W_{t}\big[A_{i,j,K}^{n}\big]\right\} \left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\right\} \\ &&{\kern6pt} +\,2\Delta z\Delta t\frac{D_{\phi }}{r_{i}}\left\{\frac{b}{a}P_{\phi }\left\{ W_{t}\left[A_{i,j,1}^{n}\right]\right\} \left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right]\right.\\ &&{\kern6pt}+\left.\sum\limits_{k=2}^{K-1}P_{\phi }\left\{ W_{t}[A_{i,j,k}^{n}]\right\} \left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right]\right. \\ &&{\kern6pt}+\left.\frac{b^{*}}{a^{*}}P_{\phi }\left\{ W_{t}\left[A_{i,j,K}^{n}\right]\right\} \left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\right\} \\ &&{\kern6pt} +\,\Delta tr_{i}D_{z}\left\{\vphantom{\sum\limits_{k=2}^{K-1}}\nabla _{\bar{z} }\left[A_{i,j,2}^{n+1}+A_{i,j,2}^{n}\right]\left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right]\right.\\ &&{\kern52pt} \!+\,2\Delta z\sum\limits_{k=2}^{K-1}P_{z}\!\left\{ W_{t}\left[A_{i,j,k}^{n}\right]\right\}\left[A_{i,j,k}^{n+1}\!+\!A_{i,j,k}^{n}\right] \!\!\!\!\!\!\!\!\! \\ &&{\kern52pt} -\left.\nabla _{\bar{z} }\left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\vphantom{\sum\limits_{k=2}^{K-1}}\right\} \\ &&{\kern6pt} -\,\Delta z\Delta tr_{i}\left\{\vphantom{\sum\limits_{k=2}^{K-1}}\frac{b}{a}C_{i,j,1}^{n+\frac{1}{2} }\left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right]^{2}\right.\\ &&{\kern52pt} +\,\sum\limits_{k=2}^{K-1}C_{i,j,k}^{n+ \frac{1}{2}}\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right]^{2}\\ &&{\kern52pt} +\left.\frac{b^{*}}{a^{*}} C_{i,j,K}^{n+\frac{1}{2}}\left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]^{2}\vphantom{\sum\limits_{k=2}^{K-1}}\right\} \\ &&{\kern6pt} +\,2\Delta z\Delta tr_{i}\left\{\vphantom{\sum\limits_{k=2}^{K-1}}\frac{b}{a}f_{i,j,1}^{n+\frac{1}{2} }\left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right]\right.\\ &&{\kern58pt} +\,\sum\limits_{k=2}^{K-1}f_{i,j,k}^{n+\frac{1}{2}}\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right] \\ &&{\kern58pt} +\left.\frac{b^{*}}{a^{*}}f_{i,j,K}^{n+\frac{1}{2} }\left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\vphantom{\sum\limits_{k=2}^{K-1}}\right\}. \label{eqA12} \end{array}$$

(29)

Since the third term on the RHS of Eq. (29) can be simplified to

$$ \begin{array}{rll} &&{\kern-6pt} 2\Delta z\sum\limits_{k=2}^{K-1}P_{z}\left\{ W_{t}\left[A_{i,j,k}^{n}\right]\right\} \left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right]\\ &&{\kern12pt}+\nabla_{\bar{z}}\left[A_{i,j,2}^{n+1}+A_{i,j,2}^{n}\right]\left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right] \\ &&{\kern12pt}-\nabla _{\bar{z}}\left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right] \\ &&=\sum\limits_{k=2}^{K-1}\left\{\nabla _{\bar{z} }\left[A_{i,j,k+1}^{n+1}+A_{i,j,k+1}^{n}\right]-\nabla _{\bar{z}}\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right]\right\}\\ &&{\kern6pt} \times\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right] +\nabla _{\bar{z} }\left[A_{i,j,2}^{n+1}+A_{i,j,2}^{n}\right]\\ [4pt] &&{\kern6pt}\times \left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right]-\nabla _{\bar{z}}\left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\\ [5pt] &&{\kern6pt} \times\left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right] \\ [5pt] &&=\sum\limits_{k=3}^{K}\nabla _{\bar{z} }\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right]\left[A_{i,j,k-1}^{n+1}+A_{i,j,k-1}^{n}\right]\vphantom{\sum\limits_{k=2}^{K-1}}\\ [4pt] &&{\kern15pt} -\sum\limits_{k=2}^{K-1}\nabla _{\bar{z}}\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right]\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right] \\ [5pt] &&{\kern6pt} +\,\nabla _{\bar{z} }\left[A_{i,j,2}^{n+1}+A_{i,j,2}^{n}\right]\left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right]\\ [5pt] &&{\kern6pt} -\,\nabla _{\bar{z}}\left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right] \\ [5pt] &&=-\Delta z\sum\limits_{k=2}^{K}\nabla _{\bar{z} }\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right]^{2}, \label{eqA13} \end{array}$$

(30)

and $C_{i,j,k}^{n+\frac{1}{2}}[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}]^{2}\geq 0$, one may drop the third and fourth terms on the RHS of Eq. (29) and simplify it to

$$ \begin{array}{rll} &&{\kern-6pt} 2\Delta zr_{i}\left\{\vphantom{\sum\limits_{k=2}^{K-1}}\frac{b}{a} \left(\left[A_{i,j,1}^{n+1}\right]^{2}-\left[A_{i,j,1}^{n}\right]^{2}\right)\right.\\ &&{\kern24pt} +\,\sum\limits_{k=2}^{K-1}\left(\left[A_{i,j,k}^{n+1}\right]^{2}-\left[A_{i,j,k}^{n}\right]^{2}\right)\\ [4pt] &&{\kern24pt} +\left.\frac{b^{*}}{a^{*}}\left(\left[A_{i,j,K}^{n+1}\right]^{2}-\left[A_{i,j,K}^{n}\right]^{2}\right)\vphantom{\sum\limits_{k=2}^{K-1}}\right\} \\ [4pt] &&\leq 2\Delta z\Delta tD_{r}\left\{\vphantom{\sum\limits_{k=2}^{K-1}P_{r}}\frac{b}{a}P_{r}\left\{ W_{t}[A_{i,j,1}^{n}]\right\} \left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right]\right.\\ &&{\kern56pt} +\,\sum\limits_{k=2}^{K-1}P_{r}\left\{ W_{t}\left[A_{i,j,k}^{n}\right]\right\} \left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right] \\ &&{\kern56pt} +\left.\frac{b^{*}}{a^{*}}P_{r}\left\{ W_{t}[A_{i,j,K}^{n}]\right\} \left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\vphantom{\sum\limits_{k=2}^{K-1}P_{r}}\right\} \\ &&{\kern6pt}+2\Delta z\Delta t\frac{D_{\phi }}{r_{i}}\left\{\vphantom{\sum\limits_{k=2}^{K-1}}\frac{b}{a}P_{\phi }\left\{W_{t}\left[A_{i,j,1}^{n}\right]\right\} \left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right]\right.\\ &&{\kern66pt} +\,\sum\limits_{k=2}^{K-1}P_{\phi }\left\{ W_{t}\left[A_{i,j,k}^{n}\right]\right\} \left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right] \!\!\!\!\!\!\! \\ &&{\kern66pt}+\frac{b^{*}}{a^{*}}P_{\phi }\left\{ W_{t}\left[A_{i,j,K}^{n}\right]\right\}\\ &&{\kern66pt}\times\left. \left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\vphantom{\sum\limits_{k=2}^{K-1}}\right\} \\ &&{\kern6pt}+\,2\Delta z\Delta tr_{i}\left\{\vphantom{\sum\limits_{k=2}^{K-1}}\frac{b}{a}f_{i,j,1}^{n+\frac{1}{2} }\left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right]\right. \\ &&{ \kern58pt} +\,\sum\limits_{k=2}^{K-1}f_{i,j,k}^{n+\frac{1 }{2}}\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right] \\ &&{\kern58pt}+\left.\frac{b^{*}}{a^{*}}f_{i,j,K}^{n+\frac{1}{2} }\left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\vphantom{\sum\limits_{k=2}^{K-1}}\right\}. \label{eqA14} \end{array}$$

(31)

Summing Eq. (3131) over i and j, where 1 ≤ i ≤ I − 1 and 0 ≤ j ≤ J − 1, and multiplying the result by ΔrΔΦ, we obtain

$$ \begin{array}{rll} &&{\kern-6pt} 2\Delta r\Delta \phi \Delta z\sum\limits_{i=1}^{I-1}\sum\limits_{j=0}^{J-1}r_{i}\left\{\vphantom{\sum\limits_{k=2}^{K-1}}\frac{b}{a} \left(\left[A_{i,j,1}^{n+1}\right]^{2}-\left[A_{i,j,1}^{n}\right]^{2}\right)\right.\\ &&{\kern82pt} +\,\sum\limits_{k=2}^{K-1}\left(\left[A_{i,j,k}^{n+1}\right]^{2}-\left[A_{i,j,k}^{n}\right]^{2}\right) \\ &&{\kern82pt} +\left.\frac{b^{*}}{a^{*}}\left(\left[A_{i,j,K}^{n+1}\right]^{2}-\left[A_{i,j,K}^{n}\right]^{2}\right)\vphantom{\sum\limits_{k=2}^{K-1}}\right\} \\ &&\leq 2\Delta r\Delta \phi \Delta z\Delta tD_{r}\\ &&{\kern6pt} \times\,\sum\limits_{i=1}^{I-1}\sum\limits_{j=0}^{J-1}\left\{\frac{b}{a} P_{r}\left\{ W_{t}\left[A_{i,j,1}^{n}\right]\right\} \left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right]\vphantom{\sum\limits_{k=2}^{K-1}}\right. \\ &&{\kern54pt} +\sum\limits_{k=2}^{K-1}P_{r}\left\{ W_{t}\left[A_{i,j,k}^{n}\right]\right\}\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right]\\ &&{\kern54pt} +\left.\frac{b^{*}}{a^{*}}P_{r}\left\{ W_{t}\left[A_{i,j,K}^{n}\right]\right\} \left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\vphantom{\sum\limits_{k=2}^{K-1}}\right\} \\ &&{\kern6pt}+2\Delta r\Delta \phi \Delta z\Delta t\\ &&{\kern6pt} \times\sum\limits_{i=1}^{I-1}\sum\limits_{j=0}^{J-1}\frac{D_{\phi }}{r_{i}}\left\{\vphantom{\sum\limits_{k=2}^{K-1}}\frac{b}{a}P_{\phi }\left\{ W_{t}\left[A_{i,j,1}^{n}\right]\right\} \left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right] \right. \\ &&{\kern65pt}+\,\sum\limits_{k=2}^{K-1}P_{\phi }\left\{ W_{t}\left[A_{i,j,k}^{n}\right]\right\}\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right]\!\!\!\!\!\!\!\!\!\!\!\!\\ &&{\kern65pt} +\,\frac{b^{*}}{a^{*}}P_{\phi }\left\{W_{t}\left[A_{i,j,K}^{n}\right]\right\} \\ &&{\kern65pt}\times\left. \left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\vphantom{\sum\limits_{k=2}^{K-1}}\right\} \\ &&{\kern6pt}+\,2\Delta r\Delta \phi \Delta z\Delta t\\ &&{\kern6pt} \times\,\sum\limits_{i=1}^{I-1}\sum\limits_{j=0}^{J-1}r_{i}\left\{\vphantom{\sum\limits_{k=2}^{K-1}}\frac{b}{a}f_{i,j,1}^{n+\frac{1}{2}}\left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right]\right.\\ &&{\kern60pt} +\,\sum\limits_{k=2}^{K-1}f_{i,j,k}^{n+\frac{1}{2}}\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right] \\ &&{\kern60pt}+\left.\frac{b^{*}}{a^{*}}f_{i,j,K}^{n+\frac{1}{2} }\left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\vphantom{\sum\limits_{k=2}^{K-1}}\right\}. \label{eqA15} \end{array}$$

(32)

Using a similar argument in Eq. (30) together with Eq. (28), we have

$$ \begin{array}{rrl} &&{\kern-6pt} 2\Delta r\sum\limits_{j=0}^{J-1}\sum\limits_{i=1}^{I-1}P_{r}\left\{ W_{t}\left[A_{i,j,k}^{n}\right]\right\} \left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right]\\ &&=-\Delta r\sum\limits_{j=0}^{J-1}\sum\limits_{i=1}^{I}r_{i-\frac{1}{2}}\nabla _{\bar{r }}\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right]^{2}, \label{eqA16a} \end{array}$$

(33a)

$$ \begin{array}{rll} &&{\kern-6pt} 2\Delta \phi \sum\limits_{j=0}^{J-1}P_{\phi }\left\{ W_{t}\left[A_{i,j,k}^{n}\right]\right\} \left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right]\\ &&=-\Delta \phi \sum\limits_{j=0}^{J-1}\nabla _{\bar{\phi} }\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right]^{2}, \label{eqA16b} \end{array}$$

(33b)

for any 1 ≤ i ≤ I − 1 and 1 ≤ k ≤ K. Substituting Eqs. (33a) and (33b) into Eq. (32), we may drop the first two terms on the RHS of Eq. (32) and simplify it to

$$\begin{array}{rrl} &&{\kern-6pt} 2\Delta r\Delta \phi \Delta z\sum\limits_{i=1}^{I-1}\sum\limits_{j=0}^{J-1}r_{i}\left\{\vphantom{\sum\limits_{k=2}^{K-1}}\frac{b}{a} \left(\left[A_{i,j,1}^{n+1}\right]^{2}-\left[A_{i,j,1}^{n}\right]^{2}\right)\right.\\ &&{\kern80pt} +\,\sum\limits_{k=2}^{K-1}\left(\left[A_{i,j,k}^{n+1}\right]^{2}-\left[A_{i,j,k}^{n}\right]^{2}\right) \\ &&{\kern80pt} +\left.\frac{b^{*}}{a^{*}}\left(\left[A_{i,j,K}^{n+1}\right]^{2}-\left[A_{i,j,K}^{n}\right]^{2}\right)\vphantom{\sum\limits_{k=2}^{K-1}}\right\} \\ &&\leq 2\Delta r\Delta \phi \Delta z\Delta t\sum\limits_{i=1}^{I-1}\sum\limits_{j=0}^{J-1}\left\{\vphantom{\sum \limits_{k=2}^{K-1}}\frac{b}{a}f_{i,j,1}^{n+ \frac{1}{2}}\left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right]\right.\\ &&{\kern100pt} +\sum \limits_{k=2}^{K-1}f_{i,j,k}^{n+\frac{1}{2}}\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right] \\ &&{\kern100pt} +\left.\frac{b^{*}}{a^{*}}f_{i,j,K}^{n+\frac{1}{2} }\left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\vphantom{\sum \limits_{k=2}^{K-1}}\right\}. \\ \label{eqA17} \end{array}$$

(34)

Using Cauchy–Schwartz’s inequality ($2xy\leq \varepsilon x^{2}+\frac{1}{ \varepsilon }y^{2})$ and (x + y)² ≤ 2x ² + 2y ², we have

$$ \begin{array}{rll} &&{\kern-6pt} 2f_{i,j,1}^{n+\frac{1}{2}}\left[A_{i,j,1}^{n+1}+A_{i,j,1}^{n}\right]\\ &&\leq \frac{a}{b}\left[ f_{i,j,1}^{n+\frac{1}{2}}\right]^{2} +2\frac{b}{a}\left\{ \left[A_{i,j,1}^{n+1}\right]^{2}+\left[A_{i,j,1}^{n}\right]^{2}\right\} , \label{eqA18a} \end{array}$$

(35a)

$$ \begin{array}{rll} &&{\kern-6pt}2f_{i,j,K}^{n+\frac{1}{2}}\left[A_{i,j,K}^{n+1}+A_{i,j,K}^{n}\right]\\ &&\leq \frac{a^{*}}{ b^{*}}\left[f_{i,j,K}^{n+\frac{1}{2}}\right]^{2}+2\frac{b^{*}}{a^{*}}\left\{ \left[A_{i,j,K}^{n+1}\right]^{2}+\left[A_{i,j,K}^{n}\right]^{2}\right\} , \label{eqA18b} \end{array}$$

(35b)

$$ \begin{array}{rll} &&{\kern-6pt}2f_{i,j,k}^{n+\frac{1}{2}}\left[A_{i,j,k}^{n+1}+A_{i,j,k}^{n}\right]\\ [4pt] &&\leq \left[f_{i,j,k}^{n+ \frac{1}{2}}\right]^{2}+2\left\{ \left[A_{i,j,k}^{n+1}\right]^{2}+\left[A_{i,j,k}^{n}\right]^{2}\right\} , \label{eqA18c} \end{array}$$

(35c)

for any 2 ≤ k ≤ K − 1. Substituting Eqs. (35a)–(35c) into Eq. (34) and then denoting

$$\begin{array}{rll} A\left( n\right) &=&2\Delta r\Delta \phi \Delta z\sum\limits_{i=1}^{I-1}\sum\limits_{j=0}^{J-1}r_{i}\left\{\frac{b}{a}\left[A_{i,j,1}^{n+1}\right]^{2}\!+\!\sum\limits_{k=2}^{K-1}\left[A_{i,j,k}^{n+1}\right]^{2}\right.\!\!\!\!\!\!\!\\ &&{\kern88pt} +\left.\frac{b^{*} }{a^{*}}\left[A_{i,j,K}^{n+1}\right]^{2}\vphantom{\sum\limits_{k=2}^{K-1}}\right\} \label{eqA19a} \end{array}$$

(36a)

and

$$\begin{array}{rll} F\left( n\right) &=&\Delta r\Delta \phi \Delta z\sum\limits_{i=1}^{I-1}\sum\limits_{j=0}^{J-1}r_{i}\left\{\frac{b}{a}\left[ f_{i,j,1}^{n+\frac{1}{2}}\right]^{2}+\sum\limits_{k=2}^{K-1}\left[f_{i,j,k}^{n+\frac{1}{2}}\right]^{2}\right.\\ &&{\kern84pt} +\left.\frac{b^{*}}{a^{*}}\left[f_{i,j,K}^{n+\frac{1}{2}}\right]^{2}\vphantom{\sum\limits_{k=2}^{K-1}}\right\}, \label{eqA19b} \end{array}$$

(36b)

we can further simplify Eq. (34) to

$$ \begin{array}{rll} &&{\kern-6pt} A(n+1)\\ &&\leq \frac{1+\Delta t}{1-\Delta t}A\left( n\right) +\frac{\Delta t}{1-\Delta t}F\left( n\right) \\ &&\leq \frac{1+\Delta t}{1-\Delta t}\left[\frac{1+\Delta t}{1-\Delta t}A\left( n-1\right) +\frac{\Delta t}{1-\Delta t}F\left(n-1\right) \right]\\ &&{\kern6pt} +\,\frac{\Delta t}{1-\Delta t}F\left( n\right) \\ &&\leq \left(\frac{1+\Delta t}{1-\Delta t}\right)^{n+1}A\left( 0\right)+\frac{ \Delta t}{1-\Delta t} \\ &&{\kern6pt} \times\left[1+\frac{1+\Delta t}{1-\Delta t}+\cdots +\left(\frac{ 1+\Delta t}{1-\Delta t}\right)^{n}\right]\max_{0\leq n_{1}\leq n}F\left( n_{1}\right) \\ &&\leq \left(\frac{1+\Delta t}{1-\Delta t}\right)^{n+1}\left[A\left( 0\right) +\max_{0\leq n_{1}\leq n}F\left( n_{1}\right) \right]. \label{eqA20} \end{array}$$

(37)

Using inequalities $\left( 1+\varepsilon \right) ^{n}\leq e^{n\varepsilon },\varepsilon >0$, and ( 1 − ε) ^− 1 ≤ $e^{2\varepsilon },0<\varepsilon <\frac{1}{2},$ we obtain

$$ \begin{array}{rll} A\left( n+1\right) &\leq & e^{3n\Delta t}\left[A\left( 0\right) +\max_{0\leq n_{1}\leq n}F\left( n_{1}\right) \right]\\ &\leq& e^{3t_{0}}\left[A\left(0\right) +\max_{0\leq n_{1}\leq n}F\left( n_{1}\right) \right], \label{eqA21} \end{array}$$

(38)

for any $0\leq \left( n+1\right) \Delta t\leq t_{0}$, implying that the scheme is unconditionally stable with respect to the initial condition and source term.□