Derivation of cable equation by multiscale analysis for a model of myelinated axons

The paper concerns the multiscale modeling of a myelinated axon. Taking into account the microstructure with alternating myelinated parts and nodes Ranvier, we derive a nonlinear cable equation describing the potential propagation along the axon. We assume that the myelin is not a perfect insulator, and assign a low (asymptotically vanishing) conductivity in the myelin. Compared with the case when myelin is assumed to have zero conductivity, an additional potential arises in the limit equation. The coefficient in front of the effective potential contains information about the geometry of the myelinated parts.

however, several grouped fascicles, each of them containing many axons. The jump of the potential across the membrane of each individual axon can be modelled in the framework of the Hodgkin-Huxley model, but the alternating myelinated and unmyelinated parts of the membrane present an obvious problem for those attempting to describe its macroscopic response to the electrical stimulation. In order to model and simulate the respons of biological tissues to electrical stimulation one needs to know how signals propagate along single neurons and, as the next step, how they influence each other in a bundle of axons.
The signal propagation along a neuron is modelled by a cable equation, usually derived by modeling dendrites and axons as cylinders composed of segments with capacitances and resistances combined in parallel ( [1], [8], [9], [11], [10]). The coefficients in such equation depend on the membrane resistances and capacitance of Ranvier nodes and internodes (myelinated parts), as well as on the length of nodes and internodes. There are several works where formal two-scale expansion is applied to a one-dimensional model in order to show that a myelinated neuron can be approximated by a homogeneous cable ( [11], [12]), but these results do not take into account the microstructure of the fibers, and the geometry of the myelin sheath in particular, as well as they do not justify the formal approximation.
There are many results where the homogenization is applied to cardiac tissue: [13], [14], [15], [16] Cardiac muscle is however fundamentally different from nerve tissue because the heart is a syncytium. The intracellular space of each cardiac cell is coupled to its neighbor's through intercellular channels. Thus, current can flow from the interior of one cell to the interior of another without crossing a cell membrane.
The present work presents a rigorous derivation of a nonlinear cable equation for signal propagation along a myelinated neuron. We assume that the conductivity of the myelin sheath is small, but not zero, that leads to the appearance of a potential in the limit equation. The potential depends on the geometry of the myelin sheath.
The paper is organized as follows. In Section 2 we formulate the problem and present the main result in Theorem 2.1. The rest of the paper is devoted to the proof of Theorem 2.1. In Sections 3 we derive a priori estimates for the potential u ε and its jump across the Ranvier nodes. In Section 4 we construct an auxiliary test function which is used when passing to the limit in Section 5.

Problem setup
Let us consider a myelinated axon sparsely suspended in an extracellular medium. We assume that the axon has a periodic structure, containing myelinated and unmyelinated parts (nodes of Ranvier) as illustrated on Figure 1. Figure 1. Simplified geometry of the cross-section of a myelinated axon and the periodicity cell Y .
A periodicity cell will be denoted by Y = (− 1 2 , 1 2 ) × D R0 where D R0 is the disk in R 2 with the radius R 0 (see Figure 1). Y consists of an intracellular part Y i = (0, 1)×D r0 , an extracellular medium Y e , and the myelin sheath Y m as shown in Fugure 1 (a detailed description of the domain is given in Section 4). We denote by Γ mi (Γ me ) the interface between Y m and Y i (Y e ). Γ m = Γ mi ∪Γ me is the myelinated part of the interface, and Γ is the unmyelinated one (surface of a Ranvier node). The lateral boundary of Y is denoted by Σ (we will assume periodicity in y 1 ). We assume that the boundary of the myelin part Γ m is Lipschitz continuous. The periodicity cell is then scaled by a small parameter ε > 0 and translated along the x 1 -axis to form a thin periodic cylinder (thickness of order ε) suspended in the extracellular medium (thickness of order ε) with alternating myelinated and unmyelinated parts on the lateral boundary.
In what follows we denote Let Ω i,ε = (0, L) × (εD r0 ) denote the intracellular domain, Ω e,ε denote the extracellular domain, Ω m,ε denote the myelin part, Γ ε be the unmyelinated part of the boundary, and Γ m,ε be the myelinated one. For simplicity L consists of integer number of periods.
The whole domain Ω ε = (0, L) × ε(− 1 2 , 1 2 ) 2 is the union of the extracellular, intracellular and myelin domains, and the Ranvier nodes: Let u i ε , u e ε , u m ε denote the electrical potential in the intracellular, extracellular and myelin domains, respectively. We assume that the electric potential satisfies homogeneous Neuman boundary conditions on the lateral boundary Σ ε and homogeneous Dirichlet boundary conditions on the bases . The transmembrane potential is the jump of the potential through the axon's membrane. We denote it by [u ε ] = u i ε − u e ε . Let the conductivity be a piecewise constant function and u ε denote the potential u ε = u l ε in Ω i,ε , l = i, e, m. The potential distribution in Ω ε is described by the following system of equations: We study the asymptotic behavior of u ε , as ε → 0, and derive a one-dimensional effective equation describing the action potential propagation along the axon.
On the Ranvier nodes we assume the continuity of currents (2), and the Hodgkin-Huxley dynamics for the transmembrane potential (3). Following the Hodgkin-Huxley model, the applied current through the membrane is a sum of the capacitive current c m ∂ t [u ε ], where c m is the membrane capacitance per unit area, and the ionic current I ion ([u ε ], g ε ) through the ion channels. In the classical Hodgkin-Huxley model there are three types of channel: a sodium channel (Na), a potassium channel (K), and a leakage channel. The conductances of the various ionic fluxes are regulated by the vector of gating variables g ε .
We assume the homogeneous Dirichlet boundary condition for u e ε and for u i ε on the bases of the domain, when x 1 = 0 and x 1 = L; on the lateral boundary of Ω ε we assume the homogeneous Neumann boundary condition; ν is the unit normal exterior to Ω e,ε on Σ and Γ me , and exterior to Ω i,ε on Γ ε and Γ mi . Note that ν on Γ is orthogonal to the x 1 -axes, that is its first component is zero.
We assume that (H1) The function I ion (v, g) is linear w.r.t v and has a form where g ,j is the jth component of g, v r,j is the jthe component of the resting potential v r , and H j is positive, bounded, and Lipschitz continuous The constant v r is the reference constant voltage, and g ε is a gate variable vector with positive components 0 < (g ε ) j < 1, j = 1, m.
Remark 1. When measuring the respons of a neuron to the external stimulation, one wants to exclude appearance of the action potential in the absence of the external stimulation. To this end one can control the initial state of ionic channels (initial condition for the gate variables) in order to guarantee zero potential at the initial moment. This motivates the choice of zero initial condition for the transmembrane potential v ε .
We will use test function φ ∈ L ∞ (0, T ; H 1 (Ω ε \ Γ ε )), ∂ t φ ∈ L 2 (0, T ; L 2 (Γ ε )) such that φ = 0 for x 1 = 0 and x 1 = L. The jump of φ across the Ranvier nodes is denoted by The weak formulation corresponding to (1-7) is given by: Find such that u ε = 0 for x 1 = 0 and x 1 = L, for any test functions φ ∈ L ∞ (0, T ; H 1 (Ω ε \ Γ ε )), φ = 0 for x 1 = 0 and x 1 = L, and for almost all t ∈ (0, T ) The vector of gate variables g ε solves the following ordinary differential equation Since HH is linear with respect to g ε , we can solve the last ODE and obtain g ε as a function (integral functional) of the jump [u ε ]: Substituting this expression into (15) we obtain the weak formulation of (1)-(7) in terms of the potential u ε and its jump v ε = [u ε ] across Γ ε : The main result of the paper is given in the following theorem.
Theorem 2.1. The transmembrane potential [u ε ] and the vector of gating variables g ε converge uniformly with respect to t in C(0, T ; L 2 (Γ ε )) to the unique solution (v 0 , g 0 ) of the following one-dimensional problem: The effective coefficient a eff is given by where the 1-periodic in y 1 function N solves an auxiliary cell problem The constant Λ depends on the geometry of the myelin sheath (see Figure 4) and the conductivities, and is given by Remark 2. The effective coefficient a eff can be interpreted as the conductivity of the bulk medium corresponding to the conductivity of the intra-and extracellular domains connected in series. The effective potential Λ is a decreasing function of the angles ϕ A , ϕ B and it goes to zero when the angles approach π.
Remark 3. Note that, since the equation for g 0 is linear in g 0 , we can solve it explicitly ) is also a Lipschitz function. In this way the effective problem is one nonlinear diffusion equation v 0 (0, x 1 ) = 0, To prove Theorem 2.1 we first derive a priori estimates in Section 3 (Lemma 3.2), then we prove the two-scale convergence of u ε and its gradient (Lemma 3.5) and the convergence of [u ε ] in appropriate spaces (Lemma 3.6). Finally, in Section 5 we pass to the limit in the weak formulation and derive the limit problem (10). Section 4 is devoted to the construction of an auxiliary function, the main ingredient of the test function used when passing to the limit in the weak formulation.

A priori estimates
Proof. The existence of a mild solution follows from the classical semigroup theory (see, for example, [5]). For the existence of more regular solutions see [4], [3]. We present just an idea of the proof.
Denote v ε = [u ε ] and let us rewrite (1)- (7) in the form where the operator A : into the solution u ε and then to the normal derivative σ ε ∇u ε · ν. To construct such an operator we fix is closed and densely defined. Due to the Poincaré inequality, the quadratic form is negative and thus the resolvent set of A contains R + . Futherrmore, for λ > 0 and v L 2 (Γε) = 1 we have that implies that A is the infinitesimal generator of a strongly continuous semigroup of contractions (see Theorem 3.1 in [5]). Since ). This is done by deriving a priori estimates as in Lemma 3.2.
Then the following estimates hold: Proof. Let us multiply (1) by u ε , integrate by parts over Ω ε \ Γ ε and divide the resulting identity by ε 2 (the scaling factor of the order of measure of the thin domain Ω ε ): Integrating the last equality with respect to t we get Using Lipschitz continuity of I ion and applying Grönwall's inequality we obtain the following estimate for v ε : for some constants C, C 1 and γ 0 > 0 in dependent of ε. Estimate (i) is proved.
From (18) and (i) we derive an integral estimate for ∇u ε : Let us now multiply (1) by ∂ t u ε and integrate by parts over Ω ε \ Γ ε : To find ∇w ε (x) = ∇u ε t=0 we solve the following elliptic problem It is clear that ∇w ε = 0. The Grönwalls inequality applied in (19) yields (ii). Estimates (19) and (ii) imply that Since u ε satisfies the homogeneous Dirichlet boundary condition for x 1 = 0, Friedrichs's inequality is valid for u ε in Ω i,ε and Ω e,ε which gives us (ii). In order to obtain an L 2 -bound for u ε in Ω m,ε we use the Poincaré inequality inequality in each myelin part εY m,k and then sum them up to obtain an estimate in ∪ k εY m,k = Ω m,ε . Namely, for u ∈ H 1 (εY m,k ), letū m ε,k denote the mean value over kthe interface between the intracellular domain and myelin Γ mi,k u m ε,k = 1 |εΓ mi | εΓ mi,k u ds.
We derive with the help of the Poincarï¿ 1 2 inequality Due to the continuity of traces of u ε , Combining (21)-(23) we obtain Adding up εY m,j and taking into account (20) yields the estimate for the which completes the proof.
Proof of Lemma 3.6. Let us cover Ω ε into a union of overlapping cells εY k as depicted in Figure 2. We recall that Γ εY k εY k+1 εY e,k εY e,k+1 εΓ k Figure 2. Overlapping cells Y k covering Ω ε .
We start by estimating the difference between the mean values of [u ε ] over εΓ k and εΓ k+1 . Letū l ε,k = 1 |εΓ| εΓ k u l ε ds, l = i, e.
For each εY l,k , l = i, e, we have owing to the Poicarï¿ 1 2 inequality, with C independent of ε. Considering traces on Γ k by simple scaling argument one has Then the difference between two averagesū ε,k andū ε,k+1 is estimated as follows Adding up in k the above estimates we obtain an estimate in Ω l,ε : Introduce the following notation Then (24) and (25) yield Bounds (26) show that [u ε ] in each cell εY k is close to a constantv ε,k , and the difference betweenv ε,k andv ε,k+1 is small. Now we construct a piecewise linear functionṽ ε (t, x 1 ) interpolating valuesv ε,k linearly and show that Indeed, (27), (28) follow directly from (26): Estimate (28) is proved in a similar way using (26): Let us prove (29). Differentiatingv ε,k with respect to t, using the the Cauchy-Schwarz inequality yields Similarly to (30), estimate (29) follows from the last bound and (ii) in Lemma 3.2. Estimate (i) in the current lemma follows from (26). The uniform convergence on (0, T ) of the constructed piecewise linear approximation is given by the Arzelà-Ascoli theorem.
Theorem 3.7 (Arzelà-Ascoli theorem). Let (X, d) be a compact metric space. Then a set F ⊂ C 0 (X; E) is precompact (any sequence has a converging subsequence converging uniformly in X to f ∈ C 0 (0, T ; E), not necessarily in F) provided (1) F(x) in precompact in E, for each x ∈ X.
(2) F is equicontinuous at each x ∈ X, that is for all γ > 0 there exists The first condition is guaranteed forṽ ε due to (28), while the equicontinuity property follows from the bounds (29): Applying Arzelà-Ascoli theorem completes the proof.

Auxiliary minimization problem
We assume that the domains Y i , Y m , Y e are given in cylindric coordinates (x 1 , r, φ) by (x 1 , r) ∈ Y i , (x 1 , r) ∈ Y m , (x 1 , r) ∈ Y e . Y m is a simply connected domain whose boundary is naturally divided into two parts Γ mi = ∂Y m ∩ ∂Y i and Γ me = ∂Y m ∩ ∂Y e . The first part is the segment {r 0 } × (a, b), while the second one is a smooth curve which never intersects or touches Y i except at endpoints A = (a, r 0 ) and B = (b, r 0 ), and locally near these points it is given by r = r a (x 1 ) Figure 3. Cross-section of the periodicity cell in the neighborhood of a Ranvier node. and r = r b (x 1 ). Moreover, we assume that r a and r b are C 2 -functions whose derivatives do not vanish at points a and b.
Let σ δ be given by

Consider the minimization problem
where the infimum is taken over 1-periodic in x 1 -variable functions θ ∈ H 1 (Y \ Γ), [θ] denotes the jump of θ across Γ, [θ] = θ i − θ e , θ i and θ e being limit values (traces) of θ on Γ from Y i and Y e , correspondingly. It is easy to see that the infimum in (31) is attained on a function θ δ which is defined up to a multiplicative and an additive constant, and θ δ satisfies Moreover, thanks to the radial symmetry θ δ = θ δ (x 1 , r) and Lemma 4.1.
(i) The infimum in (31) admits the bound then and the following uniform in δ > 0 bound holds: Proof. (i) We begin by constructing an approximation of θ δ away from points A and B. There exists a function Θ ∈ C 2 loc (Y \ Γ ) such that 0 ≤ Θ ≤ 1, and Θ = 1 in where ∇ 2 Θ denotes norm of the Hessian of Θ. Since |∇Θ| blows up at points A and B with the rate 1/dist(x , {A}) and 1/dist(x , {B}), any such a function Θ does not belong to H 1 (Y \ Γ), hence it is to be corrected near endpoints A and B of Γ . For simplicity we assume that in a neighborhood of points A and B the boundary of domain Y m is formed by two rays with angles ϕ A and ϕ B . Consider the δ-neighborhood D δ (B) of the point B and pass to polar coordinates (ρ, ϕ) with the center at B. Note that for sufficiently small δ the set Y m ∩ D δ (B) is a circular sector given by 0 < ϕ < ϕ B and 0 < ρ < δ. We set and α δ solving the transcendental equation There is a unique solution α δ of (42) on (0, 1/(2δ)) and it is asymptotically given by Note that θ B δ is continuous on R 2 \ R − and accordingly. To glue Θ, θ A δ and θ B δ together introduce a cut-off function χ ∈ C ∞ (R) such that χ(ρ) = 0 for ρ ≥ 1 and χ(ρ) = 1 for ρ ≤ 1/2. Set The last two integrals are similar and we consider only the second one: Also, where Λ is given by (ii) Convergences in (38) easily follow from the bound (36) and the Poincare inequality. To prove (39) multiply the equation in (34) by θ δ |θ δ | p−2 , p ≥ 2, and integrate over Y \ Γ to find, after integrating by parts, Therefore we have where Γ i and Γ e denote opposite sides of the surface Γ. Thus for p ≥ 2 it holds with C independent of δ and p ≥ 2. This yields H 1 -bounds for |θ δ | p/2 in Y i and Y e , that in turn lead to bounds for traces of |θ δ | p/2 on Γ i and Γ e : Since H 1/2 (Γ) is continuously embedded in L 2q (Γ) for some q > 1 (optimal q = 2) we have (45) It follows from (36) and (37) that Then iterative use of (45) yields for every integer k ≥ 0. The series ∞ 1 log(C 2 q j )/q j converges, hence θ δ L ∞ (Γi∪Γe) ≤ C. Finally by the maximum principle θ δ satisfies the same L ∞ -bound on Y \ Γ.
Next we show that the bound (36) for λ δ is in fact precise to the leading order.
Lemma 4.2. The following asymptotic result holds: where Λ is given by (44).
The following pointwise bounds hold in these domains: Thus, Observe that θ δ on S i,1 and S e,1 is sufficiently close to its mean values over Y i and Y e , correspondingly. Namely, by Hardy's inequality This leads to the following It remains to calculate the integral in the right hand side integrating by parts This completes the proof of the Lemma.
Next we show that, θ δ being normalized by (37), one has To this end we use the test functionθ δ constructed in the proof of Lemma 4.1 to write where we have used Lemma 4.2 together with the fact that θ δ minimizes (31), and calculations from the proof of Lemma 4.1. Representingθ δ asθ δ = (θ δ − θ δ ) + θ δ and expanding the left hand side of (58) we get where we have used (48) to derive the last inequality. Straightforward calculations show that 1 √ δ σ δ ∇θ δ 0, which in conjunction with (60) yields (57). We summarise properties of θ δ in Integrating by parts the second integral in (63) containing ∇ y θ ε 2 x ε and using (34) and Lemma 4.2, we have To pass to the two-scale limit in (64) we use (iv) in Lemma 3.5 and (i) in Lemma 4.4 and get Ye σ e (e 1 ∂ x1 u e 0 + ∇ y w e ) · (e 1 ∂ x1 U e (t, x 1 ) + ∇ y U 1 (t, x 1 , y)) dy dx 1 dt.