We first introduce the notation for this model. A dimension defines some physical characteristics, for example, length [L], time [T ], charge [C], poten-tial [V ].
In mathematics, the steady-state flow of ions can be approximated by the DD or Poisson-Nernst—Planck (PNP) model, that is, by partial differential equations for conservation of each ion species and Poisson’s equation for the electrostatic potential φ (x):
∇ · js = 0, js= zsµsEns− Ds∇ns
1 L3· T
(2.1)
−∇ · (∇φ) = qNN +
s
qsns, E= −∇φ
C L3
(2.2) where for each ion species s, js is the number current density (particle flux), e is the proton charge, qs is the ionic charge (for example, qK+ = +1e, qCl− = −1e), zs = qs/e, µs is the mobility coefficient, Ds is the diffusion coefficient, is the dielectric coefficient, and qN is the the protein charge (for example, qN = −4e). The total electric current density (charge flux) is
jelec =
s
qsjs.
C L2· T
(2.3) The physical parameters , µs , and Ds are functions of x.
Figure 1. Diagram of computational region for channel, membrane, and baths.
As an example of a physiologically important channel, we will focus our attention here on the K channel illustrated in Fig. 1 (the KcsA channel structure shown is derived from X-ray crystallography). K channels play a central role in electrical signaling in the nervous system. A typical nerve cell has hundreds of thousands of K channels. The K channel is selective;
i.e., it allows K+ ions to flow freely between the interior and exterior of the cell, but not (for example) Na+ or Ca++ ions. (Cl− ions are prevented from flowing through the K channel by the electrostatic field in the channel.) We will model the channel plus regions of the bath illustrated in Fig. 2 out to a distance where the ion densities and the electrostatic potential take on their asymptotic values in the baths: ns = Nbi , φ = 0 at the left boundary, and ns = Nbi , φ = V at the right boundary.
Figure 2. K channel, membrane, and interior and exterior baths.
Boundary conditions should not be applied at the ends of the channel, due to the fact that boundary layers in ionic charge density and in the electrosta-tic potential develop there, as illustrated in Figs. 4 and 6. The ionic charge densities and the electrostatic potential reach their equilibrium far-field val-ues approximately two Debye lengths into the baths (the Debye lengths for K+ and Cl− are on the order of 1 nm in the baths).
In the one-dimensional approximation to the K channel problem with the channel 3.5 nm long, the baths are represented by conical funnels (Fig. 3) extending 5 nm into the baths and opening at 45◦ angles on either side of the z axis (the opening angle has only a weak effect on the computed solution and total current).
The boundary condition (BC) types in Fig. 2 are defined by ns = Nbi , φ = 0 (interior bath far-field BC)
ns = Nbi , φ = V (exterior bath far-field BC) ns = Nbi , ∂φ∂r = 0 (ambient bath BC)
n·∇ns = 0, n·∇φ = 0 (no-flux BC), where n is a unit normal vector to the boundary.
Figure 3. Diagram of 1D computational region for channel and bath funnels.
In order to have a self-adjoint expression as motivated by [3], we introduce new variable us for ns = exp (csφ) us, where csis a constant to be determined.
With the us, equation (2.1) is reformulated as
js = zsµsEexp (csφ) us− Ds(exp (csφ) ∇us+ csexp (csφ) ∇φus)
= [zsµs(−∇φ) exp (csφ) − Dscsexp (csφ) ∇φ] us
−Dsexp (csφ) ∇us. (2.4)
An useful choice of the unknown constant cs is to eliminate the gradient of us leading to a divergence formulation for the new variable. We thus slove the following equation for cs :
(zsµs(−∇φ) exp (csφ) − Dscsexp (csφ) ∇φ) us = 0 (2.5) which yields cs = −zsDµss = −VzsT, where VT is thermal voltage. Then the new
variable us can therefore be defined as Note that this expression is purely mathematical and very similar to that of the Slotboom variables (22) and (23) in [3]. The equations (2.1) − (2.3) are then reduced to
∇ · js = 0, js= −Dsexp (csφ) ∇us
As noted in [8], the channel problem illustrated in Fig. 1 has an ap-proximate cylindrical symmetry; that is, the functions us, φ depends to an excellent approximation only on r and z, where cylindrical coordinates are denoted by (r, θ, z), with the z axis along the length of the channel. (With cylindrical symmetry, no-flux boundary conditions are imposed along the z axis as well.) If we neglect the dependence of us and φ on r in the chan-nel and the near bath regions, which is a good approximation if the cross
sectional area A (z) varies slowly with z, then we obtain a one-dimensional approximation to the channel problem.
To write down the 1D equations, recall that the definition of the diver-gence of a vector field u [14] is
∇ · u = divu = lim
where V is the volume of an arbitrary shaped region in R3 that includes the point x, S(V ) is the surface of that volume, and the integral is a surface integral with n being the outward normal to that surface. For u(x) =u (z)ez
and V =A (z) ∆z, (2.10) reduces to
by Taylor expanding in z. (Note that the gradient is unchanged: ∇φ = dφ/dz for φ (x) = φ (z).)
Using this definition, we can obtain that
∇ · (−Dsexp (csφ) ∇us) = 1
The steady-state drift diffusion equations thus become in the 1D approx-imation
−1 where the simulation region which includes the channel plus baths runs from 0 to L = lc+2lb and us(x) = us(z) /A (z) , φ (x) = φ (z) , N (x) = N (z) /A (z) , Nbi(x) = Nbi(z) /A (z) . For both the parabolic transport equation (2.12) and the elliptic Poisson’s equation (2.13), Dirichlet boundary conditions on φ and us are imposed at the left and right boundaries.
3 Quantum-corrected drift diffusion model
Motivated by [5], we propose a quantum corrected model for the ionic chan-nel. Based on the density gradient theory which was developed by observing that the electron gas is energetically sensitive not only to its density but also to the gradient of the density. It captures the nonlocality of quantum mechanics to the lowest-order of h2 where h2 is the reduced Planck constant and can be rigorously derived from quantum mechanics [3] [4]. Specifically, a third order derivative term of quantum correction is added to the carrier current density as
where the coefficients bs = 12mh2∗
se are the material parameters measuring the strength of the gradient effects in the gas, m∗s are the ionic effective masses (Chlorine effective mass = 35.453 * atomic effective mass, Potassium effective mass = 39.0983 * atomic effective mass). To alleviate the difficulty in discretization caused by this higher order term, additional variables called the quantum potentials
φqs≡ −2zsbs
∆√ns
√ns
(3.2) have been introduced in [12] and thus can be lumped with the classical drift term to obtain
˜js = −zsqsµsns∇
φ + φqs
− qsDs∇ns. (3.3) We thus have a complete set of PDEs (2.1)—(2.2) and (3.2) describing both DD and DG models with the state variables φ, us and φqs.
PDEs in self-adjoint form are analytically as well as numerically appealing [3] [4]. We now consider the self-adjoint formulation of the above QCDD model and, for this purpose, introduce the following new variables
ζs=√ns. (3.4)
Assuming a Maxwell—Boltzmann energy distribution of carriers, we have the quantum correction expressions of the carriers
ns= exp cs
φ + φqs
us = ζ2s (3.5)
and rewrite the quantum potentials as φqs = 1
c ln ζ2s
− 1
c ln (us) − φ. (3.6)
For Eq. (2.2) we have
Substituting (3.5) into the current equation (3.3), we obtain
˜js = −zsqsµsns∇
Our new model for both DG and DD equations with the state variables φ, us and ζs and their associated boundary conditions (BCs) is re-organized as follows: The boundary conditions are changed accordingly to
exp
By the definition of divergence, the quantum-correction drift diffusion equations thus become in the 1D approximation
−1 where the simulation region which includes the channel plus baths runs from 0 to L = lc+2lb and us(x) = us(z) /A (z) , φ (x) = φ (z) , N (x) = N (z) /A (z) , Nbi(x) = Nbi(z) /A (z) , ζs(x) = ζs(z) /A (z) . For the parabolic transport equation (3.14), the elliptic Poisson’s equation (3.13) and density gradient equation (3.15), Dirichlet boundary conditions on φ, us and ζs are imposed at the left and right boundaries.
Note that in this 1D approximation, the current density (and thus the ion velocity) and the electric field will lie in the z direction. To simplify the mathematical analysis of the 1D models (2.12) − (2.13) and (3.13) − (3.15) ,
the cross sectional area A (z) may be set to a constant in the channel itself, but not in the baths.