{'fraction_non_alphanumeric': 0.06633499170812604, 'fraction_numerical': 0.032366901126551155, 'mean_word_length': 3.613030862569243, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 36, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multiple solutions of steady-state Poisson-Nernst-Planck equations with steric effects

Experiments measuring currents through single protein channels show unstable currents. Channels switch between 'open' or 'closed' states in a spontaneous stochastic process called gating. Currents are either (nearly) zero or at a definite level, characteristic of each type of protein, independent of time, once the channel is open. The steady state Poisson-Nernst-Planck equations with steric effects (PNP-steric equations) describe steady current through the open channel quite well, in a wide variety of conditions. Here we study the existence of multiple solutions of steady state PNP-steric equations to see if they themselves, without modification or augmentation, can describe two levels of current. We prove that there are two steady state solutions of PNP-steric equations for (a) three types of ion species (two types of cations and one type of anion) with a positive constant permanent charge, and (b) four types of ion species (two types of cations and their counter-ions) with a constant permanent charge but no sign condition. The excess currents (due to steric effects) associated with these two steady state solutions are derived and expressed as two distinct formulas. Our results indicate that PNP-steric equations may become a useful model to study spontaneous gating of ion channels. Spontaneous gating is thought to involve small structural changes in the channel protein that perhaps produce large changes in the profiles of free energy that determine ion flow. Gating is known to be modulated by external structures. Both can be included in future extensions of our present analysis.


Introduction
The Poisson-Nernst-Planck (PNP) equations, a well-known model of ion transport, play a crucial role in the study of many physical and biological phenomena (cf. [3,4,7,8,12,14,16,17,31,38,39,43,47]). Such an important model can be represented by where N is the number of ion species, c i is the distribution function, J P N P i is the flux density, D i is the diffusion constant, and z i is the valence of the ith ion species, respectively. Besides, φ is the electrostatic potential, ε is the dielectric constant, ρ 0 is the permanent (fixed) charge density of the system, k B is the Boltzmann constant, T is the absolute temperature and e is the elementary charge. Due to ionic sizes, steric repulsion may appear in crowded ions of several biological systems like DNAs, ribosomes and ion channels. When ions are crowded in a narrow channel, the PNP equations become unreliable because the ion-size effect becomes important, but the PNP equations represent ions as point particles without size (cf. [1,5,20,21,27,32,35,45]).
To include ion size effects, Eisenberg and Liu modified PNP equations into a complicated system of differential-integral equations with singular integrals that simulate successfully the selectivity of important types of calcium and sodium ion channels (cf. [29]). However, the singular integrals form an extremely singular kernel because of the Lennard-Jones (LJ) potential. Numerical efficiency and theoretical analysis disappear when forced to deal with such singularities (cf. [19,30]). To simplify the model, we truncate the (spatial) frequency range of the LJ potential, find a simpler energy functional from the leading order terms of the energy expansion with suitable scales. We derive the Poisson-Nernst-Planck equations with steric effects called PNP-steric equations (cf. [36]) where flux J i is and g ij = g ji ∼ ij (a i + a j ) 12 is a nonnegative constant depending on ion radii a i , a j and the energy coupling constant ij of the i-th and j-th species ions, respectively (cf. [26]). Note that equations (1.2)-(1.4) can be regarded as a system of reaction-diffusion equations with nonlinear cross-diffusion terms being similar to [9]. Amazingly, these equations are an effective model to simulate the selectivity of ion channels (cf. [26]).
Comparing ( g ij c j is the excess chemical potential of ion species i due to steric effects. Consequently, the excess current I ex = N i=1 z i eJ ex i due to steric effects becomes z i e k B T D i g ij c i ∇c j . (1.5) We shall use the formula (1.5) to calculate the excess currents for multiple solutions of the 1D steady-state PNP-steric equations. We are motivated by the hope-but cannot dare expect-that one solution will correspond to a closed state and the other to an open state, as found in experiments [15] and in simulations [33]. Of course, the current measured through the open state corresponds to the total current, not just the excess currents. The existence of multiple steady (equilibrium) states is important to study transitions between such states which may be related to the gating (switching between open and closed states) and selectivity of ion channels. Multiple steady states can be investigated by finding multiple solutions of the 1D steady-state PNP equations for two types of ion species with three regions of piecewise constant permanent charge under the assumption that the Debye number is large [16]. More general theorems related to multiple solutions of the 1D steady-state PNP equations involving multiple types of ions with multiple regions of piecewise constant permanent charge are discussed in [37]. With only a constant permanent charge, there is only a unique solution of the 1D steady-state PNP equations for multiple types of ions [39,46]. Instead of the 1D steady-state PNP equations, here we study multiple solutions of the 1D steady-state PNP-steric equations with spatially constant permanent charges.
For simplicity, we consider domain as a 1D interval (−1, 1) for (1.2)-(1.4) and set J i = 0, i = 1, · · · , N to get the steady-state PNP-steric equations. Then by (1.4), g ij c j   = 0 for x ∈ (−1, 1) , i = 1, · · · , N , which can be satisfied if g ij c j = 0 for i = 1, · · · , N , (1.6) holds true. Letφ = e k B T φ andg ij = 1 k B T g ij for i, j = 1, · · · , N . Then (1.3) and (1.6) can be transformed into and For notational convenience, we may remove tilde (∼) and denote (1.7) and (1.8) as Equations like (1.9) have been used to interpret bioelectric phenomena in many papers since they were adopted by Hodgkin, Huxley, and Cole (cf. [13,28]). Here we consider the following boundary condition given by where φ 0 (1), φ 0 (−1) are constants and η ε is a non-negative constant. Here φ 0 (±1) and φ (±1) are the extrachannel and intrachannel electrostatic potentials at the channel boundaries, respectively. The coefficient η ε ∼ ε0 εm is governed by the ratio of ε 0 the dielectric constant of the electrolyte solution and ε m the dielectric constant of the membrane (cf [48]). Note that (1.11) is of the Robin boundary condition if η ε > 0; and of the Dirichlet boundary condition if η ε = 0. The Robin boundary condition includes polarization (e.g. dielectric) charges in the bath and/or electrodes which the Dirichlet boundary condition does not. Such charges, induced by and dependent on the electric field play a prominent role in the art of real experiments, because they are important determinants of the background noise and stability of high speed recordings. The theoretical reasons for these practical realities have not been investigated to the best of our knowledge.
As N = 2, the existence, uniqueness and the solution's asymptotic behavior of (1.9)-(1.11) are investigated under non-symmetry breaking condition 0 ≤ g 12 = g 21 ≤ √ g 11 g 22 which implies that solution (c 1 , c 2 ) of (1.9) is uniquely determined by φ (cf. [34]). Hence (1.9) and (1.10) can be reduced to a single differential equation of φ. However, as the symmetry breaking condition g 12 = g 21 > √ g 11 g 22 holds true, solution (c 1 , c 2 ) of (1.9) may not be uniquely determined by φ. In Section 2, we introduce new variables ξ, Σ and transform (1.9) into a quadratic polynomial which can be solved precisely to get explicit formulas and represent two branches of solution curves. Using these explicit formulas, we can then define biological conductance (for that condition) as the biologists do and perform the comparison using formulas like (1.12)-(1.15). Note that the symbol g is used for conductance (units siemens) in biology and this is not equivalent to our g ij . In this paper, we want to study multiple solutions of (1.9)-(1.11) for the cases of N = 3, 4, and g 12 = g 21 , g 34 = g 43 sufficiently large such that symmetry breaking condition g 12 = g 21 > √ g 11 g 22 , g 34 = g 43 > √ g 33 g 44 holds true.
Natural biological solutions always contain at least three species (sodium, potassium, and chloride, and usually calcium). Experiments are often done, however, with just two species (say sodium chloride) along with traces of hydrogen ion, and perhaps other contaminants. Gating occurs in simplified unnatural situations and so we hope to study mathematical solutions in corresponding situations in a separate paper. Now we state the main result of this paper as follows: Theorem 1.1. Let N = 3, z 2 = −z 1 = q ≥ 1, z 3 > 0 and ρ 0 > 0 be a constant. Assume that g 11 = g 22 = g > 0 is fixed and g i3 = g 3i = 0 for i = 1, 2, 3. Then as g 12 = g 21 = z > 0 is sufficiently large, the system of equations (1.9)-(1.11) has two uniformly bounded (to ε) solutions φ A ε and φ B ε such that φ A ε (x) → φ A1,0 and φ B ε (x) → φ B1,0 for x ∈ (−1, 1) as ε → 0, where φ A1,0 and φ B1,0 are two distinct constants.
In most of the "cation" (e.g., sodium, potassium, and calcium) channels, ρ 0 is a negative number. There are regions ('rings') of negative charge and some channels (sodium channel DEKA) have a ring of positive charge as well. Here we assume the positive sign of ρ 0 which may produce the values φ A1,0 and φ B1,0 (see Figure 4 in Section 3.1), and the proof of Theorem 1.1 is given in Section 3.1.
For solutions φ A ε and φ B ε , the corresponding excess currents defined in (1.5) may be denoted as I ex A and I ex B , respectively. Under the same hypotheses of Theorem 1.1 for three ion species, we may use the new variable Σ to derive the following formulas (see Section 5.1): (1.14) and The difference between I ex A,M and I ex B,N may also give various ion flows related to currents observed in channels as they switch (i.e., gate) from one level of current to another.
The rest of this paper is organized as follows: We may solve system (1.9) of algebraic equations with N = 2, z 2 = −z 1 = q ≥ 1 and g 11 = g 22 > 0 in Section 2. Theorem 1.1 and 1.2 are proven in Section 3. The proofs of Lemma 4.1 and 4.5 are given in Section 4, and formulas (1.12)-(1.15) are derived in Section 5.
2 Solutions of (1.9) with N = 2, z 2 = −z 1 = q ≥ 1 and g 11 = g 22 In this section, we study equation (1.9) with N = 2, z 2 = −z 1 = q ≥ 1 and g 11 = g 22 = g which can be denoted as follows: where z = g 12 and g = g 11 = g 22 are positive constants. Physically, g ij ∼ ij (a i + a j ) 12 , where a i is the ion radius of i-th ion species with concentration c i , and ij > 0 is the energy coupling constant between i-th and j-th ion species for i = 1, 2. Note that (2.1) and (2.2) are formulated as a system of algebraic equations. We want to solve these equations and get solutions for (c 1 , c 2 ) as a function of φ. Adding (2.1) and (2.2), we get Now we introduce new variables as follows: Multiplying σ by c 1 , we get a quadratic polynomial of c 1 as follows: and hence by c 1 c 2 = ξ, (c 1 , c 2 ) can be expressed as Hence the solution (c 1 , c 2 ) of (2.1) and (2.2) may be described by two curves A and B parameterized by the total concentration Σ and denoted as and Here Σ z > 0 a critical total concentration is the unique positive solution of Σ 2 = 4e −(g+z)Σ such that concentrations c 1 and c 2 are equal to 1 2 Σ z as the total concentration Σ = Σ z . Then Take (2.6) and (2.7) into (2.1), and let φ A = φ on curve A, and φ B = φ on curve B, respectively. Then and when the total concentration Σ is larger than Σ z . Note that curve A and B are joined only at a single point (c 1 , Then φ A and φ B can be regarded as one variable φ and c 1 − c 2 may become a strictly monotone increasing function of φ. The result is stated as follows: Proof. Suppose 0 < z ≤ g. Then by (2.11) and (2.12), we have (2.14) Here we have used 0 We may combine φ A and φ B as one variable φ (see Figure 1) defined as follows: Hence by (2.14) and inverse function theorem, Σ can be denoted as Σ = Σ(φ) and become a singlevalued function of φ with domain being the entire space R and range [Σ z , ∞) such that Σ(0) = Σ z and (2.13) hold true. The derivative of Σ with respect to φ is Therefore, c 1 − c 2 is strictly monotone increasing to φ and we complete the proof.
When z = g 12 is increased, for example when the ion is divalent like calcium, the profiles of φ A and φ B may lose monotonicity and become oscillatory. It is well known in experiments that calcium has profound and complex effects on the current voltage relations of channels (cf. [2,22]). Suppose z > 1 + g 2 > 0. Then z 2 − g 2 > 1 and there exists a unique Σ c > 0 (because (1 + gΣ) e (g+z)Σ is strictly monotone increasing to Σ > 0) depending on Σ z such that We shall prove that Σ c may be located in the domain of φ A and φ B i.e. Σ c > Σ z > 0 if z is sufficiently large (see Proposition 2.2). By (2.11) and (2.12) . Then Σ c is a unique (global) minimal point of φ A and a unique (global) maximal point of φ B , respectively (see Figure 2). Moreover, by (2.10), (2.16) Figure 3). and Similarly, the inverse image of function φ B consists of another two functions Σ B1 : Thus by (2.8), we may consider two functions of (c 1 − c 2 ) • Σ A1 and (c 1 − c 2 ) • Σ B1 as follows: and Here we have used (2.38) and (2.39).
) are strictly monotone increasing to φ (see (2.21) and (2.22)), then we may use (2.8) to get . Now we claim that if z is sufficiently large, then Σ c > Σ z > 0 i.e. Σ c is located in the domain of φ A and φ B as follows: Proof. Firstly, we claim that g c is well-defined. For any z > 0, we may define a function Hence there exists a unique Σ c,z > 0 such that f z (Σ c,z ) = 0. Let Σ/// z > 0 be the unique solution of and then f z (Σ z ) < 0 as z sufficiently large. Since f z (Σ c,z ) = 0 and f z (Σ z ) < 0 as z sufficiently large, then by (2.27), we have Σ c,z > Σ z as z sufficiently large. Consequently, the set is nonempty and the value g c = inf z∈Z z (defined in (2.25)) is well-defined. Note that the existence of To complete the proof of Proposition 2.2, we need the following result: Claim 1. Suppose f z0 (Σ z0 ) = 0 and Σ z0 > 0 for some z 0 > 1 + g 2 . Then there exist z l , z r > 1 + g 2 and z l < z 0 < z r such that f z (Σ z ) > 0 for z ∈ (z l , z 0 ) and f z (Σ z ) < 0 for z ∈ (z 0 , z r ).
It is obvious that Now we want to prove that where g c = inf z∈Z z. Due to the continuity of f z , (2.29) implies that the set Z is open. Suppose the set Z has two components. Then without loss of generality, we may assume that there exists z a > g c such that Z = (g c , z a ) ∪ (z a , ∞). Hence f za (Σ za ) = 0 and f z (Σ z ) < 0 for z ∈ (g c , z a ) ∪ (z a , ∞). However, Claim 1 implies that f z (Σ z ) > 0 for z ∈ (z l , z a ) which contradicts to f z (Σ z ) < 0 for z ∈ (g c , z a ). Thus the proof of (2.34) is done. On the other hand, Claim 1 also implies that Otherwise, by (2.33), there exists z b ∈ ( 1 + g 2 , g c ) such that f z b (Σ z b ) = 0. Then as for (2.32), we (i) The proof of Proposition 2.2 shows that f z (Σ z ) > 0 for 0 < z < g c and f z (Σ z ) < 0 for z > g c (see (2.34) and (2.35)). Hence by the continuity of f z , f gc (Σ gc ) = 0.
3 Proof of Theorem 1.1 and 1.2

Proof of Theorem 1.1
In this section, we study multiple solutions of the system of equations (1.9)-(1.11) with N = 3 and the following assumptions: Then we may get solutions of (1.9) by solving and let c 3 = e −z3φ .
Proof. Without loss of generality, we may assume the function f satisfying (F1). Replacing φ by φ + c, we may assume c = 0 and f (0) = 0 in the whole proof for notational convenience. Since the domain of the function f is only [A, M ], then we firstly extend it smoothly to the entire real line R in order to use the standard direct method to get the existence of solution φ. Hence we may temporarily assume the function f as a smooth and strictly monotone increasing function on R. Actually, such an assumption can be ignored because of (4.5).
To prove Lemma 4.1, we need the following Proposition: x a ) and increasing in (x a , 1).
Notice that both Case I and II produce contradiction. Similarly, the condition φ(1) > 0 and φ (1) ≤ 0 and the other condition φ(1) ≤ 0 also result in contradiction, respectively. Therefore, we may complete the proof of Claim I.
We divide three cases to complete the proof as follows: Case I. The solution φ = φ(x) is monotone decreasing to x i.e. φ (x) ≤ 0 for x ∈ (−1, 1).
Using the equation εφ = f (φ) and the condition (a), we have Case II. The solution φ = φ(x) is monotone increasing to x i.e. φ (x) ≥ 0 for x ∈ (−1, 1). As for the argument of Case I, we obtain Case III. The solution φ = φ(x) has a local minimum point at x 0 ∈ (−1, 1) such that φ (x 0 ) = 0 and φ (x 0 ) > 0. Note that since εφ = f (φ) ≥ f (A) > 0, it is impossible to have any local maximum point in (−1, 1). By the equation εφ = f (φ) and the condition (a), we have On the other hand, Therefore, as ε > 0 sufficiently small, (4.8) and (4.9) provide a contradiction and we may complete the proof of Lemma 4.5.