- Email: [email protected]
Boundary Value Problems of Advection–Diffusion–Reaction Type
CHAPTER
FINITE ELEMENT APPROXIMATIONS OF INITIAL VALUE/BOUNDARY VALUE PROBLEMS OF ADVECTION– DIFFUSION–REACTION TYPE
24
CHAPTER POINTS • In this chapter we consider the initial value/boundary value problem of advection–diffusion–reaction type representing the mathematical model of a balance law with a linear net production term.
• The numerical discretization of the initial value/boundary value problem is conducted using the θ-method for approximating the time derivative and the mixed GFEM for the spatial approximation.
• The stability and convergence properties of the numerical formulation are analyzed and computational examples are provided to illustrate the theoretical conclusions.
• The computational procedure is applied to the simulation of the propagation of an action potential along an axon and of ion electrodynamics in a transmembrane ion channel.
In Chapter 23 we have illustrated the Galerkin finite element formulation for the discretization of a boundary value problem with advective, diffusive, and reactive terms. In this chapter we consider the initial value/boundary value problem of advection–diffusion–reaction type that represents the mathematical model of a balance law with a linear net production term. We first approximate the time derivative of the initial value/boundary value problem using the θ -method introduced in Section 3.5.2.2, whereas the spatial approximation of the problem is conducted using the mixed GFEM introduced in Chapter 23. Then, we analyze the properties of stability and convergence of the numerical formulation and we provide a series of computational examples to support the theoretical conclusions. We conclude the chapter by applying the computational procedure to the simulation of two biophysical problems. The first problem is the propagation of an action potential along an axon using the cable equation model introduced in Chapter 18. The second problem is the electrodynamics of ions in a transmembrane ion channel using the Poisson–Nernst–Planck model introduced in Section 13.6.3. A Comprehensive Physically Based Approach to Modeling in Bioengineering and Life Sciences. https://doi.org/10.1016/B978-0-12-812518-2.00037-8 Copyright © 2019 Elsevier Inc. All rights reserved.
579
580
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
24.1 THE TIME-DEPENDENT ADVECTION–DIFFUSION–REACTION PROBLEM In this section we address the study of the one-dimensional advection–diffusion–reaction problem (22.1) with respect to the temporal coordinate. To simplify the presentation, we set γ∂ = 0, β∂ = 0, and α∂ = 1, corresponding to homogeneous Dirichlet boundary conditions at the endpoints of . Therefore, our considered initial value/boundary value problem reads as follows: Find u(x, t) for all (x, t) ∈ QT such that ∂u ∂ J + + κu = g ∂t ∂x ∂u J = uv − μ ∂x u(x, 0) = u0 (x) u(x, t) = 0
∀(x, t) ∈ QT ,
(24.1a)
∀(x, t) ∈ QT ,
(24.1b)
∀x ∈ , ∀(x, t) ∈ ∂ × (0, T ).
(24.1c) (24.1d)
As an example, the initial value/boundary value problem (24.1) can be regarded as the heat equation in the form of a balance law. In this physical picture, the dependent variable u has the meaning of the temperature of the medium in which the heat conduction process takes place, whereas J is the heat flux density. The initial value/boundary value problem (24.1) represents the mathematical model of the spatial and temporal diffusion and passive transport of heat in a medium of length L with thermal conductivity μ and heat transport velocity v, subject to a heat production g, a heat consumption rate κ, an initial thermal distribution u0 , and a boundary temperature equal to zero for all t ∈ (0, T ). Property 24.1 (Asymptotic stability). Set κ = 0 and g = 0 in (24.1a). Then we have lim u(x, t) = 0
t→∞
∀x ∈ .
(24.1e)
Relation (24.1e) expresses the fact that the initial value/boundary value problem (24.1) is asymptotically stable. Asymptotic stability of the initial value/boundary value problem is a consequence of the fact that μ is strictly positive and expresses the physical property that the system modeled by (24.1) is dissipative. Let us now use the θ -method introduced in Section 22.2 for the discretization of the initial value/boundary value problem (24.1) with respect to spatial and temporal variables. In the following, we denote by C a positive constant independent of h and t whose value may not be the same at each occurrence. Using the notation of Section 22.2 and Section 23.1.3, the linear algebraic system that has to be solved at each time level tk , k = 0, . . . , NT − 1, reads Aθ uk+1 = fk+1 θ ,
(24.1f)
1 M, t
(24.1g)
where Aθ = θ (Kdmix + Kamix + Krmix ) +
fk+1 = −(1 − θ )(Kdmix + Kamix + Krmix )uk + θ
1 Muk + θ gk+1 + (1 − θ )gk , t
(24.1h)
24.1 THE TIME-DEPENDENT ADVECTION–DIFFUSION–REACTION MODEL
581
having introduced the mass matrix M = ij
L
i, j = 1, . . . , Nh .
ϕj (x)ϕi (x) dx,
(24.1i)
0
The Matlab function 29.50 computes the mass matrix (24.1i) using the Cavalieri–Simpson quadrature formula (23.32c). We note that the mass matrix (24.1i) is a special instance of the reaction matrix (23.30c) corresponding to setting σ (x) = 1 for all x ∈ . Therefore, the lumping procedure illustrated in Section 23.4.1 for the diagonalization of the reaction matrix can be applied also to the mass matrix. This leads to = hδij , i, j = 1, . . . , Nh , (24.1j) M ij
where δij is the Kronecker symbol. The Matlab function 29.51 computes the lumped mass matrix (24.1j). In the above relations, θ is a parameter to be chosen in the interval [0, 1]. In particular, we note that: 1. If θ = 0, relation (24.1f) becomes an explicit formula to compute the solution at time level tk+1 using information only at the level tk . The corresponding scheme is called forward Euler or explicit Euler method. 2. If θ > 0, the family of schemes associated with (24.1f) takes the name of implicit methods because to compute uk+1 we need to solve a linear system using information at both levels tk and tk+1 . Definition 24.1 (Stability). Set κ = 0 and g = 0 in (24.1a). The θ -method (24.1f) is said to be stable (or also, asymptotically stable) if the sequence uk is such that lim uk = 0.
k→∞
(24.1k)
Condition (24.1k) is the discrete counterpart of (24.1e) and is not automatically inherited by the θ-method. In this sense, the following condition needs to be satisfied (for the proof, see [259, Chapter 11] and [258, Chapter 13]). Theorem 24.1 (Stability of the θ -method). If θ ≥ 1/2 the θ -method verifies (24.1k) for every value of t. In this case we say that the θ-method is unconditionally stable. If 0 ≤ θ < 1/2 the θ-method verifies (24.1k) only if t ≤
Ch2 . μ0 (1 − 2θ )
(24.1l)
In this case we say that the θ-method is conditionally stable. Conditional stability can be a severe computational constraint for the θ-method because for a small value of h (which is typically needed to achieve a reasonable spatial accuracy of the approximation) the corresponding value of t to make the θ -method stable turns out to be even smaller. Therefore, the use of an implicit method is in order thanks to its unconditional stability. The choice of such a method is governed by the next important result.
582
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
Theorem 24.2 (Convergence of the θ-method). Under suitable regularity assumptions on the data in (22.1), we have uk − ukh L2 () ≤ u0 − u0h L2 () + C(u0 , f, u) h2 + t p(θ) (24.1m) ∀k ∈ [1, NT ], where p(1/2) = 2 and p(θ) = 1 for all θ = 1/2. Theorem 24.2 tells us that: 1. if lim u0 − u0h L2 () = 0, then the θ -method is convergent; h→0
2. the θ -method is convergent with order 2 with respect to both spatial and time coordinates only in the case of the Crank–Nicolson scheme. All the other implicit schemes with θ > 1/2 converge linearly with respect to t. Quadratic convergence with respect to the spatial coordinate is a consequence of (23.91).
24.1.1 CONVERGENCE EXAMPLES IN THE TIME-DEPENDENT CASE In the following example we study the initial value/boundary value problem (24.1) using the θ-method in the case where μ = 1, v = 0, κ = 0, = (0, 1), IT = [0, 1], u0 (x) = sin(2πx), and g(x, t) = (4π 2 − λ)e−λt sin(2πx),
x ∈ , and t ∈ IT ,
with λ > 0, in such a way that the exact solution is u(x, t) = e−λt sin(2πx).
FIGURE 24.1 The time-dependent problem, convergence of the θ-method as a function of h with t = 1/100. The spatial mesh size is h = (10 · k)−1 , k = 1, . . . , 20; λ is equal to 1.
Figs. 24.1 and 24.2 illustrate the results obtained by solving the initial value/boundary value problem with the backward Euler and Crank–Nicolson methods. The parameter λ is set equal to 1. T Comparison between the two schemes is carried out by computing the quantity uNT − uN h L2 () . Fig. 24.1 is obtained by running the Matlab script 28.47 and performs a study of convergence with
24.2 SIMULATION OF THE PROPAGATION OF AN ACTION POTENTIAL
583
FIGURE 24.2 The time-dependent problem, convergence of the θ-method as a function of t with h = 1/500, the temporal mesh size is t = (10 · k)−1 , k = 1, . . . , 20.
respect to h having fixed t = 1/100, whereas the spatial grid size is allowed to vary in the interval [1/10, 1/200] by setting h = (10 · k)−1 , k = 1, . . . , 20. Conversely, Fig. 24.2 is obtained by running the Matlab script 28.48 and performs a study of convergence with respect to t having fixed h = 1/500, whereas the time grid size is allowed to vary in the interval [1/10, 1/200] by setting t = (10 · k)−1 , k = 1, . . . , 20. The two convergence studies conducted above clearly indicate that the Crank–Nicolson method is far more accurate than the backward Euler method with respect to both h and t. In particular, we note that, no matter how small h, the error of the backward Euler method in Fig. 24.1 is dominated by the contribution O( t). Since this latter quantity is fixed (equal to 0.01 in T the present case) the error norm uNT − uN h L2 () does not tend to zero as h tends to zero.
24.2 SIMULATION OF THE PROPAGATION OF AN ACTION POTENTIAL ALONG AN AXON In this section we numerically investigate the time-dependent GFEM illustrated in Section 24.1 applied to the simulation of the propagation of an action potential along an ideal axon of length L = 10 cm and radius a = 2.5 mm. The computational procedure consists of the use of (i) the functional iteration techniques illustrated in Section 21.5.1; (ii) the time discretization methods illustrated in the present chapter, and (iii) the mixed GFEM illustrated in Chapter 23.
24.2.1 MODEL DATA The axon membrane thickness is tm = 10 nm, the load output resistance Rout is set equal to +∞, so that the terminal of the axon at x = L is, from the electrical standpoint, an open circuit and no electric current is allowed to flow out of it. (See Table 24.1.)
584
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
Table 24.1 Values of model parameters in the simulation of the propagation of an action potential along an axon using the cable equation model. The temperature of the system under investigation is = 298.15 K. Transmembrane conduction current is assumed to be driven by potassium, sodium, and chloride ions. Model parameter
Description
Value
r εm t0 , t1 , t2 , t3 Tf in Rout a L tm cm V1 , V2
Dielectric relative constant Time characteristic of the axon potential Final simulation time Load output resistance Radius of the ideal axon Length of the ideal axon Membrane thickness of the ideal axon Unit area capacitance of the ideal axon Electric characteristic of the axon potential
80 [·] 0,5,7.5,10 ms 15 ms +∞ 2.5 · 10−3 m 0.1 m 10−8 m 10−2 F m−2 28.73, −91.27 mV
Table 24.2 Values of intracellular and extracellular ion concentrations in the study of propagation of an action potential along an axon in the human body (cf. [107, Table 1.1].) Model parameter Description K+
(in) CK (in) CNa (in) CCl CK(out) (out) CNa (out) CCl
concentration in the intracellular region
Value 150 mM
Na+ concentration in the intracellular region
50 mM
Cl−
10 mM
concentration in the intracellular region
K+ concentration in the extracellular region Na+
5 mM
concentration in the extracellular region 460 mM
Cl− concentration in the extracellular region
125 mM
The resting potential of the cell can be computed by using the Goldman equation (17.24b) in which the ion permeabilities PK , PNa , and PCl satisfy relation (17.39b). Setting T = 300 K and using the values of the intracellular and extracellular concentrations listed in Table 24.2, we obtain ψeq = −71.27 mV.
(24.2)
24.2.2 ACTION POTENTIAL The input voltage stimulus V in is the piecewise step function illustrated in Fig. 24.3. This function is a schematic representation of the action potential which develops in an excitable cell such as a neuron or a cardiac myocyte. The applied voltage is equal to the resting potential of the cell ψeq from time t = t0 = 0 ms until time t1 = 5 ms. In these conditions, the intracellular potential is lower than the extracellular potential so that the “−” pole is the intracellular site and the “+” pole is the extracellular site. At time t = t1 , the applied voltage is suddenly raised to the value V1 = ψeq + 100 mV = +28.73 mV and is kept equal to V1 until time t2 = 7.5 ms. In these conditions, a change of polarity takes place because the intracellular potential is higher than the extracellular potential so that the “+” pole is the intracellular
24.2 SIMULATION OF THE PROPAGATION OF AN ACTION POTENTIAL
585
FIGURE 24.3 Input voltage stimulus. The piecewise constant function shown in the figure schematically represents an action potential. The jump between the resting potential to the positive value V1 = 28.73 mV corresponds to the depolarization phase. The jump from V1 to the value V2 = −91.27 mV corresponds to the hyperpolarization phase.
site and the “−” pole is the extracellular site. This change of polarity is the result of a flow of a positive ion charge, such as Na+ , that was attracted to the cytoplasm by the “+” pole in the resting condition. This phase of the action potential is referred to as depolarization (see [107, Chapter 1]). At time t = t2 , the voltage stimulus is again suddenly changed to the value V2 = ψeq − 20 mV = −91.27 mV and is kept equal to V2 until time t3 = 10 ms. In these conditions, another change of polarity takes place because the intracellular potential is more negative than the resting potential so that the “−” pole is the intracellular site and the “+” pole is the extracellular site. This change of polarity is the result of a flow of a negative ion charge, such as Cl− , that was attracted to the cytoplasm by the “+” pole in the depolarization condition. This phase of the action potential is referred to as hyperpolarization (see [107, Chapter 1]). At time t = t3 , the voltage stimulus is again suddenly changed to the value of the resting potential ψeq and is kept equal to ψeq until the final time t = Tf in = 15 ms.
24.2.3 OVERVIEW In Chapter 17 we illustrated two different models to characterize the transmembrane current, namely, the linear resistor formulation (17.44) and the nonlinear Goldman–Hodgkin–Katz (GHK) formulation (17.51). In the following paragraphs, we utilize the functional iteration (21.16) and its time and spatial discretization to compare the solutions obtained using the two models. The numerical implementation of the computational algorithm is performed by the Matlab script 28.49. The Matlab function 29.28 implements the solution of the problem using the linear resistor formulation, whereas the Matlab function 29.29 implements the solution of the problem using the GHK formulation. In both functions, the time discretization of (18.2) is performed using the backward Euler method illustrated in Section 3.5.2.1, whereas the spatial discretization is performed using the mixed GFEM illustrated in Chapter 23. A uniform time step t = (Tend − t0 )/250 = 6 · 10−5 s and a uniform spatial discretization parameter h = L/100 = 0.1 cm are used in the numerical computations. The Matlab function 29.27 computes the voltage stimulus graphically represented in Fig. 24.3.
586
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
Simulation With the Linear Resistor Model In this paragraph we perform the numerical simulation of the propagation of the action potential elicited by the voltage stimulus of Fig. 24.3, by representing the axon transmembrane current density through the linear resistor model (17.44). We consider three values of the longitudinal resistivity, and we set SW (case (a)), ρ = 10+3 ρ SW (case (b)), and ρ = 10−3 ρ SW (case (c)). To use the linear ρin = ρin in in in in resistor model, we need to determine the ion specific conductance gα (unit: S m−2 ), α ∈ {K, Na, Cl}. With this aim, we use the definition of local ion conductivity (13.75d) in the following average form: σα = q|zα |μel α
nα (0) + nα (tm ) nin + nout α = q|zα |2 Dα α , 2tm 2tm
α ∈ {K, Na, Cl} ,
(24.3)
out where nin α and nα are the intra and extracellular number densities of the αth ionic species and Einstein’s relation (13.69) is used to express the ion electric mobility μel α in favor of the ion diffusion coefficient Dα . Results are shown in Figs. 24.4, 24.5, and 24.6.
FIGURE 24.4 Simulation of axon potential evolution using the cable equation model (18.2) supplied by the linear resistor model SW = 2 · 10−1 m, for the transmembrane current. In this simulation the axoplasm resistivity ρin is set equal to ρin −1 corresponding to an axoplasm conductivity σin = 5 S m . Left panel: view from x = L. Right panel: view from x = 0.
Case (a). Fig. 24.4 shows the axon membrane potential ψm (x, t) for x ∈ [0, L] and for t ∈ [0, Tend ], computed by solving (21.16) with the linear resistor model for the axon transmembrane current density. SW = 2 · 10−1 m, the resistivity of salt The value of the axoplasm resistivity ρin is set equal to ρin ◦ water at 20 C. This value is very close to the data shown in [83], where the experimentally measured SW . Results indicate the resistivity of the extruded squid axon turned out to be equal to (1.4 ± 0.2)ρin following: • At t = 0, the membrane potential undergoes a fast increase in the neighborhood of x = 0 and reaches a constant value close to 0 mV along the axon length because of a flux of positive charge inside the axon.
24.2 SIMULATION OF THE PROPAGATION OF AN ACTION POTENTIAL
587
FIGURE 24.5 Simulation of axon potential evolution using the cable equation model (18.2) supplied by the linear resistor model SW = 2 · 102 m, for the transmembrane current. In this simulation the axoplasm resistivity ρin is set equal to 103 ρin corresponding to an axoplasm conductivity σin = 5 · 10−3 S m−1 . Left panel: view from x = L. Right panel: view from x = 0.
FIGURE 24.6 Simulation of axon potential evolution using the cable equation model (18.2) supplied by the linear resistor SW = model for the transmembrane current. In this simulation the axoplasm resistivity ρin is set equal to 10−3 ρin −4 3 −1 2 · 10 m, corresponding to an axoplasm conductivity σin = 5 · 10 S m . For each t ∈ [0, Tf in ], the spatial distribution of the membrane potential along the axon is constant because of the elevated longitudinal conductivity.
• As t increases, the spatial distribution of the membrane potential remains uniform but with a smaller value compared to t = 0 because of a flux of part of the accumulated positive charge out of the intracellular axon region. • At t = t1 , the membrane potential rises abruptly at x = 0 and then monotonically decreases along the axon length, passing from positive values to negative values, because of a positive accumulated charge inside the axon near x = 0 and a negative accumulated charge for larger values of x. The membrane potential exhibits an opposite behavior at t = t2 .
588
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
• At t = t3 , the membrane potential returns abruptly to its resting value, which is less negative than V2 . This induces an accumulation of positive charge inside the intracellular region of the axon and the membrane potential becomes uniform along the axon length, reaching a value of about −55 mV. • For larger values of t, the accumulated positive charge flows out of the axon, thereby giving rise to more negative values of the membrane potential. After a certain time instant, all the accumulated positive charge has flown out of the axon, so that the membrane potential becomes flat along the whole axon length and equal to the boundary condition at x = 0. • For all t ∈ [0, Tf in ], the spatial distribution of the membrane potential becomes flat at the end of the axon because of the no-flux boundary condition at x = L. Case (b). Fig. 24.5 shows the axon membrane potential computed by solving (21.16) with the linear resistor model and with an axoplasm resistivity increased by a factor of 103 with respect to the previous case. Results indicate that the membrane potential undergoes a practically linear decrease with respect to time for each x ∈ [0, L]. This behavior is due to the elevated value of the longitudinal axon resistivity which considerably slows down the transmembrane flow of positive and negative charge for each x ∈ [0, L]. We also note that the axon is far from resting conditions at t = Tf in because the membrane potential is about −50 mV. Resting conditions would be reached only for a significantly larger value of the observational time Tf in . Case (c). Fig. 24.6 shows the axon membrane potential ψm (x, t) for x ∈ [0, L] and for t ∈ [0, Tend ] in SW = the case where the axoplasm resistivity is reduced by a factor of 103 , so that we have ρin = 10−8 ρin −9 2 · 10 m. Results indicate that for each value of t ∈ [0, Tend ] the spatial distribution of the membrane potential along the axon is constant because the elevated longitudinal conductivity gives rise to an almost instantaneous transport of the value of the Dirichlet boundary condition ψm (0, t) = V in (t).
Simulation With the GHK Model In this paragraph we perform the numerical simulation of the propagation of the action potential elicited by the voltage stimulus of Fig. 24.3, by representing the axon transmembrane current density through the GHK model (17.51). We consider the same conditions as in Paragraph ‘Simulation With the Linear SW (case (a)), ρ = 10+3 ρ SW (case (b)), and ρ = 10−3 ρ SW Resistor Model’, that is, we set ρin = ρin in in in in (case (c)). Computations show that four iterations were needed to satisfy the convergence check (21.17) with a tolerance toll = 10−5 . Results are shown in Figs. 24.7, 24.8, and 24.9. Model predictions in case (c) cannot be distinguished from the predictions given by the linear resistor model in the same conditions. Model predictions in cases (a) and (b) are qualitatively similar to the corresponding simulations with the linear resistor model. In case (a), we see that: • the increase of membrane potential at t = 0 in the neighborhood of x = 0 is less sharp than predicted by the linear resistor model; • the membrane potential at x = L is less negative than predicted by the linear resistor model; • the membrane potential distribution at t = Tf in is not as flat as predicted by the linear resistor model. This indicates that, according to the GHK model, the axon has not reached resting conditions as in the case of the linear resistor model. In particular, the value of the membrane potential predicted by the GHK model at x = L is less negative than that predicted by the linear resistor model. Similar considerations apply for case (b). In particular, we see that:
24.2 SIMULATION OF THE PROPAGATION OF AN ACTION POTENTIAL
589
FIGURE 24.7 Simulation of axon potential evolution using the cable equation model (18.2) supplied by the GHK model for the SW = 2 · 10−1 m, corretransmembrane current. In this simulation the axoplasm resistivity ρin is set equal to ρin sponding to an axoplasm conductivity σin = 5 S m−1 . Left panel: view from x = L. Right panel: view from x = 0.
FIGURE 24.8 Simulation of axon potential evolution using the cable equation model (18.2) supplied by the GHK model for the SW = 2 · 102 m, corretransmembrane current. In this simulation the axoplasm resistivity ρin is set equal to 103 ρin −3 −1 sponding to an axoplasm conductivity σin = 5 · 10 S m . Left panel: view from x = L. Right panel: view from x = 0.
• for each x ∈ [0, L] the decay in time of the membrane potential is significantly smaller than predicted by the linear resistor model; • the resulting membrane potential at x = L predicted by the GHK model is much less negative than in the case of the linear resistor model, so that a much larger observational time is needed to reach resting conditions in the axon.
590
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
FIGURE 24.9 Simulation of axon potential evolution using the cable equation model (18.2) supplied by the GHK model for the SW = 2 · 10−4 m, transmembrane current. In this simulation the axoplasm resistivity ρin is set equal to 10−3 ρin 3 −1 corresponding to an axoplasm conductivity σin = 5 · 10 S m . For each t ∈ [0, Tf in ], the spatial distribution of the membrane potential along the axon is constant because of the elevated longitudinal conductivity.
24.2.4 COMPARISON AND CONCLUDING REMARKS The same considerations developed in Section 17.7.4.3 can be applied to discuss the quality and reliability of model predictions obtained in Paragraphs ‘Simulation With the Linear Resistor Model’ and ‘Simulation With the GHK Model’ using the linear resistor model and the GHK model, respectively. In particular, the accuracy of the computed action potential profiles may be warranted through the following two approaches: • an error analysis against an exact solution, and/or • comparison with experimental data, whenever available. As far as the first approach, we refer to [223], whereas for the second approach we refer, e.g., to [317].
24.3 POISSON–NERNST–PLANCK SIMULATION OF AN IDEAL ONE-DIMENSIONAL ION CHANNEL In this section we apply the Poisson–Nernst–Planck model illustrated in Section 13.6.3 to the simulation of an ideal one-dimensional ion channel of length L = 10 · 10−7 cm in which only the Ca++ and Cl− ions are considered. We warn the reader that the scope of this section is not the mathematical modeling of a realistic ion channel as addressed in Chapter 26, but, rather, the numerical investigation of the time-dependent GFEM (see Section 24.1) and the application of the functional iteration (21.22). To devise the computational procedure, we apply (i) the functional iteration techniques illustrated in Section 21.5.2, (ii) the time discretization methods illustrated in the present chapter, and (iii) the mixed GFEM illustrated in Chapter 23. All the computations reported in this section are conducted by running the Matlab code bio_1d, which can be downloaded at [210]. Fig. 24.10 shows the scheme of the one-dimensional ion channel used in the simulation of the Poisson–Nernst–Planck model.
24.3 POISSON–NERNST–PLANCK SIMULATION OF AN IDEAL
591
FIGURE 24.10 Scheme of the one-dimensional ion channel.
24.3.1 MODEL DATA In the ion channel represented in Fig. 24.10, Ca++ and Cl− ions flow under the action of diffusion and electric forces (see Section 16.2). The channel is represented by a medium having the same dielectric constant of water, so that in the Poisson equation (13.74c) we can set εm = ε0 εr,w , where ε0 is the dielectric permittivity of vacuum and εr,w = 80 is the relative dielectric permittivity of water. The presence of a surface charge density that characterizes the structure of any realistic ion channel is replaced by a volumetric negative charge density ρf ixed = −9.6 · 1018 cm−3 . Ion transport parameters such as diffusivity, electrical mobility, and electric charge are reported in Table 24.3. The simulation time spans the interval IT = [tinitial , tf inal ] with tinitial = 1 · 10−10 s and tf inal = 10 · 10−7 s and the temporal discretization is obtained by dividing IT l into a uniform partition of NT = 200 elements, so that the temporal step size is t = (tf inal − tinitial )/NT . The spatial partition Th used by the GFEM consists of a uniform subdivision of the channel length into Mh = 500 elements, so that h = L/Mh . Table 24.3 Parameters used in the simulation of an ideal one-dimensional ion chan-
nel in which Ca++ and Cl− ions are flowing. The temperature of the system under investigation is = 298.15 K.
Parameter Description
Value
r εm tinitial tf inal L ρf ixed zCa++ DCa++ μCa++ zCl− DCl− μCl− Mh NT
80 [·] 1 · 10−10 s 1 · 10−7 s 10 · 10−7 cm −9.6 · 1018 cm−3 +2 · 1.6019 · 10−19 C 1 · 10−6 cm2 s−1 7.735 · 10−5 cm2 V−1 s−1 −1.6019 · 10−19 C 0 1 · 10−6 cm2 s−1 3.868 · 10−5 cm2 V−1 s−1 500 [·] 200 [·]
Dielectric relative constant Initial simulation time Final simulation time Length of the ideal one-dimensional channel Fixed charge concentration in the channel Electric charge of the Ca++ ion Diffusivity of the Ca++ ion Electric mobility of the Ca++ ion Electric charge of the Cl− ion Diffusivity of the Cl− ion Electric mobility of the Cl− ion Number of elements in the finite element partioning Number of time intervals used in the θ-method
The dynamic behavior of the ion channel is investigated using the time-dependent Poisson–Nernst– Planck system (13.74). In the present case we have Mion = 2, so that the Poisson–Nernst–Planck model is constituted by five equations: the Poisson equation (13.74c), two Eqs. (13.74a) expressing the balance of mass for each ion, and two Eqs. (13.74b) expressing the balance of linear momentum for
592
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
each ion. Dirichlet boundary conditions are applied at both endpoints of the channel to Eqs. (13.74c) and (13.74a) and are assumed to be constant in time. Table 24.4 reports the data for the boundary and initial conditions used in the simulations. The channel is subjected to a bias condition similar to that considered in Chapter 26 by applying at x = L a negative voltage of −2.5 · 10−3 V and at x = 0 a voltage equal to 0 V, in such a way that the externally applied membrane potential is ψm = 0 − (−2.5 · 10−3 ) = 2.5 · 10−3 V. A concentration gradient is forced between the internal and external cell regions for both ions but in opposite directions. This gradient is obtained by using similar biophysical values for the ion concentrations at x = 0 and x = L as calculated in Section 26.2. Table 24.4 Initial and boundary conditions applied to the simulation of an ideal onedimensional ion channel in which Ca++ and Cl− ions are flowing.
Initial or boundary condition Description ψ0 ψL ψ(t0 ) Ca++ 0 Ca++ L Ca++ (t0 ) Cl− 0 Cl− L Cl− (t0 )
Electric potential at x = 0 and in [tinitial , tf inal ] Electric potential at x = L and in [tinitial , tf inal ] Initial value of the electric potential in = [0, L] Ca++ concentration at x = 0 and in [tinitial , tf inal ] Ca++ concentration at x = L and in [tinitial , tf inal ] Initial value of the Ca++ concentration in = [0, L] Cl− concentration at x = 0 and in [tinitial , tf inal ] Cl− concentration at x = L and in [tinitial , tf inal ] Initial value of the Cl− concentration in = [0, L]
Value 0V −2.5 · 10−3 V 0V 4 mM 1.33 mM 1.33 mM 107 mM 131 mM 107 mM
24.3.2 SIMULATION RESULTS In this section we report the simulation results of ion transport in the one-dimensional channel described above. Despite the geometric simplicity of the considered one-dimensional channel, the test reveals some interesting model implications and numerical observations. Fig. 24.11 reports the convergence behavior of the Gummel map within the functional iteration (21.22) obtained by solving the Poisson–Nernst–Planck model (13.74) up to a tolerance toll = 10−8 for each discretization time step (cf. the termination condition, Eq. (21.27)). Except for the first time step, in the remaining time intervals the number of Gummel iterations to reach convergence up to the given tolerance remains constant. The evolution of the electric quantities is reported in Fig. 24.12, where the black dashed lines refer to the profiles computed at the final simulation time t = tf inal . The left panel shows the electric potential evolution; the values at the domain endpoints are kept constant thanks to the Dirichlet boundary conditions while within the channel the electric potential is found to have always a nonlinear shape. As a consequence, the electric field is not constant as clearly visualized in the right panel of Fig. 24.12; in the first part of the channel, the electric field is directed to the right (extracellular region) while in the final part of the channel the electric field is directed to the left (intracellular region). We note that in the center of the channel the electric field strength is much lower than in the regions close to the boundary. This finding is in contrast with the assumption of a constant electric field made in Section 17.6.2 to derive the GHK model. Thus, we see that the present investigation shows how the assumptions made to derive the GHK model to calculate the transmembrane current may fail to be verified in all condi-
24.3 POISSON–NERNST–PLANCK SIMULATION OF AN IDEAL
593
FIGURE 24.11 Number of iterations required by the Gummel map to reach convergence in the functional iteration (21.22) up to a tolerance toll = 10−8 for each discretization time step.
FIGURE 24.12 Simulation of the evolution of electric quantities in an ideal one-dimensional ion channel of length L = 15 · 10−7 cm in which only Ca++ and Cl− ions are considered. Results have been obtained with the Poisson– Nernst–Planck model (13.74). The black dashed lines refer to the profiles at the final simulation time t = tf inal . (Initial conditions are omitted). Left panel: electric potential. Right panel: electric field.
tions. This is a general indication, and, as a consequence, the application of any mathematical model to describe a given experimental condition needs always to be carefully verified. Fig. 24.13 reports the evolution of ion concentrations. The black dashed lines refer to the profiles computed at the final simulation time t = tf inal . The evolution of the Ca++ concentration is shown in the left panel. Looking at the electric field close to the left boundary of the channel (x = 0), we see that the drift component of the Ca++ current density is directed towards the right boundary (x = L) because of the positive sign of the electric field. Vice versa, close to the right boundary (x = L) the electric field is negative and pushes the Ca++ ions away from it. This combination of effects results in an accumulation of Ca++ ions within the center of the first half of the channel.
594
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
The evolution of the Cl− concentration is shown in the right panel of Fig. 24.13. Due to the negative sign of the Cl− ion charge, the electric field pushes the chloride ions towards the boundaries (at x = 0 and x = L), giving rise to the formation of two pronounced boundary layers and reducing the chloride concentration in the center of the channel region.
FIGURE 24.13 Simulation of the evolution of ion concentrations in an ideal one-dimensional ion channel of length L = 15 · 10−7 cm in which only the Ca++ and Cl− ions are considered. Results have been obtained with the Poisson– Nernst–Planck model (13.74). The black dashed lines refer to the profiles at the final simulation time t = tf inal (initial conditions are omitted). Left panel: Ca++ concentration. Right panel: Cl− concentration.
Among the several indications we have obtained in the discussion of this section, we can conclude that the adoption of a more sophisticated model than the lumped formulations illustrated in Chapter 17 has the disadvantage that it requires advanced mathematical methods for its solution, but, at the same time, it may help the user to better investigate the underlying biophysical mechanisms and verify a posteriori the correctness of the assumptions made to derive the reduced-order model approach.
24.4 CONCLUDING REMARKS BIO-WARNING 24.1 The discussion conducted in Sections 24.2 and 24.3 suggests that: 1. different constitutive relations for the transmembrane current density model yield different results; and 2. a different approach to ion channel modeling results in a critical revision of the adopted assumptions of the simplified mathematical representation of ion channels. These considerations open up a number of fundamental questions in the mathematical modeling of problems in life sciences and bioengineering, for example:
24.4 CONCLUDING REMARKS
595
• To what extent can we trust simulation predictions? • Are we missing important biophysical mechanisms in the description of the problem? • What is the accuracy of the values of model parameters? Our goal here is not to provide an answer to the above questions. Rather, we merely invite the reader to go back to the founding concepts distilled in Chapter 1 and bear in mind that a genuine verification of model predictions can only be made by: 1. comparing model results with analytical solutions, whenever available; 2. comparing model results with experimental data; 3. modifying model parameters and/or the mathematical formulation, if the discrepancies between numerical results and experimental data are too large; 4. then, going back to item 1 of this list and continuing the process of model refinement until a satisfactory agreement is reached between model predictions and available data.
FINITE ELEMENT APPROXIMATIONS OF INITIAL VALUE/BOUNDARY VALUE PROBLEMS OF ADVECTION– DIFFUSION–REACTION TYPE
24
CHAPTER POINTS • In this chapter we consider the initial value/boundary value problem of advection–diffusion–reaction type representing the mathematical model of a balance law with a linear net production term.
• The numerical discretization of the initial value/boundary value problem is conducted using the θ-method for approximating the time derivative and the mixed GFEM for the spatial approximation.
• The stability and convergence properties of the numerical formulation are analyzed and computational examples are provided to illustrate the theoretical conclusions.
• The computational procedure is applied to the simulation of the propagation of an action potential along an axon and of ion electrodynamics in a transmembrane ion channel.
In Chapter 23 we have illustrated the Galerkin finite element formulation for the discretization of a boundary value problem with advective, diffusive, and reactive terms. In this chapter we consider the initial value/boundary value problem of advection–diffusion–reaction type that represents the mathematical model of a balance law with a linear net production term. We first approximate the time derivative of the initial value/boundary value problem using the θ -method introduced in Section 3.5.2.2, whereas the spatial approximation of the problem is conducted using the mixed GFEM introduced in Chapter 23. Then, we analyze the properties of stability and convergence of the numerical formulation and we provide a series of computational examples to support the theoretical conclusions. We conclude the chapter by applying the computational procedure to the simulation of two biophysical problems. The first problem is the propagation of an action potential along an axon using the cable equation model introduced in Chapter 18. The second problem is the electrodynamics of ions in a transmembrane ion channel using the Poisson–Nernst–Planck model introduced in Section 13.6.3. A Comprehensive Physically Based Approach to Modeling in Bioengineering and Life Sciences. https://doi.org/10.1016/B978-0-12-812518-2.00037-8 Copyright © 2019 Elsevier Inc. All rights reserved.
579
580
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
24.1 THE TIME-DEPENDENT ADVECTION–DIFFUSION–REACTION PROBLEM In this section we address the study of the one-dimensional advection–diffusion–reaction problem (22.1) with respect to the temporal coordinate. To simplify the presentation, we set γ∂ = 0, β∂ = 0, and α∂ = 1, corresponding to homogeneous Dirichlet boundary conditions at the endpoints of . Therefore, our considered initial value/boundary value problem reads as follows: Find u(x, t) for all (x, t) ∈ QT such that ∂u ∂ J + + κu = g ∂t ∂x ∂u J = uv − μ ∂x u(x, 0) = u0 (x) u(x, t) = 0
∀(x, t) ∈ QT ,
(24.1a)
∀(x, t) ∈ QT ,
(24.1b)
∀x ∈ , ∀(x, t) ∈ ∂ × (0, T ).
(24.1c) (24.1d)
As an example, the initial value/boundary value problem (24.1) can be regarded as the heat equation in the form of a balance law. In this physical picture, the dependent variable u has the meaning of the temperature of the medium in which the heat conduction process takes place, whereas J is the heat flux density. The initial value/boundary value problem (24.1) represents the mathematical model of the spatial and temporal diffusion and passive transport of heat in a medium of length L with thermal conductivity μ and heat transport velocity v, subject to a heat production g, a heat consumption rate κ, an initial thermal distribution u0 , and a boundary temperature equal to zero for all t ∈ (0, T ). Property 24.1 (Asymptotic stability). Set κ = 0 and g = 0 in (24.1a). Then we have lim u(x, t) = 0
t→∞
∀x ∈ .
(24.1e)
Relation (24.1e) expresses the fact that the initial value/boundary value problem (24.1) is asymptotically stable. Asymptotic stability of the initial value/boundary value problem is a consequence of the fact that μ is strictly positive and expresses the physical property that the system modeled by (24.1) is dissipative. Let us now use the θ -method introduced in Section 22.2 for the discretization of the initial value/boundary value problem (24.1) with respect to spatial and temporal variables. In the following, we denote by C a positive constant independent of h and t whose value may not be the same at each occurrence. Using the notation of Section 22.2 and Section 23.1.3, the linear algebraic system that has to be solved at each time level tk , k = 0, . . . , NT − 1, reads Aθ uk+1 = fk+1 θ ,
(24.1f)
1 M, t
(24.1g)
where Aθ = θ (Kdmix + Kamix + Krmix ) +
fk+1 = −(1 − θ )(Kdmix + Kamix + Krmix )uk + θ
1 Muk + θ gk+1 + (1 − θ )gk , t
(24.1h)
24.1 THE TIME-DEPENDENT ADVECTION–DIFFUSION–REACTION MODEL
581
having introduced the mass matrix M = ij
L
i, j = 1, . . . , Nh .
ϕj (x)ϕi (x) dx,
(24.1i)
0
The Matlab function 29.50 computes the mass matrix (24.1i) using the Cavalieri–Simpson quadrature formula (23.32c). We note that the mass matrix (24.1i) is a special instance of the reaction matrix (23.30c) corresponding to setting σ (x) = 1 for all x ∈ . Therefore, the lumping procedure illustrated in Section 23.4.1 for the diagonalization of the reaction matrix can be applied also to the mass matrix. This leads to = hδij , i, j = 1, . . . , Nh , (24.1j) M ij
where δij is the Kronecker symbol. The Matlab function 29.51 computes the lumped mass matrix (24.1j). In the above relations, θ is a parameter to be chosen in the interval [0, 1]. In particular, we note that: 1. If θ = 0, relation (24.1f) becomes an explicit formula to compute the solution at time level tk+1 using information only at the level tk . The corresponding scheme is called forward Euler or explicit Euler method. 2. If θ > 0, the family of schemes associated with (24.1f) takes the name of implicit methods because to compute uk+1 we need to solve a linear system using information at both levels tk and tk+1 . Definition 24.1 (Stability). Set κ = 0 and g = 0 in (24.1a). The θ -method (24.1f) is said to be stable (or also, asymptotically stable) if the sequence uk is such that lim uk = 0.
k→∞
(24.1k)
Condition (24.1k) is the discrete counterpart of (24.1e) and is not automatically inherited by the θ-method. In this sense, the following condition needs to be satisfied (for the proof, see [259, Chapter 11] and [258, Chapter 13]). Theorem 24.1 (Stability of the θ -method). If θ ≥ 1/2 the θ -method verifies (24.1k) for every value of t. In this case we say that the θ-method is unconditionally stable. If 0 ≤ θ < 1/2 the θ-method verifies (24.1k) only if t ≤
Ch2 . μ0 (1 − 2θ )
(24.1l)
In this case we say that the θ-method is conditionally stable. Conditional stability can be a severe computational constraint for the θ-method because for a small value of h (which is typically needed to achieve a reasonable spatial accuracy of the approximation) the corresponding value of t to make the θ -method stable turns out to be even smaller. Therefore, the use of an implicit method is in order thanks to its unconditional stability. The choice of such a method is governed by the next important result.
582
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
Theorem 24.2 (Convergence of the θ-method). Under suitable regularity assumptions on the data in (22.1), we have uk − ukh L2 () ≤ u0 − u0h L2 () + C(u0 , f, u) h2 + t p(θ) (24.1m) ∀k ∈ [1, NT ], where p(1/2) = 2 and p(θ) = 1 for all θ = 1/2. Theorem 24.2 tells us that: 1. if lim u0 − u0h L2 () = 0, then the θ -method is convergent; h→0
2. the θ -method is convergent with order 2 with respect to both spatial and time coordinates only in the case of the Crank–Nicolson scheme. All the other implicit schemes with θ > 1/2 converge linearly with respect to t. Quadratic convergence with respect to the spatial coordinate is a consequence of (23.91).
24.1.1 CONVERGENCE EXAMPLES IN THE TIME-DEPENDENT CASE In the following example we study the initial value/boundary value problem (24.1) using the θ-method in the case where μ = 1, v = 0, κ = 0, = (0, 1), IT = [0, 1], u0 (x) = sin(2πx), and g(x, t) = (4π 2 − λ)e−λt sin(2πx),
x ∈ , and t ∈ IT ,
with λ > 0, in such a way that the exact solution is u(x, t) = e−λt sin(2πx).
FIGURE 24.1 The time-dependent problem, convergence of the θ-method as a function of h with t = 1/100. The spatial mesh size is h = (10 · k)−1 , k = 1, . . . , 20; λ is equal to 1.
Figs. 24.1 and 24.2 illustrate the results obtained by solving the initial value/boundary value problem with the backward Euler and Crank–Nicolson methods. The parameter λ is set equal to 1. T Comparison between the two schemes is carried out by computing the quantity uNT − uN h L2 () . Fig. 24.1 is obtained by running the Matlab script 28.47 and performs a study of convergence with
24.2 SIMULATION OF THE PROPAGATION OF AN ACTION POTENTIAL
583
FIGURE 24.2 The time-dependent problem, convergence of the θ-method as a function of t with h = 1/500, the temporal mesh size is t = (10 · k)−1 , k = 1, . . . , 20.
respect to h having fixed t = 1/100, whereas the spatial grid size is allowed to vary in the interval [1/10, 1/200] by setting h = (10 · k)−1 , k = 1, . . . , 20. Conversely, Fig. 24.2 is obtained by running the Matlab script 28.48 and performs a study of convergence with respect to t having fixed h = 1/500, whereas the time grid size is allowed to vary in the interval [1/10, 1/200] by setting t = (10 · k)−1 , k = 1, . . . , 20. The two convergence studies conducted above clearly indicate that the Crank–Nicolson method is far more accurate than the backward Euler method with respect to both h and t. In particular, we note that, no matter how small h, the error of the backward Euler method in Fig. 24.1 is dominated by the contribution O( t). Since this latter quantity is fixed (equal to 0.01 in T the present case) the error norm uNT − uN h L2 () does not tend to zero as h tends to zero.
24.2 SIMULATION OF THE PROPAGATION OF AN ACTION POTENTIAL ALONG AN AXON In this section we numerically investigate the time-dependent GFEM illustrated in Section 24.1 applied to the simulation of the propagation of an action potential along an ideal axon of length L = 10 cm and radius a = 2.5 mm. The computational procedure consists of the use of (i) the functional iteration techniques illustrated in Section 21.5.1; (ii) the time discretization methods illustrated in the present chapter, and (iii) the mixed GFEM illustrated in Chapter 23.
24.2.1 MODEL DATA The axon membrane thickness is tm = 10 nm, the load output resistance Rout is set equal to +∞, so that the terminal of the axon at x = L is, from the electrical standpoint, an open circuit and no electric current is allowed to flow out of it. (See Table 24.1.)
584
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
Table 24.1 Values of model parameters in the simulation of the propagation of an action potential along an axon using the cable equation model. The temperature of the system under investigation is = 298.15 K. Transmembrane conduction current is assumed to be driven by potassium, sodium, and chloride ions. Model parameter
Description
Value
r εm t0 , t1 , t2 , t3 Tf in Rout a L tm cm V1 , V2
Dielectric relative constant Time characteristic of the axon potential Final simulation time Load output resistance Radius of the ideal axon Length of the ideal axon Membrane thickness of the ideal axon Unit area capacitance of the ideal axon Electric characteristic of the axon potential
80 [·] 0,5,7.5,10 ms 15 ms +∞ 2.5 · 10−3 m 0.1 m 10−8 m 10−2 F m−2 28.73, −91.27 mV
Table 24.2 Values of intracellular and extracellular ion concentrations in the study of propagation of an action potential along an axon in the human body (cf. [107, Table 1.1].) Model parameter Description K+
(in) CK (in) CNa (in) CCl CK(out) (out) CNa (out) CCl
concentration in the intracellular region
Value 150 mM
Na+ concentration in the intracellular region
50 mM
Cl−
10 mM
concentration in the intracellular region
K+ concentration in the extracellular region Na+
5 mM
concentration in the extracellular region 460 mM
Cl− concentration in the extracellular region
125 mM
The resting potential of the cell can be computed by using the Goldman equation (17.24b) in which the ion permeabilities PK , PNa , and PCl satisfy relation (17.39b). Setting T = 300 K and using the values of the intracellular and extracellular concentrations listed in Table 24.2, we obtain ψeq = −71.27 mV.
(24.2)
24.2.2 ACTION POTENTIAL The input voltage stimulus V in is the piecewise step function illustrated in Fig. 24.3. This function is a schematic representation of the action potential which develops in an excitable cell such as a neuron or a cardiac myocyte. The applied voltage is equal to the resting potential of the cell ψeq from time t = t0 = 0 ms until time t1 = 5 ms. In these conditions, the intracellular potential is lower than the extracellular potential so that the “−” pole is the intracellular site and the “+” pole is the extracellular site. At time t = t1 , the applied voltage is suddenly raised to the value V1 = ψeq + 100 mV = +28.73 mV and is kept equal to V1 until time t2 = 7.5 ms. In these conditions, a change of polarity takes place because the intracellular potential is higher than the extracellular potential so that the “+” pole is the intracellular
24.2 SIMULATION OF THE PROPAGATION OF AN ACTION POTENTIAL
585
FIGURE 24.3 Input voltage stimulus. The piecewise constant function shown in the figure schematically represents an action potential. The jump between the resting potential to the positive value V1 = 28.73 mV corresponds to the depolarization phase. The jump from V1 to the value V2 = −91.27 mV corresponds to the hyperpolarization phase.
site and the “−” pole is the extracellular site. This change of polarity is the result of a flow of a positive ion charge, such as Na+ , that was attracted to the cytoplasm by the “+” pole in the resting condition. This phase of the action potential is referred to as depolarization (see [107, Chapter 1]). At time t = t2 , the voltage stimulus is again suddenly changed to the value V2 = ψeq − 20 mV = −91.27 mV and is kept equal to V2 until time t3 = 10 ms. In these conditions, another change of polarity takes place because the intracellular potential is more negative than the resting potential so that the “−” pole is the intracellular site and the “+” pole is the extracellular site. This change of polarity is the result of a flow of a negative ion charge, such as Cl− , that was attracted to the cytoplasm by the “+” pole in the depolarization condition. This phase of the action potential is referred to as hyperpolarization (see [107, Chapter 1]). At time t = t3 , the voltage stimulus is again suddenly changed to the value of the resting potential ψeq and is kept equal to ψeq until the final time t = Tf in = 15 ms.
24.2.3 OVERVIEW In Chapter 17 we illustrated two different models to characterize the transmembrane current, namely, the linear resistor formulation (17.44) and the nonlinear Goldman–Hodgkin–Katz (GHK) formulation (17.51). In the following paragraphs, we utilize the functional iteration (21.16) and its time and spatial discretization to compare the solutions obtained using the two models. The numerical implementation of the computational algorithm is performed by the Matlab script 28.49. The Matlab function 29.28 implements the solution of the problem using the linear resistor formulation, whereas the Matlab function 29.29 implements the solution of the problem using the GHK formulation. In both functions, the time discretization of (18.2) is performed using the backward Euler method illustrated in Section 3.5.2.1, whereas the spatial discretization is performed using the mixed GFEM illustrated in Chapter 23. A uniform time step t = (Tend − t0 )/250 = 6 · 10−5 s and a uniform spatial discretization parameter h = L/100 = 0.1 cm are used in the numerical computations. The Matlab function 29.27 computes the voltage stimulus graphically represented in Fig. 24.3.
586
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
Simulation With the Linear Resistor Model In this paragraph we perform the numerical simulation of the propagation of the action potential elicited by the voltage stimulus of Fig. 24.3, by representing the axon transmembrane current density through the linear resistor model (17.44). We consider three values of the longitudinal resistivity, and we set SW (case (a)), ρ = 10+3 ρ SW (case (b)), and ρ = 10−3 ρ SW (case (c)). To use the linear ρin = ρin in in in in resistor model, we need to determine the ion specific conductance gα (unit: S m−2 ), α ∈ {K, Na, Cl}. With this aim, we use the definition of local ion conductivity (13.75d) in the following average form: σα = q|zα |μel α
nα (0) + nα (tm ) nin + nout α = q|zα |2 Dα α , 2tm 2tm
α ∈ {K, Na, Cl} ,
(24.3)
out where nin α and nα are the intra and extracellular number densities of the αth ionic species and Einstein’s relation (13.69) is used to express the ion electric mobility μel α in favor of the ion diffusion coefficient Dα . Results are shown in Figs. 24.4, 24.5, and 24.6.
FIGURE 24.4 Simulation of axon potential evolution using the cable equation model (18.2) supplied by the linear resistor model SW = 2 · 10−1 m, for the transmembrane current. In this simulation the axoplasm resistivity ρin is set equal to ρin −1 corresponding to an axoplasm conductivity σin = 5 S m . Left panel: view from x = L. Right panel: view from x = 0.
Case (a). Fig. 24.4 shows the axon membrane potential ψm (x, t) for x ∈ [0, L] and for t ∈ [0, Tend ], computed by solving (21.16) with the linear resistor model for the axon transmembrane current density. SW = 2 · 10−1 m, the resistivity of salt The value of the axoplasm resistivity ρin is set equal to ρin ◦ water at 20 C. This value is very close to the data shown in [83], where the experimentally measured SW . Results indicate the resistivity of the extruded squid axon turned out to be equal to (1.4 ± 0.2)ρin following: • At t = 0, the membrane potential undergoes a fast increase in the neighborhood of x = 0 and reaches a constant value close to 0 mV along the axon length because of a flux of positive charge inside the axon.
24.2 SIMULATION OF THE PROPAGATION OF AN ACTION POTENTIAL
587
FIGURE 24.5 Simulation of axon potential evolution using the cable equation model (18.2) supplied by the linear resistor model SW = 2 · 102 m, for the transmembrane current. In this simulation the axoplasm resistivity ρin is set equal to 103 ρin corresponding to an axoplasm conductivity σin = 5 · 10−3 S m−1 . Left panel: view from x = L. Right panel: view from x = 0.
FIGURE 24.6 Simulation of axon potential evolution using the cable equation model (18.2) supplied by the linear resistor SW = model for the transmembrane current. In this simulation the axoplasm resistivity ρin is set equal to 10−3 ρin −4 3 −1 2 · 10 m, corresponding to an axoplasm conductivity σin = 5 · 10 S m . For each t ∈ [0, Tf in ], the spatial distribution of the membrane potential along the axon is constant because of the elevated longitudinal conductivity.
• As t increases, the spatial distribution of the membrane potential remains uniform but with a smaller value compared to t = 0 because of a flux of part of the accumulated positive charge out of the intracellular axon region. • At t = t1 , the membrane potential rises abruptly at x = 0 and then monotonically decreases along the axon length, passing from positive values to negative values, because of a positive accumulated charge inside the axon near x = 0 and a negative accumulated charge for larger values of x. The membrane potential exhibits an opposite behavior at t = t2 .
588
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
• At t = t3 , the membrane potential returns abruptly to its resting value, which is less negative than V2 . This induces an accumulation of positive charge inside the intracellular region of the axon and the membrane potential becomes uniform along the axon length, reaching a value of about −55 mV. • For larger values of t, the accumulated positive charge flows out of the axon, thereby giving rise to more negative values of the membrane potential. After a certain time instant, all the accumulated positive charge has flown out of the axon, so that the membrane potential becomes flat along the whole axon length and equal to the boundary condition at x = 0. • For all t ∈ [0, Tf in ], the spatial distribution of the membrane potential becomes flat at the end of the axon because of the no-flux boundary condition at x = L. Case (b). Fig. 24.5 shows the axon membrane potential computed by solving (21.16) with the linear resistor model and with an axoplasm resistivity increased by a factor of 103 with respect to the previous case. Results indicate that the membrane potential undergoes a practically linear decrease with respect to time for each x ∈ [0, L]. This behavior is due to the elevated value of the longitudinal axon resistivity which considerably slows down the transmembrane flow of positive and negative charge for each x ∈ [0, L]. We also note that the axon is far from resting conditions at t = Tf in because the membrane potential is about −50 mV. Resting conditions would be reached only for a significantly larger value of the observational time Tf in . Case (c). Fig. 24.6 shows the axon membrane potential ψm (x, t) for x ∈ [0, L] and for t ∈ [0, Tend ] in SW = the case where the axoplasm resistivity is reduced by a factor of 103 , so that we have ρin = 10−8 ρin −9 2 · 10 m. Results indicate that for each value of t ∈ [0, Tend ] the spatial distribution of the membrane potential along the axon is constant because the elevated longitudinal conductivity gives rise to an almost instantaneous transport of the value of the Dirichlet boundary condition ψm (0, t) = V in (t).
Simulation With the GHK Model In this paragraph we perform the numerical simulation of the propagation of the action potential elicited by the voltage stimulus of Fig. 24.3, by representing the axon transmembrane current density through the GHK model (17.51). We consider the same conditions as in Paragraph ‘Simulation With the Linear SW (case (a)), ρ = 10+3 ρ SW (case (b)), and ρ = 10−3 ρ SW Resistor Model’, that is, we set ρin = ρin in in in in (case (c)). Computations show that four iterations were needed to satisfy the convergence check (21.17) with a tolerance toll = 10−5 . Results are shown in Figs. 24.7, 24.8, and 24.9. Model predictions in case (c) cannot be distinguished from the predictions given by the linear resistor model in the same conditions. Model predictions in cases (a) and (b) are qualitatively similar to the corresponding simulations with the linear resistor model. In case (a), we see that: • the increase of membrane potential at t = 0 in the neighborhood of x = 0 is less sharp than predicted by the linear resistor model; • the membrane potential at x = L is less negative than predicted by the linear resistor model; • the membrane potential distribution at t = Tf in is not as flat as predicted by the linear resistor model. This indicates that, according to the GHK model, the axon has not reached resting conditions as in the case of the linear resistor model. In particular, the value of the membrane potential predicted by the GHK model at x = L is less negative than that predicted by the linear resistor model. Similar considerations apply for case (b). In particular, we see that:
24.2 SIMULATION OF THE PROPAGATION OF AN ACTION POTENTIAL
589
FIGURE 24.7 Simulation of axon potential evolution using the cable equation model (18.2) supplied by the GHK model for the SW = 2 · 10−1 m, corretransmembrane current. In this simulation the axoplasm resistivity ρin is set equal to ρin sponding to an axoplasm conductivity σin = 5 S m−1 . Left panel: view from x = L. Right panel: view from x = 0.
FIGURE 24.8 Simulation of axon potential evolution using the cable equation model (18.2) supplied by the GHK model for the SW = 2 · 102 m, corretransmembrane current. In this simulation the axoplasm resistivity ρin is set equal to 103 ρin −3 −1 sponding to an axoplasm conductivity σin = 5 · 10 S m . Left panel: view from x = L. Right panel: view from x = 0.
• for each x ∈ [0, L] the decay in time of the membrane potential is significantly smaller than predicted by the linear resistor model; • the resulting membrane potential at x = L predicted by the GHK model is much less negative than in the case of the linear resistor model, so that a much larger observational time is needed to reach resting conditions in the axon.
590
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
FIGURE 24.9 Simulation of axon potential evolution using the cable equation model (18.2) supplied by the GHK model for the SW = 2 · 10−4 m, transmembrane current. In this simulation the axoplasm resistivity ρin is set equal to 10−3 ρin 3 −1 corresponding to an axoplasm conductivity σin = 5 · 10 S m . For each t ∈ [0, Tf in ], the spatial distribution of the membrane potential along the axon is constant because of the elevated longitudinal conductivity.
24.2.4 COMPARISON AND CONCLUDING REMARKS The same considerations developed in Section 17.7.4.3 can be applied to discuss the quality and reliability of model predictions obtained in Paragraphs ‘Simulation With the Linear Resistor Model’ and ‘Simulation With the GHK Model’ using the linear resistor model and the GHK model, respectively. In particular, the accuracy of the computed action potential profiles may be warranted through the following two approaches: • an error analysis against an exact solution, and/or • comparison with experimental data, whenever available. As far as the first approach, we refer to [223], whereas for the second approach we refer, e.g., to [317].
24.3 POISSON–NERNST–PLANCK SIMULATION OF AN IDEAL ONE-DIMENSIONAL ION CHANNEL In this section we apply the Poisson–Nernst–Planck model illustrated in Section 13.6.3 to the simulation of an ideal one-dimensional ion channel of length L = 10 · 10−7 cm in which only the Ca++ and Cl− ions are considered. We warn the reader that the scope of this section is not the mathematical modeling of a realistic ion channel as addressed in Chapter 26, but, rather, the numerical investigation of the time-dependent GFEM (see Section 24.1) and the application of the functional iteration (21.22). To devise the computational procedure, we apply (i) the functional iteration techniques illustrated in Section 21.5.2, (ii) the time discretization methods illustrated in the present chapter, and (iii) the mixed GFEM illustrated in Chapter 23. All the computations reported in this section are conducted by running the Matlab code bio_1d, which can be downloaded at [210]. Fig. 24.10 shows the scheme of the one-dimensional ion channel used in the simulation of the Poisson–Nernst–Planck model.
24.3 POISSON–NERNST–PLANCK SIMULATION OF AN IDEAL
591
FIGURE 24.10 Scheme of the one-dimensional ion channel.
24.3.1 MODEL DATA In the ion channel represented in Fig. 24.10, Ca++ and Cl− ions flow under the action of diffusion and electric forces (see Section 16.2). The channel is represented by a medium having the same dielectric constant of water, so that in the Poisson equation (13.74c) we can set εm = ε0 εr,w , where ε0 is the dielectric permittivity of vacuum and εr,w = 80 is the relative dielectric permittivity of water. The presence of a surface charge density that characterizes the structure of any realistic ion channel is replaced by a volumetric negative charge density ρf ixed = −9.6 · 1018 cm−3 . Ion transport parameters such as diffusivity, electrical mobility, and electric charge are reported in Table 24.3. The simulation time spans the interval IT = [tinitial , tf inal ] with tinitial = 1 · 10−10 s and tf inal = 10 · 10−7 s and the temporal discretization is obtained by dividing IT l into a uniform partition of NT = 200 elements, so that the temporal step size is t = (tf inal − tinitial )/NT . The spatial partition Th used by the GFEM consists of a uniform subdivision of the channel length into Mh = 500 elements, so that h = L/Mh . Table 24.3 Parameters used in the simulation of an ideal one-dimensional ion chan-
nel in which Ca++ and Cl− ions are flowing. The temperature of the system under investigation is = 298.15 K.
Parameter Description
Value
r εm tinitial tf inal L ρf ixed zCa++ DCa++ μCa++ zCl− DCl− μCl− Mh NT
80 [·] 1 · 10−10 s 1 · 10−7 s 10 · 10−7 cm −9.6 · 1018 cm−3 +2 · 1.6019 · 10−19 C 1 · 10−6 cm2 s−1 7.735 · 10−5 cm2 V−1 s−1 −1.6019 · 10−19 C 0 1 · 10−6 cm2 s−1 3.868 · 10−5 cm2 V−1 s−1 500 [·] 200 [·]
Dielectric relative constant Initial simulation time Final simulation time Length of the ideal one-dimensional channel Fixed charge concentration in the channel Electric charge of the Ca++ ion Diffusivity of the Ca++ ion Electric mobility of the Ca++ ion Electric charge of the Cl− ion Diffusivity of the Cl− ion Electric mobility of the Cl− ion Number of elements in the finite element partioning Number of time intervals used in the θ-method
The dynamic behavior of the ion channel is investigated using the time-dependent Poisson–Nernst– Planck system (13.74). In the present case we have Mion = 2, so that the Poisson–Nernst–Planck model is constituted by five equations: the Poisson equation (13.74c), two Eqs. (13.74a) expressing the balance of mass for each ion, and two Eqs. (13.74b) expressing the balance of linear momentum for
592
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
each ion. Dirichlet boundary conditions are applied at both endpoints of the channel to Eqs. (13.74c) and (13.74a) and are assumed to be constant in time. Table 24.4 reports the data for the boundary and initial conditions used in the simulations. The channel is subjected to a bias condition similar to that considered in Chapter 26 by applying at x = L a negative voltage of −2.5 · 10−3 V and at x = 0 a voltage equal to 0 V, in such a way that the externally applied membrane potential is ψm = 0 − (−2.5 · 10−3 ) = 2.5 · 10−3 V. A concentration gradient is forced between the internal and external cell regions for both ions but in opposite directions. This gradient is obtained by using similar biophysical values for the ion concentrations at x = 0 and x = L as calculated in Section 26.2. Table 24.4 Initial and boundary conditions applied to the simulation of an ideal onedimensional ion channel in which Ca++ and Cl− ions are flowing.
Initial or boundary condition Description ψ0 ψL ψ(t0 ) Ca++ 0 Ca++ L Ca++ (t0 ) Cl− 0 Cl− L Cl− (t0 )
Electric potential at x = 0 and in [tinitial , tf inal ] Electric potential at x = L and in [tinitial , tf inal ] Initial value of the electric potential in = [0, L] Ca++ concentration at x = 0 and in [tinitial , tf inal ] Ca++ concentration at x = L and in [tinitial , tf inal ] Initial value of the Ca++ concentration in = [0, L] Cl− concentration at x = 0 and in [tinitial , tf inal ] Cl− concentration at x = L and in [tinitial , tf inal ] Initial value of the Cl− concentration in = [0, L]
Value 0V −2.5 · 10−3 V 0V 4 mM 1.33 mM 1.33 mM 107 mM 131 mM 107 mM
24.3.2 SIMULATION RESULTS In this section we report the simulation results of ion transport in the one-dimensional channel described above. Despite the geometric simplicity of the considered one-dimensional channel, the test reveals some interesting model implications and numerical observations. Fig. 24.11 reports the convergence behavior of the Gummel map within the functional iteration (21.22) obtained by solving the Poisson–Nernst–Planck model (13.74) up to a tolerance toll = 10−8 for each discretization time step (cf. the termination condition, Eq. (21.27)). Except for the first time step, in the remaining time intervals the number of Gummel iterations to reach convergence up to the given tolerance remains constant. The evolution of the electric quantities is reported in Fig. 24.12, where the black dashed lines refer to the profiles computed at the final simulation time t = tf inal . The left panel shows the electric potential evolution; the values at the domain endpoints are kept constant thanks to the Dirichlet boundary conditions while within the channel the electric potential is found to have always a nonlinear shape. As a consequence, the electric field is not constant as clearly visualized in the right panel of Fig. 24.12; in the first part of the channel, the electric field is directed to the right (extracellular region) while in the final part of the channel the electric field is directed to the left (intracellular region). We note that in the center of the channel the electric field strength is much lower than in the regions close to the boundary. This finding is in contrast with the assumption of a constant electric field made in Section 17.6.2 to derive the GHK model. Thus, we see that the present investigation shows how the assumptions made to derive the GHK model to calculate the transmembrane current may fail to be verified in all condi-
24.3 POISSON–NERNST–PLANCK SIMULATION OF AN IDEAL
593
FIGURE 24.11 Number of iterations required by the Gummel map to reach convergence in the functional iteration (21.22) up to a tolerance toll = 10−8 for each discretization time step.
FIGURE 24.12 Simulation of the evolution of electric quantities in an ideal one-dimensional ion channel of length L = 15 · 10−7 cm in which only Ca++ and Cl− ions are considered. Results have been obtained with the Poisson– Nernst–Planck model (13.74). The black dashed lines refer to the profiles at the final simulation time t = tf inal . (Initial conditions are omitted). Left panel: electric potential. Right panel: electric field.
tions. This is a general indication, and, as a consequence, the application of any mathematical model to describe a given experimental condition needs always to be carefully verified. Fig. 24.13 reports the evolution of ion concentrations. The black dashed lines refer to the profiles computed at the final simulation time t = tf inal . The evolution of the Ca++ concentration is shown in the left panel. Looking at the electric field close to the left boundary of the channel (x = 0), we see that the drift component of the Ca++ current density is directed towards the right boundary (x = L) because of the positive sign of the electric field. Vice versa, close to the right boundary (x = L) the electric field is negative and pushes the Ca++ ions away from it. This combination of effects results in an accumulation of Ca++ ions within the center of the first half of the channel.
594
CHAPTER 24 FINITE ELEMENT APPROXIMATIONS OF MODELS
The evolution of the Cl− concentration is shown in the right panel of Fig. 24.13. Due to the negative sign of the Cl− ion charge, the electric field pushes the chloride ions towards the boundaries (at x = 0 and x = L), giving rise to the formation of two pronounced boundary layers and reducing the chloride concentration in the center of the channel region.
FIGURE 24.13 Simulation of the evolution of ion concentrations in an ideal one-dimensional ion channel of length L = 15 · 10−7 cm in which only the Ca++ and Cl− ions are considered. Results have been obtained with the Poisson– Nernst–Planck model (13.74). The black dashed lines refer to the profiles at the final simulation time t = tf inal (initial conditions are omitted). Left panel: Ca++ concentration. Right panel: Cl− concentration.
Among the several indications we have obtained in the discussion of this section, we can conclude that the adoption of a more sophisticated model than the lumped formulations illustrated in Chapter 17 has the disadvantage that it requires advanced mathematical methods for its solution, but, at the same time, it may help the user to better investigate the underlying biophysical mechanisms and verify a posteriori the correctness of the assumptions made to derive the reduced-order model approach.
24.4 CONCLUDING REMARKS BIO-WARNING 24.1 The discussion conducted in Sections 24.2 and 24.3 suggests that: 1. different constitutive relations for the transmembrane current density model yield different results; and 2. a different approach to ion channel modeling results in a critical revision of the adopted assumptions of the simplified mathematical representation of ion channels. These considerations open up a number of fundamental questions in the mathematical modeling of problems in life sciences and bioengineering, for example:
24.4 CONCLUDING REMARKS
595
• To what extent can we trust simulation predictions? • Are we missing important biophysical mechanisms in the description of the problem? • What is the accuracy of the values of model parameters? Our goal here is not to provide an answer to the above questions. Rather, we merely invite the reader to go back to the founding concepts distilled in Chapter 1 and bear in mind that a genuine verification of model predictions can only be made by: 1. comparing model results with analytical solutions, whenever available; 2. comparing model results with experimental data; 3. modifying model parameters and/or the mathematical formulation, if the discrepancies between numerical results and experimental data are too large; 4. then, going back to item 1 of this list and continuing the process of model refinement until a satisfactory agreement is reached between model predictions and available data.