Hysteresis and Linear Stability Analysis on Multiple Steady-State Solutions to the Poisson--Nernst--Planck equations with Steric Interactions
aa r X i v : . [ phy s i c s . c o m p - ph ] S e p Hysteresis and Linear Stability Analysis on Multiple Steady-StateSolutions to the Poisson–Nernst–Planck equations with StericInteractions
Jie Ding ∗ Hui Sun † Shenggao Zhou ‡ September 18, 2020
Abstract
In this work, we numerically study linear stability of multiple steady-state solutions to atype of steric Poisson–Nernst–Planck (PNP) equations with Dirichlet boundary conditions,which are applicable to ion channels. With numerically found multiple steady-state solutions,we obtain S -shaped current-voltage and current-concentration curves, showing hysteretic re-sponse of ion conductance to voltages and boundary concentrations with memory effects.Boundary value problems are proposed to locate bifurcation points and predict the local bi-furcation diagram near bifurcation points on the S -shaped curves. Numerical approaches forlinear stability analysis are developed to understand the stability of the steady-state solutionsthat are only numerically available. Finite difference schemes are proposed to solve a derivedeigenvalue problem involving differential operators. The linear stability analysis reveals thatthe S -shaped curves have two linearly stable branches of different conductance levels and onelinearly unstable intermediate branch, exhibiting classical bistable hysteresis. As predicted inthe linear stability analysis, transition dynamics, from a steady-state solution on the unstablebranch to a one on the stable branches, are led by perturbations associated to the mode of thedominant eigenvalue. Further numerical tests demonstrate that the finite difference schemesproposed in the linear stability analysis are second-order accurate. Numerical approachesdeveloped in this work can be applied to study linear stability of a class of time-dependentproblems around their steady-state solutions that are computed numerically. Keywords:
Steric Poisson–Nernst–Planck Equations; Multiple Solutions; Linear StabilityAnalysis; Bistability; Current-Voltage Curve; Bifurcation Diagram
Ion channels play crucial roles in regulating various biological functions, such as exchanging ionsacross cell membranes and maintaining membrane excitability [1]. Ion permeation through chan-nels is sensitively related to the applied transmembrane voltages, ionic concentrations of cells ∗ Department of Mathematics and Mathematical Center for Interdiscipline Research, Soochow University, 1 ShiziStreet, Suzhou 215006, Jiangsu, China. Email: [email protected]. † Department of Mathematics and Statistics, California State University, Long Beach, CA, U. S. A. Email:[email protected]. ‡ Corresponding author. Department of Mathematics and Mathematical Center for Interdiscipline Research,Soochow University, 1 Shizi Street, Suzhou 215006, Jiangsu, China. Email: [email protected]. I - V ) characteristic curve of ion channels [46].It is desirable to perform stability analysis of the steady-state solutions in the sense of time evo-2ution. Linear stability analysis is a powerful tool to understand the stability of solutions underperturbations [40, 47–51]. However, most of studies on linear stability analysis are carried outaround steady states that have closed-form solutions. We develop general numerical approachesto perform linear stability analysis on solutions that are only numerically available, and applythem to numerically study the stability of steady-state solutions of the steric PNP equations inthe context of ion channels.In this work, we study linear stability of multiple steady-state solutions to a type of stericPNP equations that has been used to predict ion permeation through ion channels [35, 39]. Thecurrent-voltage ( I - V ) characteristic curve of ion channel is determined by finding the steady-statesolutions of the steric PNP equations with Dirichlet boundary conditions for ionic concentrationsand the electrostatic potential. The existence of multiple steady-state solutions gives rise to an S -shaped I - V curve with a hysteresis loop, which demonstrates the dependence of ion conductanceon varying applied voltages with memory effects. In addition, we find hysteretic response of ionconductance to boundary concentrations, indicating that ion conductance may not only dependon intracellular and extracellular ionic concentrations but also their history values.We propose boundary value problems to predict the location of bifurcation points and corre-sponding local bifurcation diagram on the S -shaped I - V curve. Numerical approaches for linearstability analysis are developed to understand the stability of the steady-state solutions that areonly numerically available. Finite difference schemes are proposed to discretize a derived eigen-value problem involving differential operators. By solving the eigenvalue problem, the linearstability analysis reveals that both the lower and upper branches of the S -shaped I - V curve arelinearly stable and the intermediate branch is linearly unstable, exhibiting a classical bistable hys-teresis diagram. Furthermore, we numerically study the transition dynamics in nonlinear regime,from a steady-state solution on the unstable branch to a one on the stable branches. As predictedby the linear stability analysis, such transition dynamics is led by perturbations associated tothe mode of the dominant eigenvalue. Numerical simulations also demonstrate that the localbifurcation diagram predicted by the second derivative of V with respect to I agrees with the I - V curve. In addition, we perform numerical tests to show the convergence of the finite differenceschemes proposed in the linear stability analysis.The rest of the paper is organized as follows. In section 2, we describe our governing equations.In section 3, we present our methods to solve for steady-state solutions and bifurcation points. Insection 4, we perform linear stability analysis on steady-state solutions. Section 5 is devoted tothe presentation of our numerical results. Finally, we draw conclusions in section 6. Consider a charged system occupying a bounded domain Ω with a smooth boundary ∂ Ω. Dirichletboundary conditions are considered on ∂ Ω with boundary data given by V . The electrostaticpotential ψ satisfies the following boundary-value problem of the Poisson equation: −∇ · ε ε r ∇ ψ = M X l =1 z l ec l + ρ f in Ω ,ψ = V on ∂ Ω , c l is the ion concentration for the l -th species, z l is the valence, M is the number of ionspecies, e is the elementary charge, ε is the vacuum permittivity, ε r is the relative dielectricconstant, and ρ f is the fixed charge density. The free-energy functional of the charged system isgiven by [35, 39, 52]: F [ c , c , · · · , c M ] = 12 Z Ω M X l =1 z l ec l + ρ f ! ψd x + β − Z Ω M X l =1 c l (cid:2) log (cid:0) Λ c l (cid:1) − (cid:3) d x + M X l =1 M X r =1 β − Z Z Ω c l ω lr c r d x − Z ∂ Ω ε ε r ∂ψ∂ n V dS, (2.1)where β is the inverse thermal energy, Λ is the thermal de Broglie wavelength, and ω ij are thesteric interaction coefficients, depending on the size of i th and j th ionic species. The coefficientscan be obtained by introducing spatially band-limited functions to truncate the (spatial) frequencyrange of the singular Lennard-Jones interactions between ions [35, 39], and are given by ω ij = ε ij ( a i + a j ) S δ , where ε ij is the energy coupling constant between the i th and j th ionic species, a i and a j are thecorresponding radii of ions, and S δ = ω d − d δ − d with ω d being the surface area of d dimensionalunit ball and δ being a small number for cutoff length. Note that the steric interaction coefficients ω ij can be alternatively defined via the second-order virial coefficients for hard spheres [41].The evolution of ion concentrations is described by the continuity equation ∂ t c l = −∇ · J l , where the flux J l is given by J l = − βD l c l ∇ (cid:18) δFδc l (cid:19) , l = 1 , , · · · , M. Here D l is the diffusion constant of l th ionic species. This gives rise to the Nernst–Planck equations ∂ t c l = D l ∇ · ∇ c l + z l eβc l ∇ ψ + c l M X r =1 ω lr ∇ c r ! , l = 1 , , · · · , M. We introduce a macroscopic length scale L , a characteristic concentration c , a characteristicscreening length λ D = q ε ε r βe c , and the following dimensionless variables [15, 25]˜ x = x L , ˜ D l = D l D , ˜ t = tD Lλ D , ˜ c l = c l c , ˜ ρ f = ρ f ec , ˜ ω lr = c ω lr , ˜ ψ = βeψ. After rescaling, we obtain by dropping all the tildes the nondimensionalized Poisson–Nernst–Planck equations with steric interactions ∂ t c l = γ l ∇ · ∇ c l + z l c l ∇ ψ + c l M X r =1 ω lr ∇ c r ! , l = 1 , , · · · , M, − ∆ ψ = κ M X l =1 z l c l + ρ f ! , (2.2)4here γ l = λ D L D l and κ = L λ . We apply the steric PNP equations (2.2) to study ion conductance through a perfect cylindershaped ion channel. By assuming radial and azimuthal symmetry, we reduce (2.2) into one-dimensional equations on a rescaled domain Ω = ( − , ∂ t c l = γ l ∂ x ∂ x c l + z l c l ∂ x ψ + c l M X r =1 w lr ∂ x c r ! , l = 1 , . . . , M, − ∂ xx ψ = κ M X l =1 z l c l + ρ f ! . (2.3)Note that the dimension reduction used here is a special case of a more general one-dimensionalreduction with a cross-sectional area factor, which allows more general geometry [21, 22].We are interested in the case that two ends of an ion channel are connected to extracellularand intracellular ionic solution reservoirs. To model an applied transmembrane potential andextracellular and intracellular bulk concentrations in reservoirs, we employ the following boundaryconditions [20, 24, 39, 53] c l ( t, −
1) = c Ll and c l ( t,
1) = c Rl , l = 1 , . . . , M,ψ ( −
1) = 0 and ψ (1) = V, (2.4)where V is an applied potential difference across the channel and the bulk concentrations c Ll and c Rl at boundaries satisfy neutrality conditions M X l =1 z l c Ll = M X l =1 z l c Rl = 0 . To describe fixed charges along a channel, we consider the following fixed charge distribution ρ f ( x ) = ρ i for x ∈ ( x i − , x i ) , (2.5)where i = 1 , , · · · , K , x = − x K = 1, and ρ = ρ K = 0. To quantify ion conductance througha channel, we define the net ionic current I = M X l =1 z l J l . (2.6)We remark that ions often have strong correlations due to crowded confinement in ion channels.Although the steric interactions can be regarded as a part of ion correlations of short range, thetreatment on ion correlations in this work is incomplete. As shown in [39], however, selective ionpermeation through channels can be predicted to some extent by the inclusion of the pairwiseinteraction term P Ml =1 P Mr =1 β − R R Ω c l ω lr c r d x , with some parameters in ω ij treated as the fittingparameters. We leave for future work a detailed study with inclusion of ion correlations [54–56].5 Steady States
We find steady-state solutions of the equations (2.3) by solving equations ∂ x ∂ x c l + z l c l ∂ x ψ + c l M X r =1 w lr ∂ x c r ! = 0 , l = 1 , . . . , M, − ∂ xx ψ = κ M X l =1 z l c l + ρ f ! (3.1)with boundary conditions (2.4). We briefly recall the numerical methods proposed in our previouswork [46] for solving steady-state solutions. We first convert the problem into a system of 2 M + 2first order ODEs with the same number of boundary conditions. The ODE system is numericallysolved with a Matlab program called BVP4C, in which an ODE system is solved with a collocationmethod on an adaptive mesh and the resulting nonlinear algebraic equations are iteratively solvedby Newton-type methods [57]. Given the boundary conditions (2.4), the problem with sign-alternating piecewise constant fixed charges (2.5) may have multiple steady-state solutions. Asthe applied voltage V varies, the corresponding current I may become multi-valued, giving rise toan S -shaped current-voltage ( I - V ) characteristic curve; cf. Figure 1. To compute an I - V curve,it is helpful to employ the strategy of numerical continuation on applied voltages to get a goodinitial guess. However, continuation on V often only finds the low-current branch and high-currentbranch, and misses the intermediate branch of an S -shaped I - V curve when multiple solutionsare present; cf. Figure 1. Moreover, the continuation advances with a very small stepsize as V approaches bifurcation points on the curve, where the Jacobian of the discretized nonlinear systembecomes close to singular. To address this issue, the voltage V , instead, is viewed as a function of I , and the equations (3.1) are solved with the boundary condition ψ (1) = V in (2.4) replaced by aprescribed current I . A complete S -shaped I - V curve can be obtained by performing continuationon I and collecting the applied voltage by V = ψ (1). See the work [46] for more details. The current-voltage ( I - V ) characteristic curve plays a crucial role in the study of ion conductancethrough ion channel or nanopores. It is of practical significance to locate bifurcation points onan S -shaped I - V curve, since the corresponding V values of these critical points are thresholdvalues for hysteresis to take place, and they are also endpoints of intervals for V in which multiplesolutions exist. If we assume that an I - V curve locally determines a differentiable function of V on I , then ∂ I V = 0 holds at bifurcation points of an S -shaped I - V curve; cf. red dots in Figure 1.More generally, it is of interest to locate a point on an I - V curve such that the slope, ∂ I V , at thepoint is given by a known value.At steady states, the current I given in (2.6) is a constant over Ω. Taking derivatives of both6ides of the equations (3.1) with respect to I , we have ∂ x ∂ x ˆ c l + z l ˆ c l ∂ x ψ + z l c l ∂ x ˆ ψ + ˆ c l M X r =1 w lr ∂ x c r + c l M X r =1 w lr ∂ x ˆ c r ! = 0 , l = 1 , · · · , M, − ∂ xx ˆ ψ = κ M X l =1 z l ˆ c l , (3.2)where ˆ ψ := ∂ I ψ and ˆ c l := ∂ I c l for l = 1 , . . . , M . To locate a point determined by a given slopevalue ∂ I V , we solve the following boundary-value problem (BVP) Equations (3.1) , Equations (3.2) , ( c , · · · , c M , ψ, ˆ c , · · · , ˆ c M , ˆ ψ ) (cid:12)(cid:12)(cid:12) x = − = ( c L , · · · , c LM , , , · · · , , , ( c , · · · , c M , ˆ c , · · · , ˆ c M , ˆ ψ ) (cid:12)(cid:12)(cid:12) x =1 = ( c R , · · · , c RM , , · · · , , ∂ I V ) . (3.3)Note that we do not specify the applied voltage V as a boundary condition in the BVP (3.3).Instead, the applied voltage can be found after solving the BVP. Since P Ml =1 z l ∂ I J l = 1, we canconvert the BVP into a system of 4 M + 3 first order ODEs with the same number of boundaryconditions, and numerically solve the ODEs again with the BVP4C program. See our previouswork [46] for more details. Of particular interest are the bifurcation points that are determinedby ∂ I V = 0 on an S -shaped I − V curve. With found bifurcation points, we are able to determinethreshold applied voltages for hysteresis, as well as intervals for V in which the steric PNPequations have multiple steady states.Consider an I - V curve that locally determines a twice differentiable function of V on I . It isdesirable to compute the second derivative ∂ II V that reveals local shape (concavity) of the curve.In particular, the second derivative enables prediction of approximate local bifurcation diagramclose to bifurcation points. Taking second order derivatives of both sides of the equations (3.1)with respect to I , we have ∂ x " ∂ x ˆˆ c l + 2 z l ˆ c l ∂ x ˆ ψ + z l ˆˆ c l ∂ x ψ + z l c l ∂ x ˆˆ ψ + M X r =1 ω lr (cid:16) c l ∂ x ˆ c r + c l ∂ x ˆˆ c r + ˆˆ c l ∂ x c r (cid:17) = 0 , − ∂ xx ˆˆ ψ = κ M X l =1 z l ˆˆ c l , (3.4)where ˆˆ ψ := ∂ II ψ and ˆˆ c l := ∂ II c l for l = 1 , . . . , M . Here c l , ˆ c l , ψ , and ˆ ψ in (3.4) are obtained bynumerically solving the BVP (3.3). The boundary conditions are given by(ˆˆ c , · · · , ˆˆ c M , ˆˆ ψ ) (cid:12)(cid:12)(cid:12) x = − = (0 , · · · , ,
0) and (ˆˆ c , · · · , ˆˆ c M ) (cid:12)(cid:12)(cid:12) x =1 = (0 , · · · , . (3.5)By the fact that M X l =1 z l ∂ II J l = 0 ,
7e can convert the equations (3.4) into a system of 2 M +1 first order ODEs with the same numberof boundary conditions, and numerically solve the ODEs with the BVP4C program. After solving,we have ∂ II V = ˆˆ ψ (1). When a bifurcation point is interested, we first solve the BVP (3.3) with ∂ I V = 0 to get c l , ˆ c l , ψ , and ˆ ψ , and then compute ∂ II V to characterize the local bifurcationdiagram close to the bifurcation point. We present a linear stability analysis on steady-state solutions of the steric PNP equations (2.3).We linearize the nonlinear equations (2.3) around a steady-state solution denoted by ¯ ψ ( x ) and¯ c l ( x ), which is obtained by solving the equations (3.1) as described in Section 3.1. We assumethat the solution ( c , · · · , c M , ψ ) to the equations (2.3) is given by c ( x, t )... c M ( x, t ) ψ ( x, t ) = ¯ c ( x )...¯ c M ( x )¯ ψ ( x ) + e λt δc ( x )... δc M ( x ) δψ ( x ) + · · · , where | δc l | ≪ | δψ | ≪
1, and dots represent terms of higher orders as | δc l | and | δψ | go tozero. We also assume that the initial condition is given by c l ( x,
0) = ¯ c l ( x ) + δc l ( x ) and ψ ( x,
0) = ¯ ψ ( x ) + δψ ( x ) . By the Dirichlet boundary conditions (2.4), we have homogeneous boundary conditions for δc l and δψ : ( δc , δc , · · · , δc M , δψ ) | x = ± = (0 , , · · · , , . (4.1)To the leading order, we recover the equations (3.1) for the steady state. Comparing the first-orderterms, we obtain a linearized system λδc l = γ l ∂ x " ∂ x δc l + z l ¯ c l ∂ x δψ + z l δc l ∂ x ¯ ψ + M X r =1 w lr ( δc l ∂ x ¯ c r + ¯ c l ∂ x δc r ) , ≤ l ≤ M, − ∂ xx δψ = κ M X l =1 z l δc l . This system can be expressed in a matrix form A δc δc ... δc M = λ δc δc ... δc M , (4.2)8here A = A A · · · A M A A · · · A M ... ... ... ... A M A M · · · A MM (4.3)with operators A ii γ i = ∂ xx − z i κ∂ x ¯ c i ∂ x ( ∂ xx ) − − z i κ ¯ c i + z i ∂ x ¯ ψ∂ x + z i ∂ xx ¯ ψ + w ii ∂ x ¯ c i ∂ x + w ii ¯ c i ∂ xx + M X r =1 w ir ( ∂ x ¯ c r ∂ x + ∂ xx ¯ c r ) for i = 1 , , · · · , M,A ij γ i = − z i z j κ∂ x ¯ c i ∂ x ( ∂ xx ) − − z i z j κ ¯ c i + w ij ∂ x ¯ c i ∂ x + w ij ¯ c i ∂ xx for i = j. To numerically compute λ , we develop a finite difference scheme to discretize the operatormatrix (4.3). We introduce a uniform mesh with grid points given by x i = − ih, i = 0 , , · · · , N + 1 , where grid spacing h = N +1 . We denote by δc l,i and δψ i the numerical approximation of δc l ( x i )and δψ ( x i ) for i = 0 , , · · · , N + 1, respectively. By the boundary conditions (4.1), we have δc l, = δc l,N +1 = 0 and δψ = δψ N +1 = 0. We take interior grid points { x i } Ni =1 as our computationalmesh, and define c δc l = ( δc l, , δc l, , · · · , δc l,N ) T , ≤ l ≤ M, c δψ = ( δψ , δψ , · · · , δψ N ) T . (4.4)The first derivative is approximated by a second-order central differencing scheme, which corre-sponds to a differentiation matrix D N +2 = 12 h − − − ( N +2) × ( N +2) . The second derivative corresponds to a differentiation matrix D N +2 . To approximate the eigen-value problem (4.2) on interior grid points, we only make use of ˜ D N and ˜ D N obtained by strippingthe first and last rows and columns of D N +2 and D N +2 , i.e., in Matlab notations,˜ D N = D N +2 (2 : N + 1 , N + 1) and ˜ D N = D N +2 (2 : N + 1 , N + 1) . The operators ∂ x and ∂ xx with homogenous Dirichlet boundary conditions can be approximatedby ˜ D N and ˜ D N , respectively. Similarly, the operator ( ∂ xx ) − is approximated by ( ˜ D N ) − . Con-9equently, the A ij blocks can be approximated as A ii γ i ≈ ˜ D N − z i κ ˜ D N ¯ c i ˜ D N ( ˜ D N ) − − z i κ ¯ c i + z i ˜ D N ¯ ψ ˜ D N + z i ˜ D N ¯ ψ + w ii ˜ D N ¯ c i ˜ D N + w ii ¯ c i ˜ D N + M X r =1 w ir ( ˜ D N ¯ c r ˜ D N + ˜ D N ¯ c r ) , ≤ i ≤ M, A ij γ i ≈ − z i z j κ ˜ D N ¯ c i ˜ D N ( ˜ D N ) − − z i z j κ ¯ c i + w ij ˜ D N ¯ c i ˜ D N + w ij ¯ c i ˜ D N , i = j. We assemble the vectors c δc l defined in (4.4) for M species of ions, and define δ c := M X l =1 e l ⊗ c δc l , where e l is a unit column vector whose the l -th element is the only nonzero element being one, and ⊗ represents a Kronecker product. We recast the matrix A in (4.2) as a matrix A : R MN → R MN given by A = M X i,j =1 (cid:0) e i ⊗ e Tj (cid:1) ⊗ A ij . Therefore, the eigenvalue problem (4.2) can be approximated by A δ c = λδ c . (4.5)Such an eigenvalue problem determines the linear stability of the steric PNP equations (2.3)around a steady state as follows: Re( λ ) > λ ) < λ ) > c n +1 l − c nl ∆ t = γ l ∂ x ∂ x c n +1 l + z l c n +1 l ∂ x ψ n +1 + c n +1 l M X r =1 ω lr ∂ x c n +1 r ! , − ∂ xx ψ n +1 = κ M X l =1 z l c n +1 l + ρ f ! , (4.6)with boundary conditions (2.4). Here we denote by c nl ( x ) and ψ n ( x ) semi-discrete approximationsof c l ( t n , x ) and ψ ( t n , x ), where t n = n ∆ t and ∆ t is a uniform time step size. Such a semi-discretescheme can be converted into an ODE system and solved with the BVP4C, as described inSection 3.1. 10 Numerical Results
In our numerical simulations, we numerically solve the equations (3.1) with certain boundaryconditions to find an I - V curve. To understand more on an S -shaped I - V curve, we solve theBVP (3.3) with ∂ I V = 0 to locate bifurcation points, and solve the equations (3.4) at bifurca-tions points with boundary conditions (3.5) to find ∂ II V , which can be further used to predictlocal bifurcation diagram close to the bifurcations points. Moreover, we solve the eigenvalueproblem (4.5) to study linear stability of steady states. In nonlinear regime, we numerically in-vestigate transition dynamics from unstable steady states to stable ones with the time integrationscheme (4.6).Unless otherwise stated, we consider a charged system that consists of binary monovalentelectrolytes, i.e., M = 2, z = 1, and z = −
1, and fixed charges ρ f ( x ) = , x ∈ [ − , − . , − , x ∈ [ − . , . , , x ∈ [0 . , . , − , x ∈ [0 . , . , , x ∈ [0 . , . We consider a steric interaction matrix W with diagonal elements w = w = 0 . w = w = 0 .
01. We take γ = γ = 1, κ =100, and boundary concentrations c L = 1, c L = 1, c R = 0 .
5, and c R = 0 . We study the current-voltage relation for voltage-induced hysteresis of ionic conductance. i.e., the I - V curve as shown in Figure 1. We first numerically solve the equations (3.1) with continuationon applied voltages. When the applied voltage increases gradually, the current ascends along alow-current branch (from D to C on the solid red curve). As the voltage approaches the criticalvalue corresponding to F , the iterations for solving the equations slow down due to a more andmore singular Jacobian matrix, and continuation has to advance with small increment. Afterthe applied voltage exceeds a critical value about V = 154 .
86, the current suddenly jumps fromthe low-current branch at the bifurcation point F to a high-current branch at H , switching theionic conductance state. In the other round of sweeping with decreasing applied voltages fromhigh values, the current follows a totally different path (from I to A on the solid blue curve)and suddenly switches at the bifurcation point G with V = 148 .
28 to E on the low-currentbranch, exhibiting a typical hysteresis loop. This demonstrates that ion conductance dependson the applied voltage as well as its history values, showing a memory effect. We remark thatthis is reminiscent of the voltage-induced gating mechanism [43], though our theory is purelydeterministic. To explore more of the intermediate region of two branches, we view the appliedvoltage as a function of the current and solve the equations (3.1) with a boundary conditionprescribing the current, instead of applied voltages. Such treatment of the boundary conditionstablizes the numerical simulation of steady-state solutions that are shown to be unstable. Asshown by the green dashed curve in Figure 1, an intermediate branch connecting the low-currentand high-current branches can be found by continuation on the current. The complete S -shapedcurve features hysteretic response of ion current to varying applied voltages.11
30 140 150 160 170 V I BACD EG H IF G F Figure 1: An S -shaped I - V curve obtained by solving the equations (3.1). Linear stability analysisreveals that the steady-state solutions on the upper (blue solid) and lower (red solid) branches arestable, and those on the intermediate (green dashed) branches are unstable. The second derivativeis given by ∂ II V = 0 .
005 at the upper bifurcation point ( G ) that is located at (148 . , . ∂ II V = − .
011 at the lower bifurcation point ( F ) that is located at(154 . , . A , B , and C , when the appliedvoltage V = 150. The zoomed-in inset presents local bifurcation diagram, in which the blackdotted parabolas are predicted by the second derivatives at bifurcation points. A c c B c c C c c -1 0 1 x
0 50 100150 -1 0 1 x -1 0 1 x Figure 2: The solution profiles of the electrostatic potential and ionic concentrations, correspond-ing to A , B , and C in the Figure 1 when V = 150.12or applied voltages in the interval (148 . , . A , B , and C , are found for V = 150 . The panel presented in Figure 2 displays profiles of the electrostaticpotential and ionic concentrations corresponding to these three points. We observe that thesolutions have large variations at the locations where the fixed charges are discontinous. The ionsdistribute mainly according to the fixed charges due to electrostatic interactions. Comparing withthe low-current state, both cations and anions have larger concentrations along the channel in thehigh-current conductance state. −300 −200 −100 0−0.500.51 I m ( λ ) Re( λ )A −1 −0.8 −0.6−101 −200 −100 0 Re( λ )B Re( λ )C −0.1 −0.05 0−101 Figure 3: The eigenvalues with the largest real parts of the eigenvalue problem (4.5) at steady-state solutions corresponding to the A , B , and C in the Figure 1. Insets are zoomed-in plots ofthe eigenvalue with the largest real part.We perform linear stability analysis on the steady-state solutions on the lower, intermediate,and upper branches of the S -shaped I - V curve by solving the eigenvalue problem (4.5). Thesign of the computed eigenvalue with the largest real part, i.e., the dominant eigenvalue, for theperturbed system determines the linear stability. We present in Figure 3 six eigenvalues with thelargest real parts of the eigenvalue problem (4.5) at steady-state solutions, corresponding to the A , B , and C in the Figure 1. We find that the real part of the dominant eigenvalue is negative forboth A and C , and is positive for B . This reveals that the steady-state solutions at A and C arelinearly stable, but that at B is linearly unstable. Analogously, we consider the real part of thedominant eigenvalue for other points on the I - V curve. Remarkably, we find that the steady-statesolutions on both the lower and upper branches shown in blue solid curves are stable, but those onthe intermediate branch shown in red dashed curves are unstable, exhibiting a classical bistablediagram. This agrees with our expectation of bistability as seen in typical hysteresis loops. Also,the instability explains why it is hard to capture the steady-state solutions on the intermediatebranch in our numerical computations [46].The stable and unstable branches connect at bifurcation points marked by red solid dots, whichare located by solving the BVP (3.3) with ∂ I V = 0. Further numerical calculations based on theequations (3.4) at bifurcations points with boundary conditions (3.5) find the second derivative of13 with respect to I : ∂ II V = − .
011 at the lower bifurcation point and ∂ II V = 0 .
005 at the upperbifurcation point. Such second derivative information helps predict local bifurcation diagram closeto the bifurcation points. The dotted parabolas shown in the inset of Figure 1 demontrate thatthe predicted parabola agrees well locally with the one calculated by using steady-state solutions. h A B CError in λ A Max
Order Error in λ B Max
Order Error in λ C Max
Order × − - 1.44 × − × − - × − × − × − × − × − × − × − × − × − λ Max , which is the eigenvalue with the largest realpart, for the eigenvalue problem at A , B , and C in the Figure 1.We perform convergence tests with various mesh resolutions to further validate the proposednumerical approaches on linear stability analysis of the steric PNP equations around steady-statesolutions. Table 1 lists the error and convergence order of λ Max , which is the eigenvalue withthe largest real part, for the eigenvalue problem at A , B , and C in the Figure 1. Table 2 lists l ∞ errors and convergence order of the eigenvector, corresponding to the eigenvalue λ Max , forthe ionic concentrations and electrostatic potential at A , B, and C . The error is obtained bycomparing against a reference solution that is computed on a highly refined mesh with h = .Since second-order central differencing is employed to discretize the differential operators in (4.3),we anticipate a second-order convergence rate of our numerical approaches. From the tables,we can see that the numerical error decreases robustly as the mesh refines, and the numericalconvergence order is indeed roughly about two for both eigenvalues and eigenvectors. The linear stability analysis is no longer applicable, when long-time dynamics of the steric PNPequations are concerned. To understand how the unstable steady-state solution evolves to stableones, we numerically investigate the nonlinear transition dynamics with the time integrationscheme (4.6). The numerical simulations start from a steady-state solution corresponding to B onthe unstable branch of the S -shaped I - V curve shown in Figure 1. Since numerical noise providesbroadband perturbations, noise associated to the mode of the eigenvalue with Re( λ B Max ) > ψ and c l as ψ ∆ ( t, x ) := ψ ( t, x ) − ψ B ( x ) , c ∆ l ( t, x ) := c l ( t, x ) − c Bl ( x ) , l = 1 , · · · , M, where ψ ( t, x ) and c l ( t, x ) are numerical solutions to (4.6) with initial conditions ψ B ( x ) and c Bl ( x )that are a steady-state solution corresponding to B in Figure 1. As shown in Figure 4 (Left), both ψ ∆ ( t, x ) and c ∆ l ( t, x ) start from zero (blue lines marked with dots), and evolve to final states (redcurves marked with diamonds), corresponding to a stable solution at C in Figure 1. In addition,14 l ∞ -error in δc Order l ∞ -error in δc Order l ∞ -error in δψ Order A × − - 3.41 × − - 2.25 × − - × − × − × − × − × − × − × − × − × − B × − - 3.13 × − - 6.64 × − - × − × − × − × − × − × − × − × − × − C × − - 1.90 × − - 2.98 × − - × − × − × − × − × − × − × − × − × − l ∞ error and convergence order of the eigenvector corresponding to λ Max at the A , B , and C in the Figure 1. -100 5 c (a) -100 5 c -1 -0.5 0 0.5 1 x -100 5 -0.050 0.05 c (b) -0.050 0.05 c -1 -0.5 0 0.5 1 x -0.040 0.04 Figure 4: (a) Snapshots of the evolution of ψ ∆ ( t, x ) and c ∆ l ( t, x ) ( l = 1 , · · · , M ) starting from theunstable steady-state solution at B (blue lines marked with dots) to a final state correspondingto C (red curves marked with diamonds). (b): A normalized eigenvector corresponding to theeigenvalue with the largest real part at B . 15igure 4 (Right) shows a normalized eigenvector corresponding to the eigenvalue with the largestreal part at B . Clearly, we find that the evolution of concentrations and the electrostatic potentialfollows exactly in the direction dictated by the eigenvector, especially in the early stage. This isconsistent with our linear stability analysis that the perturbations associated to the mode of theeigenvalue with the largest real part lead transition dynamics. -100 5 c (a) -100 5 c -1 -0.5 0 0.5 1 x -100 5 -0.050 0.05 c (b) -0.050 0.05 c -1 -0.5 0 0.5 1 x -0.040 0.04 Figure 5: (a): Snapshots of the evolution of ψ ∆ ( t, x ) and c ∆ l ( t, x ) ( l = 1 , · · · , M ) starting from theinitial condition (5.1) (blue lines marked with dots) to a final state corresponding to A (red curvesmarked with diamonds). (b): Another normalized eigenvector corresponding to the eigenvaluewith the largest real part at B .To study transition dynamics from B to A , numerical integration starts from an initial con-dition ψ (0 , x ) = ǫψ A ( x ) + (1 − ǫ ) ψ B ( x ) , c l (0 , x ) = ǫc Al ( x ) + (1 − ǫ ) c Bl ( x ) , l = 1 , · · · , M, (5.1)where ǫ = 0 . B slightly toward that at A . Figure 5 (Left) shows the evolution of ψ ∆ ( t, x ) and c ∆ l ( t, x ) from theinitial condition (5.1) to a final state at A . Figure 5 (Right) shows another normalized eigenvectorcorresponding to the eigenvalue with the largest real part at B . Notice that this normalizedeigenvector points in an opposite direction contrast to the one showed in Figure 4 (Right). Again,we can observe that the evolution of concentrations and the electrostatic potential accords withthe given eigenvector, agreeing with our linear stability analysis. In this example, we investigate the effect of boundary concentrations on the existence and stabilityof multiple steady-state solutions to the steric PNP equations. We set boundary concentrations c l ( −
1) = c Ll for l = 1 , c B I Figure 6: A current-concentration curve showing concentration-induced hysteresis. The steady-states on the blue solid curve are linearly stable, and those on the red dashed curve are linearlyunstable.the relation c (1) = c (1) for boundary concentrations at the right end where the potential V = 150. We define the variable c B := c (1). Similar to the S -shaped I - V curve shown inFigure 1, two rounds of sweeping with an increasing and decreasing c B are used to capture low-current and high-current branches in the current-concentration curve shown in Figure 6. To findthe intermediate branch, we prescribe the current I as a boundary condition and treat c B as avariable to be determined. See our previous work [46] for more details. Such a strategy helpscomplete the whole current-concentration curve with a hysteresis loop, which indicates that thechannel conductance state not only depends on intracellular and extracellular ionic concentrationsbut also their history values.In addition, we study linear stability of steady state solutions on the current-concentrationcurve by solving the eigenvalue problem (4.5). For the points A , B, and C on the curve with c B = 1, we can see from Figure 7 that the dominant eigenvalue has a negative real part for A and C , and a positive real part for B . This indicates that the steady-state solutions at A and C arelinearly stable and that at B is linearly unstable. Similarly, we check the real part of the dominanteigenvalue for the whole current-concentration curve. Results demonstrate that the steady-statesolutions on blue solid curves are linearly stable, and that on the dashed intermediate branch inred are unstable, showing a typical bistable hysteresis loop. In this work, we have studied linear stability of multiple steady-state solutions to the steric PNPequations. The current-voltage ( I - V ) characteristic curve has been determined by solving steadystates to the steric PNP equations with Dirichlet boundary conditions arising from ion channels.We have numerically found multiple steady-state solutions that lead to an S -shaped current-voltage and current-concentration curves showing hysteresis loops. The hysteretic response of ion17
300 −200 −100 0−0.500.51 I m ( λ ) Re( λ )A −5.5 −5 −4.5−101 −300 −200 −100 0 Re( λ )B Re( λ )C −0.15 −0.1 −0.05−101 Figure 7: The eigenvalues with the largest real parts of the eigenvalue problem (4.5) at steady-state solutions corresponding to the points A , B , and C in the current-concentration curve shownin Figure 6. Insets are zoomed-in plots of the eigenvalue with the largest real part.conductance to applied voltages and boundary concentrations demonstrates that ion conductancestate depends on their currently applied values as well as their history values, i.e, the memoryeffect. We have proposed boundary value problems to locate of bifurcation points and to predictlocal bifurcation diagram near bifurcation points on the S -shaped I - V curve. To understand thestability of steady-state solutions that are only numerically available, we have performed linearstability analysis that is based on finite difference discretization of differential operators. The lin-ear stability analysis has unveiled that the S -shaped I - V curve has two linearly stable brancheswith a high-conductance state and low-conductance state, and has one linearly unstable interme-diate branch. This presents a classical hysteresis diagram with bistable behavior. Furthermore,extensive numerical tests have shown that the numerical method proposed in the linear stabilityanalysis is second-order accurate. In addition, we have numerically studied the nonlinear transi-tion dynamics from unstable steady-state solutions to stable solutions. It has been found that thetransition follows the direction predicted by eigenvectors associated to the dominant eigenvalue inthe linear stability analysis. The numerical approaches developed in this work are quite generaland useful, in the sense that they can be applied to study linear stability of other time-dependentproblems around their steady-state solutions that are only numerically available. Acknowledgments.
J. Ding was supported by National Natural Science Foundation of China(No. 21773165, 11771318, and 11790274) and Postgraduate Research & Practice Innovation Pro-gram of Jiangsu Province. H. Sun was supported by Simons Foundation Collaborative Grant(No. 522790). S. Zhou was supported by the grants NSFC 21773165, Natural Science Founda-tion of Jiangsu Province, Young Elite Scientist Sponsorship Program by Jiangsu Association forScience and Technology, and National Key R&D Program of China (No. 2018YFB0204404).18 eferences [1] J. Zheng and M. Trudeau.
Handbook of ion channels . CRC Press, 2015.[2] F. Sigworth. Voltage gating of ion channels.
Q. Rev. Biophys. , 27:1–40, 1994.[3] I. Josephson and Y. Cui. Voltage- and concentration-dependent effects of lidocaine on cardiacNa channel and Ca channel gating charge movements.
Pflugers Arch. , 428:485–491, 1994.[4] M. Pustovoit, A. Berezhkovskii, and S. Bezrukova. Analytical theory of hysteresis in ionchannels: Two-state model.
J. Chem. Phys. , 125, 2006.[5] B. Das, K. Banerjee, and G. Gangopadhyay. Entropy hysteresis and nonequilibrium thermo-dynamic efficiency of ion conduction in a voltage-gated potassium ion channel.
Phys. Rev.E , 86:061915, 2013.[6] C. Tilegenovaa, D. Cortesa, and L. Cuello. Hysteresis of KcsA potassium channel’s activation-deactivation gating is caused by structural changes at the channel’s selectivity filter.
Proc.Natl. Acad. Sci. USA. , 114:3234–3239, 2017.[7] T. Andersson. Exploring voltage-dependent ion channels in silico by hysteretic conductance.
Math. Biosci. , 226:16–27, 2010.[8] D. Fologea, E. Krueger, Y. Mazur, C. Stith, Y. Okuyama, R. Henry, and G. Salamo. Bi-stability, hysteresis, and memory of voltage-gated lysenin channels.
Biochim. Biophys. Acta ,1808:2933–2939, 2011.[9] V. Barcilon, D. Chen, and R. Eisenberg. Ion flow through narrow membrane channels: PartII.
SIAM J. Appl. Math. , 52:1405–1425, 1992.[10] A Singer and J Norbury. A Poisson-Nernst-Planck model for biological ion channels-anasymptotic analysis in a three-dimensional narrow funnel.
SIAM J. Appl. Math. , 70(3):949–968, 2009.[11] P. Markowich.
The Stationary Semiconductor Device Equations . Springer-Verlag, New York,1986.[12] M. Ward, L. Reyna, and F. Odeh. Multiple steady-state solutions in a multijunction semi-conductor device.
SIAM J. Appl. Math. , 51:90–123, 1991.[13] V. Barcilon, D. Chen, R. Eisenberg, and J. Jerome. Qualitative properties of steady-statePoisson–Nernst–Planck systems: Perturbation and simulation study.
SIAM J. Appl. Math. ,57:631–648, 1997.[14] M. Bazant, K. Chu, and B. Bayly. Current-voltage relations for electrochemical thin films.
SIAM J. Appl. Math. , 65(5):1463–1484, 2005.[15] M. Kilic, M. Bazant, and A. Ajdari. Steric effects in the dynamics of electrolytes at largeapplied voltages. II. Modified Poisson–Nernst–Planck equations.
Phys. Rev. E , 75:021503,2007. 1916] A. Singer, D. Gillespie, J. Norbury, and R. Eisenberg. Singular perturbation analysis of thesteady-state Poisson–Nernst–Planck system: Applications to ion channels.
Euro. J. Appl.Math. , 19:541–560, 2008.[17] X. Wang, D. He, J. Wylie, and H. Huang. Singular perturbation solutions of steady-statePoisson–Nernst–Planck systems.
Phys. Rev. E , 89:022722, 2014.[18] L. Ji, P. Liu, Z. Xu, and S. Zhou. Asymptotic analysis on dielectric boundary effects ofmodified Poisson–Nernst–Planck equations.
SIAM J. Appl. Math. , 78:1802–1822, 2018.[19] N. Abaid, R. Eisenberg, and W. Liu. Asymptotic expansions of I-V relations via a Poisson–Nernst–Planck system.
Phys. Rev. E , 7:1507–1526, 2008.[20] W. Liu. Geometric singular perturbation approach to steady-state Poisson–Nernst–Plancksystems.
SIAM J. Appl. Math. , 65:754–766, 2005.[21] B. Eisenberg and W. Liu. Poisson-Nernst-Planck systems for ion channels with permanentcharges.
SIAM J. Math. Anal. , 38:1932–1966, 2007.[22] W. Liu. One-dimensional steady-state Poisson–Nernst–Planck systems for ion channels withmultiple ion species.
J. Differ. Equ. , 246:428–451, 2009.[23] S. Ji, W. Liu, and M. Zhang. Effects of (small) permanent charges and channel geometry onionic flows via classical Poisson–Nernst–Planck models.
SIAM J. Appl. Math. , 75:114–135,2015.[24] B. Eisenberg, W. Liu, and H. Xu. Reversal permanent charge and reversal potential: casestudies via classical Poisson–Nernst–Planck models.
Nonlinearity , 28:103–127, 2015.[25] M. Bazant, K. Thornton, and A. Ajdari. Diffuse-charge dynamics in electrochemical systems.
Phys. Rev. E , 70(2):021506, 2004.[26] M. Bazant, M. S. Kilic, B. D. Storey, and A. Ajdari. Towards an understanding of induced-charge electrokinetics at large applied voltages in concentrated solutions.
Adv. Colloid. In-terface Sci. , 152:48–88, 2009.[27] M. Kilic, M. Bazant, and A. Ajdari. Steric effects in the dynamics of electrolytes at largeapplied voltages. I. Double-layer charging.
Phys. Rev. E , 75:021502, 2007.[28] B. Li. Minimization of electrostatic free energy and the Poisson–Boltzmann equation formolecular solvation with implicit solvent.
SIAM J. Math. Anal. , 40:2536–2566, 2009.[29] B. Li. Continuum electrostatics for ionic solutions with nonuniform ionic sizes.
Nonlinearity ,22:811–833, 2009.[30] M. Ma and Z. Xu. Self-consistent field model for strong electrostatic correlations and inho-mogeneous dielectric media.
J. Chem. Phys. , 141:244903, 2014.[31] S. Zhou, Z. Wang, and B. Li. Mean-field description of ionic size effects with non-uniformionic sizes: A numerical approach.
Phys. Rev. E , 84:021901, 2011.2032] B. Li, P. Liu, Z. Xu, and S. Zhou. Ionic size effects: generalized boltzmann distributions,counterion stratification, and modified debye length.
Nonlinearity , 26(10):2899, 2013.[33] B. Lu and Y. Zhou. Poisson-Nernst-Planck equations for simulating biomolecular diffusion-reaction processes II: Size effects on ionic distributions and diffusion-reaction rates.
Biophys.J. , 100:2475–2485, 2011.[34] Y. Qiao, X. Liu, M. Chen, and B. Lu. A local approximation of fundamental measure theoryincorporated into three dimensional Poisson–Nernst–Planck equations to account for hardsphere repulsion among ions.
J. Stat. Phys. , 163:156–174, 2016.[35] T. Lin and B. Eisenberg. A new approach to the Lennard-Jones potential and a new model:PNP-steric equations.
Commun. Math. Sci. , 12:149–173, 2014.[36] Y. Hyon, B. Eisenberg, and C. Liu. A mathematical model for the hard sphere repulsion inionic solutions.
Commun. Math. Sci. , 9:459–475, 2010.[37] S. Xu, P. Sheng, and C. Liu. An energetic variational approach for ion transport.
Commun.Math. Sci. , 12:779–789, 2014.[38] B. Eisenberg, Y. Hyon, and C. Liu. Energy variational analysis EnVarA of ions in waterand channels: Field theory for primitive models of complex ionic fluids.
J. Chem. Phys. ,133:104104, 2010.[39] T. Horng, T. Lin, C. Liu, and R. Eisenberg. PNP equations with steric effects: a model ofion flow through channels.
J. Phys. Chem. B , 116(37):11422–11441, 2012.[40] N. Gavish. Poisson–Nernst–Planck equations with steric effects - non-convexity and multiplestationary solutions.
Phys. D , 368:50–65, 2018.[41] J. Zhou, Y. Jiang, and M. Doi. Cross interaction derives stratification in drying film of binarycolloidal mixtures.
Phys. Rev. Lett. , 10:108002, 2017.[42] J. Ding, Z. Wang, and S. Zhou. Positivity preserving finite difference methods for Poisson–Nernst–Planck equations with steric interactions: Application to slit-shaped nanopore con-ductance.
J. Comput. Phys. , 397:108864, 2019.[43] N. Gavish, C. Liu, and B. Eisenberg. Do bistable steric Poisson–Nernst–Planck modelsdescribe single-channel gating?
J. Phys. Chem. B , 22(20):5183–5192, 2018.[44] T. Lin and B. Eisenberg. Multiple solutions of steady-state Poisson–Nernst–Planck equationswith steric effects.
Nonlinearity , 28:2053–2080, 2015.[45] L. Hung and M. Minh. Stationary solutions to the Poisson–Nernst–Planck equations withsteric effects. arXiv preprint , 1:02456, 2015.[46] J. Ding, H. Sun, Z. Wang, and S. Zhou. Computational study on hysteresis of ion channels:Multiple solutions to steady-state Poisson–Nernst–Planck equations.
Commun. Comput.Phys. , 23(5):1549–1572, 2018. 2147] M. S. Rosin and H. Sun. Stability and energetics of bursian diodes.
Phys. Rev. E , 87:043114,2013.[48] A. Yochelis. Spatial structure of electrical diffuse layers in highly concentrated electrolytes:A modified Poisson090808Nernst090808Planck approach.
J. Phys. Chem. C , 118:5716–5724,2014.[49] A. Yochelis. Transition from non-monotonic to monotonic electrical diffuse layers: impact ofconfinement on ionic liquids.
Phys. Chem. Chem. Phys. , 16:2836–2841, 2014.[50] B. Li, H. Sun, and S. Zhou. Stability of a cylindrical solute-solvent interface: Effect ofgeometry, electrostatics, and hydrodynamics.
SIAM J. Appl. Math. , 75(3):907–928, 2015.[51] N. Gavish and A. Yochelis. Theory of phase separation and polarization for pure ionic liquids.
J. Phys. Chem. Lett. , 7(7):1121–1126, 2016.[52] X. Liu, Y. Qiao, and B. Lu. Analysis of the mean field free energy functional of electrolytesolution with non-homogenous boundary conditions and the generalized PB/PNP equationswith inhomogeneous dielectric permittivity.
SIAM J. Appl. Math. , 78:1131–1154, 2018.[53] C. L. Gardner and J. R. Jones. Electrodiffusion model simulation of the potassium channel.
J. Theor. Bio. , 291:10–13, 2011.[54] D. Gillespie, W. Nonner, and R. S. Eisenberg. Density functional theory of charged, hard-sphere fluids.
Phys. Rev. E , 68:031503, 2003.[55] M. Bazant, B. D. Storey, and A. A. Kornyshev. Double layer in ionic liquids: Overscreeningversus crowding.
Phys. Rev. Lett. , 106:046102, 2011.[56] Z. Xu, M. Ma, and P. Liu. Self-energy-modified Poisson–Nernst–Planck equations: WKBapproximation and finite-difference approaches.
Phys. Rev. E , 90(1):013307, 2014.[57] J. Kierzenka and L. Shampine. A BVP solver based on residual control and the Matlab PSE.