• 沒有找到結果。

model: PNP-steric equations

N/A
N/A
Protected

Academic year: 2022

Share "model: PNP-steric equations"

Copied!
23
0
0

加載中.... (立即查看全文)

全文

(1)

國立臺灣大學 數學系 預印本 Department of Mathematics, National Taiwan University

www.math.ntu.edu.tw/ ~ mathlib/preprint/2013- 07.pdf

A new approach to the Lennard-Jones potential and a new model: PNP-steric equations

Tai-Chia Lin and Bob Eisenberg

April 19, 2013

(2)

model: PNP-steric equations

Tai-Chia Lin

and Bob Eisenberg

Abstract

A class of approximate Lennard-Jones (LJ) potentials with a small parameter is found whose Fourier transforms have a simple asymptotic behavior as the parameter goes to zero.

When the LJ potential is replaced by the approximate LJ potential, the total energy functional becomes simple and exactly the same as replacing the LJ potential by a delta function. Such a simple energy functional can be used to derive the Poisson-Nernst-Planck equations with steric effects (PNP-steric equations), a new mathematical model for the LJ interaction in ionic solutions. Using formal asymptotic analysis, stability and instability conditions for the 1D PNP-steric equations with the Dirichlet boundary conditions for one anionic and cationic species are expressed by the valences, diffusion constants, ionic radii and coupling constants.

This is the first step to study the dynamics of solutions of the PNP-steric equations.

1 Introduction

The Lennard-Jones (LJ) potential, a well-known mathematical model for the interaction be- tween a pair of ions, has important applications in many fields of biology, chemistry and physics (cf. [46]). Such a capable model can be represented by

Ψ (x) = C1

r12 −C2

r6 for r =|x| > 0, x ∈ Rd,

where C1, C2 are positive constants related to finite ion size, and d≤ 3 is the spatial dimension (cf. [36]). The inverse twelfth-power term is the repulsive term of Ψ and the inverse sixth-power term is the attractive term of Ψ. The LJ potential Ψ can be extended to an l− m LJ potential given as follows:

Ψl,m(x) = C1

rm −C2

rl for r =|x| > 0, x ∈ Rd,

where l and m are any positive constants with m > l > d. All the mathematical arguments of the LJ potential Ψ here can easily be generalized to the l− m LJ potential for m > l > d, C1> 0, and C2≥ 0.

To compute the energy of ions interacting by the LJ potential, the following energy functional is considered:

ELJ[ci, cj] =

∫∫

Rd×Rd

Ψ (x− y) ci(x) cj(y) dxdy , (1.1) for nonnegative functions ci and cj which denote the distribution (concentration) functions of the ith and jth ion species (cf. [29]). Similar energy functionals for Coulomb interactions can be found in [9]. Note that the energy functional ELJ is for two ion species if i ̸= j, but for only one ion species if i = j. Because the LJ potential Ψ is singular at the origin, the functional ELJ is a singular integral defined by

ELJ[ci, cj] = lim

σ→0

∫∫

Rd×Rd

Ψχσ(x− y) ci(x) cj(y) dxdy , (1.2)

Department of Mathematics, National Taiwan University, Taida Institute for Mathematical Sciences (TIMS), No.1, Sec.4, Roosevelt Road, Taipei 106, Taiwan, email: tclin@math.ntu.edu.tw

Department of Molecular Biophysics Physiology Rush Medical Center, 1653 West Congress, Parkway, Chicago, IL 60612, USA, email: beisenbe@rush.edu

1

(3)

for ci, cj ∈ L2+

(Rd)

={

f ∈ L2( Rd)

: f≥ 0}

, where χσ= χσ(z) is the characteristic function of the exterior ball{

z∈ Rd: |z| > σ}

. In [16, 17, 18, 28], the energy functional ELJ with C2= 0 is used, but here both C2= 0 and C2> 0 are considered.

Conventionally, the system of Poisson-Nernst-Planck (PNP) equations, a model of ion transport, plays a crucial role in the study of many physical and biological phenomena (cf. [2, 3, 6, 7, 10, 12, 14, 15, 30, 32, 40, 41, 42, 45, 48, 52, 53, 59, 60]). However, 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, 4, 8, 11, 18, 21, 22, 26, 33, 37, 38, 50, 51, 56, 58, 61]). Hence the PNP equations need to be modified in order to describe solutions where the ion-size effect is important. Biological solutions are mixtures containing divalents in which ion size effects are always important (cf. [17, 19, 20]).

To modify the PNP equations, many efforts have been made to combine the energy functional of the PNP equations with the other exclusion terms which may come from liquid state theory and density functional theory [5, 23, 34, 38, 40, 47]. Being related to liquid state theory, the Lennard-Jones (LJ) potential is often used as an approximate model of the van der Waals force (cf. [29, 46]). However, it seems that no one had ever used the LJ potential to modify the PNP equations before the pioneering works of Eisenberg and Liu who derived the PNP equations with size effects i.e. equations (2.6)-(2.8) by combining the repulsive term of the LJ potential and the energy functional of the PNP equations (cf. [16, 28]).

Equations (2.6)-(2.8) importantly generalize the PNP equations and numerically simulate the selectivity of ion channels which can not be obtained by solving the PNP equations alone (cf. [62, 63, 64, 65, 66]). Nevertheless, due to imposing the LJ potential, equations (2.6)-(2.8) become a complicated system of differential-integral equations having no numerical efficiency (cf. [17, 28]) and allowing no theoretical result either. The goal of this paper is to approximate the LJ potential and simplify equations (2.6)-(2.8) into the PNP-steric equations, a new mathematical model for the LJ interaction in ionic solutions. Instead of singular integrals of equations (2.6)-(2.8), the PNP- steric equations are composed of the PNP equations and nonlinear differential terms with coupling constants. Numerical simulations of the PNP-steric equations are presented in [25], which shows the numerical efficiency to simulate the selectivity of ion channels previously studied in Monte Carlo simulations with results comparable to a wide range of experiments.

The main difficulty of numerical simulation of the PNP-steric equations is how to choose cou- pling constants suitably. This motivates us to study a simple case (one anion and one cation species) of the PNP-steric equations with the Dirichlet boundary conditions using suitable asymp- totic expansions to see the effect of coupling constants gij’s. The formal asymptotic analysis gives the stability and instability conditions represented by gij’s which can be regarded as the first step in a series of analyses. (see Section 4 and 5). More theoretical results will be done soon which may be useful for the choice of coupling constants in order to do further numerical simulations on ion channels.

This paper has two major parts: one is the approximation of the LJ potential using band-limited functions and the other is the stability and instability conditions for the PNP-steric equations with the Dirichlet boundary conditions using asymptotic expansions. The approximate LJ potentials are introduced in Section 2.1 and the detailed mathematical arguments are stated in Section 3.

The PNP-steric equations are derived in Section 2.2 and the stability and instability conditions for the PNP-steric equations are proved in Section 4 and 5, respectively.

2 Preliminaries

2.1 The approach to Lennard-Jones potential

Because the LJ potential Ψ does not have a Fourier transform, it is difficult to study the energy functional ELJ directly. When the spatial frequency variable ξ is bounded, the Fourier transform dΨχσ(ξ) tends to infinity as σ goes to zero, but by the Riemann-Lebesgue Lemma,

|ξ|→∞lim Ψχdσ(ξ) = 0 for all σ > 0. Hence the asymptotic behavior of the Fourier transform dΨχσ(ξ) is dominated at bounded spatial frequencies and negligible at high spatial frequencies. This gives a reason to cut off high spatial frequencies in order to see the asymptotic behavior of the Fourier transform dΨχσ(ξ). From a physical and biological point of view, it seems obvious that particularly

(4)

small spatial features cannot be used by evolution to produce biological function and so should be absent in analysis. The high spatial frequency cut-off function φσ is a band-limited function defined by

φσ(x) = (1− χσ−γ(ξ))v, (2.1)

for x, ξ∈ Rd, where v denotes the inverse Fourier transform. Obviously, ˆ

φσ(ξ) = 1− χσ−γ(ξ) =

{ 1 if |ξ| ≤ σ−γ,

0 if |ξ| > σ−γ, (2.2)

for ξ ∈ Rd, where hat denotes the Fourier transform, 0 < γ < 1 is a constant independent of σ and ξ.

Generically, band-limited functions can handle spatial information locally in frequency domains and play important roles in Fourier analysis of mathematics and have many applications to en- gineering, physics and statistics (cf. [31, 35, 55]). Physically, it is obvious that the band-limited approximation is reasonable if not inevitable. After all, derivations of PNP replace Maxwell’s equations with Poisson’s equation, and neglect ’capacitance to ground’ (capacitive coupling be- tween ionic solutions and nearby ground planes) always present in experiments. These and many other effects at high frequencies are better attenuated to zero than approximated irrationally (i.e., without known error bounds) as they are in treatments with unlimited bandwidth. Here a new ap- proach to the LJ potential Ψ uses the spatially band-limited function φσto define the approximate potential Ψ as follows:

Ψσ(z) = (Ψχσ) ⋆ φσ(z) for z∈ Rd, (2.3) where the asterisk is the standard convolution, φσis the spatially band-limited function defined in (2.1) and χσ is the characteristic function of the exterior ball{z ∈ Rd:|z| > σ}. Note that χσ(z) and 1− χσ−γ(ξ) are the characteristic functions of {

z∈ Rd:|z| > σ} and {

ξ∈ Rd:|ξ| ≤ σ−γ} , respectively, and both of them extend to the entire spaceRdas the small parameter σ goes to zero.

As σ goes to zero, χσ∼ 1, 1 − χσ−γ ∼ 1, and

Ψcσ(ξ) = dΨχσ(ξ)φcσ(ξ) = dΨχσ(ξ) [1− χσ−γ(ξ)]∼ bΨ (ξ) . which implies Ψσ∼ Ψ. Hence formally,

Rd

ci(x) (Ψσ⋆ cj) (x) dx∼

Rd

ci(x) (Ψ ⋆ cj) (x) dx = ELJ[ci, cj] This shows how we approximate the energy functional ELJ.

The approximate energy functional ELJ,σ is defined by ELJ,σ[ci, cj] =

∫∫

Rd×Rd

Ψσ(x− y) ci(x) cj(y) dxdy =

Rd

ci(x) (Ψσ ⋆ cj) (x) dx.

As σ goes to zero, the functional ELJ,σ tends to the functional ELJ if the following hypothesis holds:

(H) lim

σ→0+

[∫

RdΨχdσ(ξ) [1− χσ−γ(ξ)]cbi(ξ)cbj(ξ)dξ−

RdΨχdσ(ξ)cbi(ξ)cbj(ξ)dξ ]

= 0.

(see Proposition 3.1). Here the meaning of ‘approximate’ is different from that of conventional approximation theory. Note that the characteristic function 1−χσ−γ truncates the high frequencies

|ξ| > σ−γ, but still preserves the spatial frequencies of order|ξ| ∼ σ−θ (for all 0 < θ < γ) tending to infinity as σ goes to zero. By standard theorems of Fourier analysis (cf. [54]), the functional ELJ,σsatisfies ELJ,σ[ci, cj] =∫

Rdcbi(ξ) \Ψσ⋆ cj(ξ) dξ =

Rdcbi(ξ) cΨσ(ξ)cbj(ξ) dξ On the other hand, (2.2) and (2.3) imply cΨσ(ξ)∼ C1 ωd

12−dσd−12 as σ goes to zero, where ωd is the surface area of d dimensional unit ball (see the proof of Proposition 3.2). Hence

ELJ,σ[ci, cj]∼ C1

ωd

12− dσd−12

Rdcbi(ξ)cbj(ξ) dξ = C1

ωd

12− dσd−12

Rd

ci(x)cj(x) dx . This shows that the functional ELJ,σis asymptotically close to another functional ˜ELJ,σas follows:

ELJ,σ[ci, cj]∼ ˜ELJ,σ[ci, cj] = C1Sσ

Rd

ci(x)cj(x) dx , (2.4)

(5)

where Sσ = 12ω−dd σd−12. Thus the functional ˜ELJ,σ can be regarded as an approximate energy functional to ELJ. Note that the constant C1 comes from the repulsive term of the LJ potential, but the constant C2 (for the attractive term of the LJ potential) does not affect the leading term of the asymptotic behavior of cΨσ. This may support the work of [16, 17, 18, 28] which only use the repulsive term of the LJ potential to describe ionic interactions. However, the effects of the attractive terms remain to be investigated. They may have been selected by evolution to produce qualitative behavior of importance to biology. Note that the functional ˜ELJ,σis much simpler than the functional ELJ and it can be expressed as

E˜LJ,σ[ci, cj] =

∫∫

Rd×Rd

C1Sσδ0(x− y)ci(x) cj(y) dxdy

being the same as replacing the LJ potential Ψ by a delta function C1Sσδ0in the energy functional ELJ, where δ0(·) is the standard delta function concentrating at the origin. Such a simple energy functional ˜ELJ,σ can be used to derive the PNP-steric equations as a new model of ionic solutions.

2.2 PNP-steric equations

The Poisson-Nernst-Planck (PNP) equations consist of the Nernst-Planck equations coupled with the Poisson equation being expressed as follows:









∂ci

∂t =kDi

BT∇ ·(

ciδEδcpnpi )

, i = 1,· · · , N,

−∇ · (ε∇ϕ) = ρ0+

N i=1

zieci,

where Epnpis the energy functional of the PNP equations given by Epnp[c1,· · · , cN, ϕ] :=

Rd

[ kBT

N i=1

cilog ci+12 (

ρ0+

N i=1

zieci )

ϕ ]

dx

Here N is the number of ion species, ci is the distribution function, Di is the diffusion constant, and zi 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, kB is the Boltzmann constant, T is the absolute temperature and e is the elementary charge. More precisely, the PNP equations are denoted as









∂ci

∂t = Di∇ ·(

∇ci+kzie

BTci∇ϕ)

, i = 1,· · · , N,

−∇ · (ε∇ϕ) = ρ0+

N i=1

zieci

(2.5)

To include the hard sphere repulsion of ions, the energy functional Epnpis modified by adding the energy functional

N i,j=1

ELJ[ci, cj] with the constants C1= 12ϵij(ai+ aj)12 and C2= 0, where ciand cjare the distribution (concentration) functions of the ith and jth ion species with the radii aiand aj, respectively. Then the modified energy functional Empnp becomes

Empnp[c1,· · · , cN, ϕ] :=

Rd

( kBT

N i=1

cilog ci+1 2

( ρ0+

N i=1

zieci

) ϕ

) dx

+1 2

N i,j=1

∫∫

Rd×Rd

ϵij(ai+ aj)12

|x − y|12 ci(x) cj(y) dxdy ,

Using energy variational analysis (cf. [16, 28]), the modified PNP equations are the following equations

∂ci

∂t = kDi

BT∇ ·(

ciδEδcmpnpi )

, i = 1,· · · , N,

(6)

i.e. ∂ci

∂t +∇ · Ji = 0 , i = 1,· · · , N, (2.6) coupled with the Poisson equation

−∇ · (ε∇ϕ) = ρ0+

N i=1

zieci, (2.7)

where flux Ji is

Ji=−Di∇ci−Dici

kBTzie∇ϕ − Dici

kBT

N j=1

Rd

ϵij(ai+ aj)12

|x − y|12 cj(y) dy . (2.8) However, equations (2.6)-(2.8) are difficult to investigate theoretically because they are partial differential-integral equations with singular integrals. Moreover, due to the effect of high (Fourier) frequencies, the numerical computations of equations (2.6)-(2.8) may lose accuracy and become inefficient (cf. [28]).

Instead of equations (2.6)-(2.8), a simple model can be derived by replacing the energy func- tional

N i,j=1

ELJ[ci, cj] (with C1=12ϵij(ai+aj)12and C2= 0) by the approximate energy functional

N i,j=1

E˜LJ,σ[ci, cj] = 1

2ϵij(ai+ aj)12Sσ

Rd

ci(x)cj(x) dx

defined in (2.4). Then the energy functional Empnp can be approximated by

Eσ[c1,· · · , cN, ϕ] :=

Rd

( kBT

N i=1

cilog ci+1 2

( ρ0+

N i=1

zieci

) ϕ

) dx

+1 2

N i,j=1

ϵij(ai+ aj)12Sσ

Rd

ci(x) cj(x)dx ,

which is much simpler than the energy functional Empnp. Using energy variational analysis (cf. [16, 28]), a new model called the PNP-steric equations is expressed by

∂ci

∂t =kDi

BT∇ ·(

ciδEδcσi)

, i = 1,· · · , N,

i.e. ∂ci

∂t +∇ · Ji= 0 , (2.9)

coupled with the Poisson equation

−∇ · (ε∇ϕ) = ρ0+

N i=1

zieci, (2.10)

where flux Ji is

Ji=−Di∇ci−Dici

kBTzie∇ϕ − Dici kBTSσ

N j=1

ϵij(ai+ aj)12∇cj. (2.11)

Here the symmetry ϵij = ϵjihas been assumed for notation convenience. The PNP-steric equations are of convection-diffusion type with the following energy dissipation law:

d

dtEσ[c1,· · · , cN, ϕ] =−

Rd

N i=1

Dici

kBT|∇ (kBT log ci+ zieϕ + µi)|2dx

(7)

where µi= δcδ

i

N j,k=1

E˜LJ,σ[cj, ck] = Sσ

N j=1

ϵij(ai+ aj)12cj is the chemical potential. On the other hand, the PNP-steric equations have more nonlinear differential terms than the PNP equations so they can simulate the selectivity of ion channels efficiently (cf. [25]). Note that the selectivity of ion channels can not be found by using the PNP equations. This shows that the PNP-steric equations are more capable than the PNP equations and much simpler than equations (2.6)-(2.8).

For the case of two species ions (i.e. N = 2) with one anionic and cationic species, the index j = 1, 2 is replaced by j = n, p for notation convenience and the PNP-steric equations (2.9)-(2.11) are presented as

∂cn

∂t =Dn

[

∇ · (

∇cn+ zne kBTcn∇ϕ

)

+ ˜Sσ∇ · (gnncn∇cn+ gnpcn∇cp) ]

, (2.12)

∂cp

∂t =Dp

[

∇ · (

∇cp+ zpe kBTcp∇ϕ

)

+ ˜Sσ∇ · (gppcp∇cp+ gnpcp∇cn) ]

, (2.13)

−∇ · (ε∇ϕ) =ρ0+ znecn+ zpecp, (2.14)

for x∈ Rd, t > 0, where the function cn is for anion concentration, cp is for cation concentration, S˜σ = k1

BT Sσ, gnp = ϵ12(a1+ a2)12, gnn = ϵ11(2a1)12 and gpp = ϵ22(2a2)12. Due to the spatial dimension d≤ 3, the constant ˜Sσ∼ σd−12 becomes a large quantity tending to infinity as σ goes to zero. Let ˜ε = keB2Tε, ˜δ = 1/ ˜Sσ, τ = t/˜δ, ˜ϕ(x, τ ) = ke

BTϕ(x, t) and ˜cj(x, τ ) = cj(x, t) for j = n, p.

Then the equations (2.12)-(2.14) are transformed into 1

Dn

∂˜cn

∂τδ∇ ·(

∇˜cn+ zn˜cn∇ ˜ϕ)

+∇ · (gnn˜cn∇˜cn+ gnp˜cn∇˜cp) , (2.15) 1

Dp

∂˜cp

∂τδ∇ ·(

∇˜cp+ zpc˜p∇ ˜ϕ)

+∇ · (gpp˜cp∇˜cp+ gnp˜cp∇˜cn) , (2.16)

−∇ · (˜ε∇ ˜ϕ) =˜ρ0+ zn˜cn+ zp˜cp, (2.17) for x∈ Rd, τ > 0, where ˜ρ0=ρe0. Note that ∂c∂tj = ˜δ−1 ∂˜c∂τj for j = n, p, and τ ∼ 1 is equivalent to t∼ ˜δ = 1/ ˜Sσ being a small time scale. For notation convenience, removing tilde and replacing τ by t, the equations (2.15)-(2.17) become

1 Dn

∂cn

∂t =δ∇ · (∇cn+ zncn∇ϕ) + ∇ · (gnncn∇cn+ gnpcn∇cp) , (2.18) 1

Dp

∂cp

∂t =δ∇ · (∇cp+ zpcp∇ϕ) + ∇ · (gppcp∇cp+ gnpcp∇cn) , (2.19)

−∇ · (ε∇ϕ) =ρ0+ zncn+ zpcp, (2.20)

for x∈ Rd, t > 0.

Instead of the entire space Rd, here the spatial domain is considered as an one-dimensional (d = 1) interval (−1, 1), and then the equations (2.18)-(2.20) are changed as the following equations:

1 Dn

∂cn

∂t

∂x· (∂cn

∂x + zncn

∂ϕ

∂x )

+

∂x (

gnncn

∂cn

∂x + gnpcn

∂cp

∂x )

, for x∈ (−1, 1) , t > 0 , 1

Dp

∂cp

∂t

∂x (∂cp

∂x + zpcp∂ϕ

∂x )

+

∂x (

gppcp∂cp

∂x + gnpcp∂cn

∂x )

, for x∈ (−1, 1) , t > 0 ,

−ε∂2ϕ

∂x2 0+ zncn+ zpcp, for x∈ (−1, 1) , t > 0 .

(8)

Setting ε = 1 and ρ0= 0, these equations become 1

Dn

∂cn

∂t

∂x (∂cn

∂x + zncn∂ϕ

∂x )

+

∂x (

gnncn∂cn

∂x + gnpcn∂cp

∂x )

, for x∈ (−1, 1) , t > 0 , (2.21) 1

Dp

∂cp

∂t

∂x (∂cp

∂x + zpcp

∂ϕ

∂x )

+

∂x (

gppcp

∂cp

∂x + gnpcp

∂cn

∂x )

, for x∈ (−1, 1) , t > 0 , (2.22)

−∂2ϕ

∂x2 =zncn+ zpcp, for x∈ (−1, 1) , t > 0 ,

(2.23) where δ is a small parameter tending to zero. For simplicity, we consider the following Dirichlet boundary conditions:





cn=Bzn1, cp= Bz1

p, ϕ = B3 as x = 1 , cn=Bzn2, cp= Bz2

p, ϕ = B4 as x =−1 ,

(2.24)

where Bk, k = 1,· · · , 4 are constants. Here we only focus on the equations (2.21)-(2.23) with the boundary condition (2.24) but not (2.18)-(2.20) on the entire space. Using the formal asymptotic analysis, we may get the stability and instability conditions of the equations (2.21)-(2.23) with the boundary condition (2.24) which show the effect of coupling constants gij’s.

To see the solution of the equations (2.21)-(2.23), the following asymptotic expansions are used:

cn = cn,0+ δcn,1+ δ2cn,2+· · · , cp= cp,0+ δcp,1+ δ2cp,2+· · · , ϕ = ϕ0+ δϕ1+ δ2ϕ2+· · · . Then the zeroth order solution (cn,0, cp,0, ϕ0) satisfies

1 Dn

∂cn,0

∂t =

∂x (

gnncn,0

∂cn,0

∂x + gnpcn,0

∂cp,0

∂x )

, (2.25)

1 Dp

∂cp,0

∂t =

∂x (

gppcp,0

∂cp,0

∂x + gnpcp,0

∂cn,0

∂x )

, (2.26)

−∂2ϕ0

∂x2 = zncn,0+ zpcp,0, (2.27)

for x∈ (−1, 1) and t > 0. As gnp= 0, both (2.25) and (2.26) become the standard porous medium equations (cf. [57]) hence for gnp̸= 0, the equations (2.25) and (2.26) can be regarded as a coupled system of porous medium equations.

The steady state equations of (2.25) and (2.26) are denoted as (gnncn,0cn,0+ gnpcn,0cp,0)

(x) = 0 , (2.28)

(gppcp,0cp,0+ gnpcp,0cn,0)

(x) = 0 , (2.29)

for x ∈ (−1, 1), where a prime mark (′) is denoted as differentiation to the spatial variable x.

Electroneutrality (which means the charge of anions is equal to that of cations) holds in most biological systems (cf. [67]). A particular assumption called pointwise electroneutrality which means the charge of anions is everywhere equal to that of cations appears in the zeroth order equation [24] as a pleasingly natural physical approximation. To find solutions of (2.28)-(2.29) with pointwise electroneutrality, it is natural to assume that cn,0and cp,0have the following form cn,0(x) = zpw(x) , cp,0(x) =−znw(x) for x∈ (−1, 1) , (2.30) so the total charge of cn,0 and cp,0 becomes zero i.e. zncn,0(x) + zpcp,0(x) = 0 for x ∈ (−1, 1).

Note that zn and zp are the associated valences, the charge on an individual ion. To solve the

(9)

equations (2.28) and (2.29), the function w = w(x) must satisfy ( w2)′′

(x) = 0 for x∈ (−1, 1), and hence w = w(x) =√

K1x + K0, where K0 and K1 are constants determined by the Dirichlet boundary conditions as follows:

−zncn,0(1) = zpcp,0(1) = B1, −zncn,0(−1) = zpcp,0(−1) = B2, (2.31) and B1and B2are the positive constants for the Dirichlet boundary conditions in (2.24). Because w(x) =

K1x + K0 for x ∈ [−1, 1], then by (2.30) and (2.31), it is obvious that

K1+ K0 = w(1) = −zB1

nzp and

K0− K1= w(−1) = −zBn2zp which imply

K1= B12− B22

2zn2z2p and K0=B21+ B22 2zn2z2p . Hence

w(x) =

(B21− B22)x + B12+ B22

2|znzp| > 0 for x∈ [−1, 1] , (2.32) being non-constant if B1̸= B2. Combining (2.27) and (2.30), the electric potential ϕ0satisfies the following equation:

−ϕ′′0(x) = 0 for x∈ (−1, 1) . (2.33)

The equation (2.33) is a standard differential equation solved as follows:

ϕ0(x) = a0x + b0 for x∈ (−1, 1) , (2.34) where a0 and b0 are constants determined by the boundary conditions of ϕ0 given in (2.24). The first order solution (cn,1, cp,1, ϕ1) = (cn,1(x, t), cp,1(x, t), ϕ1(x, t)) satisfies

1 Dn

∂cn,1

∂t =

∂x (∂cn,0

∂x + zncn,0∂ϕ0

∂x )

+

∂x [

gnn∂(cn,0cn,1)

∂x + gnp (∂cp,0

∂x cn,1+ cn,0∂cp,1

∂x )]

, (2.35) 1

Dp

∂cp,1

∂t =

∂x (∂cp,0

∂x + zpcp,0

∂ϕ0

∂x )

+

∂x [

gpp

∂(cp,0cp,1)

∂x + gnp

(∂cn,0

∂x cp,1+ cp,0

∂cn,1

∂x )]

, (2.36)

−∂2ϕ1

∂x2 = zncn,1+ zpcp,1, (2.37)

for x∈ (−1, 1), t > 0. By (2.30), the equations (2.35) and (2.36) become 1

Dn

∂cn,1

∂t = zp(w+ zn0)+

∂x [

gnnzp

∂(w cn,1)

∂x + gnp

(

−znwcn,1+ zpw∂cp,1

∂x )]

, (2.38) 1

Dp

∂cp,1

∂t =−zn(w+ zp0)+

∂x [

−zngpp

∂(w cp,1)

∂x + gnp

(

zpwcp,1− znw∂cn,1

∂x )]

, (2.39) for x∈ (−1, 1), t > 0.

Let ψ = zncn,1+ zpcp,1and φ =−zncn,1+ zpcp,1i.e. cn,1= 2z1

n(ψ− φ) and cp,1=2z1

p(ψ + φ).

Then after some algebraic calculations, the equations (2.38) and (2.39) can be transformed into the equations for ψ and φ as follows:

∂ψ

∂t = (Dn− Dp) znzpw′′(x) + znzp(Dnzn− Dpzp) [w(x)ϕ0(x)] (2.40) + ( ˜dnp− ˜d)∂2[w(x)φ]

∂x2 + (¯g + ¯gnp)

∂x[w(x)ψ] + (¯g− ¯gnp)

∂x [

w(x)∂ψ

∂x ]

, and

∂φ

∂t =− (Dn+ Dp) znzpw′′(x)− znzp(Dnzn+ Dpzp) [w(x)ϕ0(x)] (2.41) + (¯gnp+ ¯g)∂2[w(x)φ]

∂x2 + ( ˜dnp− ˜d)

∂x[w(x)ψ]− ( ˜dnp+ ˜d)

∂x [

w(x)∂ψ

∂x ]

,

(10)

where ¯gnp, ˜dnp, ¯g and ˜d are constants defined as follows:



















¯

gnp =12gnp(Dpzp− Dnzn) , d˜np =12gnp(Dpzp+ Dnzn) ,

¯

g =12(Dnzpgnn− Dpzngpp) , d˜ =12(Dnzpgnn+ Dpzngpp) .

(2.42)

The equations (2.40) and (2.41) depend on valences zi, diffusion constants Di and coefficients gij∼ (ai+ aj)12(i, j = n, p) related to ionic radii. In particular, these equations can represent the difference between NaCl, KCl and CaCl2.

For notational convenience, the equations (2.40) and (2.41) can be expressed as

∂t

ψ

φ

 =

f1

f2

 +(

L11 L12

L21 L22

)  ψ

φ

 , (2.43)

where f1 and f2 are external force functions given by

f1(x) = (Dn− Dp) znzpw′′(x) + znzp(Dnzn− Dpzp) [w(x)ϕ0(x)] , (2.44) f2(x) =− (Dn+ Dp) znzpw′′(x)− znzp(Dnzn+ Dpzp) [w(x)ϕ0(x)] , (2.45) and Lij, i, j = 1, 2 are differential operators defined by

L11ψ = (¯g + ¯gnp)

∂x[w(x)ψ] + (¯g− ¯gnp)

∂x [

w(x)∂ψ

∂x ]

, (2.46)

L12φ = ( ˜dnp− ˜d)∂2[w(x)φ]

∂x2 , (2.47)

L21ψ = ( ˜dnp− ˜d)

∂x[w(x)ψ]− ( ˜dnp+ ˜d)

∂x [

w(x)∂ψ

∂x ]

, (2.48)

L22φ = (¯gnp+ ¯g)∂2[w(x)φ]

∂x2 . (2.49)

When ¯g > ¯gnp, the asymptotic stability of system (2.43) (i.e. the equations (2.40) and (2.41)) is proved in Corollary 4.2 (see Section 4). Thus the equations (2.38) and (2.39) have asymptotic stability if the condition ¯g > ¯gnp holds. However, if gnn = gpp = gnp = 0, then system (2.43) becomes ψt= f1and φt= f2which imply (ψ, φ) = (ψ0+ tf1, φ0+ tf2) for t > 0, where ψ0= ψ|t=0 and φ0= φ|t=0are initial conditions of ψ and φ, respectively. Hence the asymptotic stability of system (2.43) is gone. This shows that the effect of coupling constants gij’s changes the asymptotic stability of system (2.43).

Note that the operator

( L11 L12

L21 L22

)

=

( L11 0 0 L22

)

becomes diagonal if ˜dnp− ˜d = d˜np+ ˜d = 0 i.e. ˜dnp= ˜d = 0. By (2.42), the condition ˜dnp= ˜d = 0 is equivalent to Dpzp+ Dnzn= 0 and Dnzpgnn+ Dpzngpp = 0 being fulfilled if Dp= Dn, zp=−zn and gnn= gpp. In Section 5, the case of Dp = Dn = D > 0, zp =−zn = z > 0 and gnn= gpp = g > 0 (i.e. symmetry electrolytes) is considered in order to get the instability condition g < gnp.

3 The approximate LJ potential

In this section, we study the approximation of the functional ELJ defined in (1.2) and written as follows:

ELJ[ci, cj] = lim

σ→0

∫∫

Rd×Rd

[(Ψ12− Ψ6) χσ] (x− y) ci(x) cj(y) dxdy, (3.1)

(11)

for ci, cj ∈ L2+

(Rd)

={

f ∈ L2( Rd)

: f≥ 0}

. Here Ψ12(z) = C1|z|−12 is the repulsion term and Ψ6(z) = C2|z|−6 is the attraction term of the LJ potential, and χσ = χσ(z) is the characteristic function of the exterior ball{

z∈ Rd: |z| > σ}

satisfying

χσ(z) =



1 if |z| > σ, 0 if |z| ≤ σ.

The main idea here is to approximate \Ψmχσ(ξ) (m = 6, 12) the Fourier transform of the kernel Ψmχσ by \Ψmχσφˆσ(ξ) in the frequency space, where hat denotes Fourier transform and φσ is a band-limited function satisfying

ˆ

φσ(ξ) = 1− χσ−γ(ξ) for ξ∈ Rd. (3.2) Here 1− χσ−γ is defined by

1− χσ−γ(ξ) =



1 if |ξ| ≤ σ−γ, 0 if |ξ| > σ−γ,

where σ > 0 is a small parameter tending to zero, and 0 < γ < 1 is a constant independent of ξ and σ. Note that 1− χσ−γ denotes the characteristic function of the ball{ξ ∈ Rd :|ξ| ≤ σ−γ} expanding the entire spaceRd as σ goes to zero.

By the standard formulas of Fourier analysis, it is obvious that

Ψ\mχσ(ξ) ˆφσ(ξ) =[(Ψm\χσ) ⋆ φσ](ξ) , (3.3) for ξ ∈ Rd, where the asterisk denotes convolution. The approximate kernel of Ψm is given as follows:

Km,σ(z) = (Ψmχσ) ⋆ φσ(z) for z∈ Rd, m = 6, 12 , (3.4) and the approximate energy functional ELJ,σ given by

ELJ,σ[ci, cj] :=

∫∫

Rd×Rd

(K12,σ− K6,σ) (x− y)ci(x)cj(y)dxdy, (3.5)

for ci, cj ∈ L2+

(Rd)

. By standard theorems of Fourier analysis (cf. [54]), it is easy to check that the functions Ψmχσ and φσare of L2(Rd)∩ L(Rd) which implies Km,σ∈ L2(Rd)∩ L(Rd). The following proposition is for the approximation of ELJ,σ to ELJ.

Proposition 3.1. Assume ci, cj ∈ L2+(Rd) satisfy the following hypothesis:

(H) lim

σ→0+

[∫

RdΨχdσ(ξ) [1− χσ−γ(ξ)]cbi(ξ)cbj(ξ)dξ−

RdΨχdσ(ξ)cbi(ξ)cbj(ξ)dξ ]

= 0 , where 0 < γ < 1 is a constant independent of σ and ξ. Then

σlim→0+{ELJ,σ[ci, cj]− ELJ[ci, cj]} = 0 . (3.6) Note that the hypothesis (H) is achievable at least for functions ci’s satisfying cbi(ξ) = 0 for

|ξ| > σ−γ and i = 1,· · · , N, which means all the high frequencies of ci’s have been cut off. Due to the strong singularity of Ψ, dΨχσ(ξ) has no asymptotic behavior like the right side of (3.12) as σ goes to zero, especially for|ξ| ∼ σ−1 (see Remark 3.3). This motivates us to replace dΨχσ(ξ) by Ψχdσ(ξ) [1− χσ−γ(ξ)] which is a kind of truncation on the ξ variable and has a simple asymptotic formula (3.12). The truncation may lose the effect of high frequencies|ξ| > σ−γ(0 < γ < 1), but still involve a large part of the steric effects because it keeps the effect of frequencies like|ξ| ∼ σ−α for 0 < α < γ tending to infinity as σ goes to zero. Please note that the main goal of our works is to simplify the model of Liu and Eisenberg [16, 28]. Here we present a simplified model which is easy to study and useful to understand the selectivity of ion channels [25].

參考文獻

相關文件

The five-equation model system is composed of two phasic mass balance equations, the mixture momentum equation, the mixture total energy equation, and an evolution equation for

We have derived Whitham equations to the SGN model and show that they are strictly hyperbolic for arbitrary wave amplitudes, i.e., the corresponding periodic wave trains

For a polytomous item measuring the first-order latent trait, the item response function can be the generalized partial credit model (Muraki, 1992), the partial credit model

Especially, the additional gauge symmetry are needed in order to have first order differential equations as equations

Abstract In this paper, we consider the smoothing Newton method for solving a type of absolute value equations associated with second order cone (SOCAVE for short), which.. 1

To improve the convergence of difference methods, one way is selected difference-equations in such that their local truncation errors are O(h p ) for as large a value of p as

The finite difference equations will now be applied to solve the problem of a doubly drained clay layer, undergoing one dimensional consolidation. The assumptions made in deriving

The design of a sequential circuit with flip-flops other than the D type flip-flop is complicated by the fact that the input equations for the circuit must be derived indirectly