# A Bubble Model for the Gating of Kv Channels

Zilong Song  Robert Eisenberg  Shixin Xu  Huaxiong Huang  Department of
Mathematics and Statistics, Utah State University, 3900 Old Main Hill, Logan,
UT 84322, USA.  Department of Applied Mathematics, Illinois Institute of
Technology, Chicago, IL, 60616, USA  Department of Physiology and Biophysics,
Rush University, Chicago, IL, 60612, USA.  Duke Kunshan University, 8 Duke
Ave, Kunshan, Jiangsu, China.  Research Centre for Mathematics, Advanced
Institute of Natural Sciences, Beijing Normal University (Zhuhai), China  BNU-
HKBU United International College, Zhuhai, China  Department of Mathematics
and Statistics, York University, Toronto, ON, M3J 1P3, Canada.

###### Abstract

Voltage-gated Kv channels play fundamental roles in many biological processes,
such as the generation of the action potential. The gating mechanism of Kv
channels is characterized experimentally by single-channel recordings and
ensemble properties of the channel currents. In this work, we propose a bubble
model coupled with a Poisson-Nernst-Planck (PNP) system to capture the key
characteristics, particularly the delay in the opening of channels. The
coupled PNP system is solved numerically by a finite-difference method and the
solution is compared with an analytical approximation. We hypothesize that the
stochastic behaviour of the gating phenomenon is due to randomness of the
bubble and channel sizes. The predicted ensemble average of the currents under
various applied voltage across the channels is consistent with experimental
observations, and the Cole-Moore delay is captured by varying the holding
potential.

###### keywords:

voltage-gated channels , Poisson-Nernst-Planck system , bubble model , Cole-
Moore delay

††journal: Journal of Theoretical Biology

##  1 Introduction

Voltage-gated ion channels play fundamental roles in many biological
activities, such as signal generation and propagation in the nervous system,
pacemaker activity in the heart, and coordination of contraction in skeletal
muscle [1, 2, 3]. For example, the voltage-gated Na (Nav) and K (Kv) channels
are key players in the generation of action potential (AP) signals in the
nervous system [4], cardiac and skeletal muscle. This rapid and transient
change of membrane potential propagates long distances (meters) in the nervous
system and muscle fibers as well. The opening and closing of ion channels as
the voltage changes across the membrane determine the depolarization (positive
change in potential) and repolarization (negative change of membrane
potential) that form the propagating AP [5, 6, 7, 8].

The opening of ion channels follows the change in voltage with a delay and
that delay is an important determinant of the conduction velocity. The
conduction velocity helps determine how fast the nervous system can function.
One of the objectives of the present work is to model the delay in the opening
of single Kv channels as well as their ensemble properties. The delay in the
opening of Nav channels is particularly important in determining the
conduction velocity of the action potential. Therefore, understanding the
mechanism of delay is of great biological importance. It is not unreasonable
to expect that the delay is set by a process that is optimized as much as
possible within the constraints of physics, protein structure, and
evolutionary history [9].

Hodgkin and Huxley (HH) provided an empirical model of the generation of AP in
1952 [10]. The conductances they used are ensemble averages of those from many
channels. Understanding the molecular mechanisms that produce these
conductances and the AP is one of the main goals of biophysics for the past
seventy years. Recent advances in structural biology [11] and single-channel
recording [12] have catalyzed our understanding of the physical mechanisms
that produce these conductances. The ionic basis of selective conduction is
now understood reasonably well for sodium channels [13, 14, 15, 16].

The opening and closing of voltage-dependent channels involves many steps [17,
18, 19]. Some of the steps in the voltage-dependent gating of Kv are now known
in molecular and physical detail [20, 21, 22, 23, 24, 25, 26, 27, 28, 4, 29,
30, 31, 32]. The first step is the response of the voltage sensor to the
voltage change, and significant progress has been made in understanding the
physics of that response. It is plausible [10] that the permeability changes
depend on the presence of voltage sensors in the form of charged or dipole
particles, as suggested earlier in a different form [7, 33]. The second step
is the communication of the voltage sensor with the conduction pore of the
channel. This was revealed experimentally in the single-channel ON-OFF
currents (that occur at random intervals) measured by bilayer or patch-clamp
experiments from one channel protein at a time [34, 35]. The development of
patch-clamp experiments [36] was a breakthrough in the understanding of the
gating mechanisms and provided experimental verification at high resolution of
many studies and models.

In the patch-clamp experiments, the recordings of single Kv channels showed a
delay of currents in response to a step voltage change. The ionic current was
generated rapidly after the delay, and vanished when the channel closed
suddenly [29, 30, 4, 37]. The recordings also showed that the delays varied in
each ON-OFF experiment: the gating transitions are stochastic. The ensemble
average has a smoother transient time course for the currents (or opening and
closing of channels), which resembles the classical macroscopic currents (or
voltage-dependent conductances) in the HH model.

The delay in opening was first studied in the inaugural issue of the
Biophysical Journal [38] in the ensemble of channels. Cole and Moore were able
to control the resting potential (i.e., their holding potential) present
before the AP mechanism was turned on. The earlier work of Hodgkin and Huxley
had not addressed this issue in detail because the actual resting potential of
their squid nerve was substantially different from that used as a holding
potential [39]. Hodgkin and Huxley chose to use nerve fibers with more
positive resting potentials so their voltage clamp system could control the
voltage throughout the nerve fiber, something not easy to do [40]. Cole and
Moore found the delay in the response of the nerve fiber to a change in
voltage was much larger when the initial potential (also called the holding
potential) was more negative.

Given the importance of this delay (we call the Cole-Moore delay), it is
striking that a molecular scale biophysical explanation has not been
developed, as far as we know, until very recently [41]. Given the obvious
evolutionary disadvantage of additional delay, it seems likely that whatever
is responsible for the delay is an essential component of the ionic channels
that create the AP. We expect the cause of the Cole-Moore delay to be found in
many channel types where it has not been investigated in detail.

The amount of work on channel proteins that produce the AP has increased
spectacularly in the last decades. The most important single advance (from a
biophysical point of view) is the ON-OFF properties of the single channels,
that in ensemble produce the Cole-Moore delay. Many researchers have proposed
that the ON-OFF property arises from the collapse of a bubble. When the
single-channel current is zero, a region of the voltage sensor with low
effective dielectric constant [26] acts as a hydrophobic gasket that excludes
water and ions from that region of the protein, forming a dewetted region,
which is known as a bubble. Direct evidence for the existence of bubbles is
emerging as structural biologists exploit the magnificent capabilities of
modern techniques of x-ray crystallography and cryoelectron microscopy [20,
42].

Various modeling efforts have been devoted to understanding the gating
mechanisms of Kv channels [22]. In the early years, kinetic models (or called
Markov models) were used for channel gating, by assuming the voltage sensor
has multiple subunits which make transitions between different states [43, 44,
45]. Formally, such kinetic models have some similarity to Hodgkin and
Huxley’s, as the four $n$-gates in the HH model can be interpreted as four
independent subunits that control the gating [4]. Such models have been able
to predict some important features of the gating mechanism of Kv channels
(e.g., Shaker channel) and to fit experimental data, but could not reveal much
about the physics of the gating process. With the availability of more
structural information about channels and advances in computing power,
quantitative models using molecular dynamics (MD) have been developed in
recent decades [46, 47, 2]. MD simulations incorporating physical laws and
interactions of atoms provide insights into the movement of the voltage
sensors, intermediate states, and closure of the pore (forming a dewetted
region). However, the MD approach is limited by the timescale of the
simulations, resolving events in the timescale of $10^{-15}$s, and the total
simulation length is orders of magnitude lower than the timespan (e.g.,
$10^{-3}$ s) of experimentally or biologically relevant processes. This makes
it difficult to directly validate the MD results by using the macroscopic
currents in experiments. To overcome these limitations, alternative multiscale
or macroscopic models [48, 49, 31, 28] have been developed with reasonable
approximations. Some models are based on the formulation of Brownian dynamics,
where the voltage sensor is treated as a Brownian particle [21]. Brownian
models are able to predict macroscopic gating currents, where the free
parameters involved have been estimated based on multiscale modeling
approaches [26].

Here we take a different approach. Following the previous hypothesis of the
hydrophobic region, we construct a specific macroscopic model of a bubble
within the framework of Poisson-Nernst-Planck (PNP) systems and show how it
produces the time course of single Kv channels and the ensemble properties,
including the Cole-Moore delay. The PNP system and its variants have been
found successful in modeling and simulation of many biological processes [50,
51, 28, 52, 53, 54, 55], such as current-voltage curves through ion channels,
the selectivity of ion channels, and ion transport processes in the cell and
tissue scales. In this work, a bubble is assumed to be present in the pore (or
filter) region of the Kv channel, due to the properties of the gating sensor
and channel walls. In the bubble, ions are not present and so cannot carry
charge through it, whereas outside the bubble, ion transport is governed by
the PNP system. The model is constructed so it can easily accommodate more
specific structural information such as the shape, permanent charge (e.g., the
spatial distribution of acid and base residue side chains), and dielectric
properties of the voltage sensor and conduction pore of channels. We calculate
the properties of a single channel containing a bubble and an ensemble average
based on a simple statistical distribution of such channels to represent the
macroscopic currents usually recorded in studies of the opening and closing of
channels. This average does not depend on models [56, 57] of single-channel
kinetics. It only assumes that the opening of each channel (or voltage sensor)
is independent of the others (because channels are many Debye lengths apart,
shielded by the ions, water dipoles (and quadrupoles), and the ionic
atmosphere of proteins and lipid bilayer).

This manuscript is arranged as follows. Section 2 sets up the bubble model
within the framework of PNP systems, followed by a nondimensionalization. In
section 3, the results for a single channel are presented. The bubble model is
solved by a finite-difference method and also solved with analytical
approximations. The results for the profiles of quantities in the model and
the macroscopic currents through the channel are cross-validated by both
methods. Section 4 shows the results for ensemble properties of the Kv
channels and the Cole-Moore delay, with certain assumptions on the statistical
distributions of the bubble locations and cross-sectional area of the channel.
Finally, some concluding remarks are provided in Section 5.

##  2 A bubble model for a voltage-gated Potassium channel

###  2.1 The model setup

We consider a voltage-gated Potassium (Kv) channel in one spatial dimension,
as shown in Figure 1. The total length of channel is set as $2L$, and the
length of the middle (filter and pore) region is $2s$. The positions $x=\pm s$
are the locations of the two edges of the middle region. The bubble, which
carried negative charges with magnitude $q_{b}>0$, can occupy all or part of
it, and is centered at $x=x_{b}$. We assume that the charge is uniformly
distributed inside the bubble. The left chamber is connected to a bath
environment similar to the exterior of a cell, while the right chamber is
connected to one similar to interior of a cell. We further assume that the
right interface of the bubble is fixed at $x=s$ and the left interface
$x=s_{b}=-s+2x_{b}$ is mobile. We anticipate that when the voltage at the
right end of the channel is elevated, the bubble shrinks and moves to the
right. When the left and right interfaces coincide, the bubble vanishes.

Figure 1: Sketch of the Kv channel with the bubble in the middle region.

We consider the case with three ions species K+, Na+ and Cl- (sometimes called
the major bio-ions) outside of the bubble, and the ions can not penetrate into
the bubble. Outside of the bubble, the PNP system is used to model ion
transport

$\displaystyle-\epsilon_{0}\partial_{x}(\epsilon_{r}\partial_{x}\phi)=e_{0}(c_{1}+c_{2}-c_{3}),\quad-L<x<s_{b},~{}s<x<L$
(1) $\displaystyle\frac{\partial c_{i}}{\partial
t}=-\partial_{x}J_{i}=D_{i}\partial_{x}\left(\partial_{x}c_{i}+\frac{e_{0}z_{i}}{k_{B}T}c_{i}\partial_{x}\phi\right),\quad
i=1,2,3,$

where $c_{1}$, $c_{2}$ and $c_{3}$ are the concentrations of K+, Na+ and Cl-
with valences $z_{1}=1$, $z_{2}=1$, $z_{3}=-1$, $\phi$ is the electric
potential, $\epsilon_{r}$ is the dielectric constant, $D_{i}$ ($i=1,2,3$) are
diffusion coefficients, and $\epsilon_{0},e_{0},k_{B},T$ are constants given
in Appendix A.

The boundary conditions are given by

$\displaystyle\phi(-L,t)=0,\quad\phi(L,t)=V_{0}+V_{1}H(t-t_{1}),$ (2)
$\displaystyle c_{i}(-L,t)=c_{i}^{L},\quad c_{i}(L,t)=c_{i}^{R},\quad
i=1,2,3,$

where $V_{0}$ is the initial (holding) membrane potential when the bubble is
in equilibrium (or resting state), $V_{1}$ is the voltage jump at $t=t_{1}$,
$H(t)$ is a Heaviside function, and $c_{i}^{L}$ and $c_{i}^{R}$ ($i=1,2,3$)
are given bath concentrations at the left and right ends [30], which are
electro-neutral. In the experiment, the holding potential $V_{0}$ is not the
same as the Nernst potential of K+. The leak current is allowed to flow
through a different pathway while maintaining $V_{0}$.

Inside the bubble, we have

$\displaystyle-\epsilon_{0}\partial_{x}(\epsilon_{r}\partial_{x}\phi)=-\frac{q_{b}}{V_{b}},\quad
s_{b}(t)<x<s,$ (3)

where $V_{b}=(s-s_{b})A$ is the volume of the bubble and $A$ is the cross
sectional area for the bubble region. The dielectric constant is defined as

$\displaystyle\epsilon_{r}=\begin{cases}\epsilon_{r0},\quad[s_{b},s]\\\
\epsilon_{r1},\quad\mathrm{others}.\end{cases}$ (4)

In addition, we assume that there exists a dipole on the left interface of the
bubble $x=s_{b}$, responsible for maintaining a voltage difference on the two
sides of the bubble. Since the membrane potential is not 0 at equilibrium, the
presence of the dipole with a suitable dipole strength $p_{b}$ guarantees that
the bubble is in equilibrium initially. We can rewrite the equation of $\phi$
in a compact form in the entire domain

$\displaystyle-\epsilon_{0}\partial_{x}(\epsilon_{r}\partial_{x}\phi)=e_{0}(c_{1}+c_{2}-c_{3})-\frac{q_{b}}{V_{b}}+p_{b}\partial_{x}(\delta(x-s_{b})),~{}-L<x<L,$
(5)

with the interpretation that $c_{i}=0$ ($i=1,2,3$) in the bubble and $q_{b}=0$
outside of the bubble.

The total electric force on the bubble is

$\displaystyle\int_{A}\int_{s_{b}}^{s}-\frac{q_{b}}{V_{b}}\left(-\partial_{x}\phi\right)dxdA=q_{b}\frac{\phi(s)-\phi_{(}s_{b})}{s-s_{b}},$
(6)

then the motion of the bubble is modeled by

$\displaystyle\frac{dx_{b}}{dt}=\frac{q_{b}D_{b}}{k_{B}T}\frac{\phi(s)-\phi_{(}s_{b})}{s-s_{b}},$
(7)

where $D_{b}$ ($\ll D_{i}$) is the diffusion coefficient of the bubble. Using
the relationship $s_{b}=-s+2x_{b}$, we can rewrite equation (7) as

$\displaystyle\frac{ds_{b}}{dt}=\frac{2q_{b}D_{b}}{k_{B}T}\frac{\phi(s)-\phi_{(}s_{b})}{s-s_{b}}.$
(8)

At the two interfaces $x=s_{b}$ and $s$, the electric potential and electric
displacement are continuous, and there is no ionic flux across the bubble
interfaces. Mathematically, we set

$\displaystyle[\phi]=0,\quad[\epsilon_{r}\partial_{x}\phi]=0,\quad
J_{i}=0,\;(i=1,2,3),$ (9)

where square brackets mean the jump across the interface, e.g.,
$[\phi(s)]=\phi(s+)-\phi(s-)$. If we include the effect of dipole ($p_{b}$ in
equation (5)) on the interface, we obtain a nonzero jump $[\phi]$ at
$x=s_{b}$. When the two interfaces coincide (i.e., $s_{b}=s$), the bubble
collapses. And we assume that the dipole disappears (i.e., it is treated as an
intrinsic property of the bubble) and the interface conditions are replaced by
continuity conditions

$\displaystyle[\phi]=0,\quad[\epsilon_{r}\partial_{x}\phi]=0,\quad[c_{i}]=0,\quad[J_{i}]=0,\;(i=1,2,3).$
(10)

In summary, we have a system of equations for ion transport coupled with the
motion of the bubble, given by (5), $(\ref{eq1})_{2}$ and (8), together with
boundary and interface conditions (2,9,10). The total current is conserved in
this model, by including three different types of current, given in Appendix
B. This is a special case of the continuity of total current for Maxwell
equations [58, 59], and is also similar to the case of a PNP system for
electric eels [50].

Remark 1. If the dipole does not vanish (i.e., it is treated as property of
the channel or channel wall) after the bubble collapses, we will have nonzero
jump $[\phi]$ related to the dipole, and $[c_{i}]=0$ ($i=1,2,3$) are replaced
by continuity of electro-chemical potentials.

###  2.2 Nondimensionalization

In this subsection, we nondimensionalize our model, which will be used in the
calculations in the subsequent sections. We adopt the following scales

$\displaystyle\tilde{x}=\frac{x}{L},\quad\tilde{s}=\frac{s}{L},\quad\tilde{x}_{b}=\frac{x_{b}}{L},\quad\tilde{V}_{b}=\frac{V_{b}}{LA},$
(11)
$\displaystyle\tilde{\phi}=\frac{\phi}{k_{B}T/e_{0}},\quad\tilde{V}_{0}=\frac{V_{0}}{k_{B}T/e_{0}},\quad\tilde{V}_{1}=\frac{V_{1}}{k_{B}T/e_{0}},$
$\displaystyle\tilde{c}_{i}=\frac{c_{i}}{c_{0}},\quad\tilde{c}_{i}^{L}=\frac{c_{i}^{L}}{c_{0}},\quad\tilde{c}_{i}^{R}=\frac{c_{i}^{R}}{c_{0}},\quad\tilde{D}_{i}=\frac{D_{i}}{D_{0}},(i=1,2,3),$
$\displaystyle\tilde{D}_{b}=\frac{D_{b}}{D_{0}},\quad\tilde{p}_{b}=\frac{p_{b}}{e_{0}c_{0}L^{2}},\quad\tilde{q}_{b}=\frac{q_{b}}{e_{0}},$
$\displaystyle\tilde{t}=\frac{t}{t_{0}},\quad
t_{0}=\frac{L^{2}}{D_{0}},\quad\tilde{J}=\frac{J}{J_{0}},\quad
J_{0}=\frac{D_{0}c_{0}}{L}.$

Some typical values in the above scales and the following boundary conditions
are based on [30] and given in Appendix A.

Substituting (11) into the system in the previous subsection, we obtain a
dimensionless system for variables with tilde (like $\tilde{\phi}$). In order
to simply the notations, we drop the tilde and use the quantities (like
$\phi$) in the dimensionless system. We have the following set of equations in
nondimensional form

$\displaystyle-\epsilon\partial_{x}(\epsilon_{r}\partial_{x}\phi)=c_{1}+c_{2}-c_{3}-\frac{1}{\beta}\frac{{q}_{b}}{(s-s_{b})}+p_{b}\partial_{x}(\delta(x-s_{b})),\quad-1<x<1$
(12) $\displaystyle\frac{\partial c_{i}}{\partial
t}=-\partial_{x}J_{i}=D_{i}\partial_{x}\left(\partial_{x}c_{i}+z_{i}c_{i}\partial_{x}\phi\right),\quad
i=1,2,3,\quad-1<x<s_{b},\,s<x<1,$

with the interpretation that $c_{i}=0$ ($i=1,2,3$) in the bubble
$x\in[s_{b},s]$ and $q_{b}=0$ outside of the bubble. Here the two
dimensionless parameters are defined by

$\displaystyle\epsilon=\frac{\epsilon_{0}k_{B}T}{e_{0}^{2}c_{0}L^{2}},\quad\beta=LAc_{0}.$
(13)

The dielectric constant remains the same

$\displaystyle\epsilon_{r}=\begin{cases}\epsilon_{r0},\quad[s_{b},s]\\\
\epsilon_{r1},\quad\mathrm{others}\end{cases}$ (14)

The motion of the bubble is given by

$\displaystyle\frac{ds_{b}}{dt}=2D_{b}q_{b}\frac{\phi(s)-\phi_{(}s_{b})}{s-s_{b}}.$
(15)

Boundary conditions are given by

$\displaystyle\phi(-1,t)=0,\quad\phi(1,t)=V_{0}+V_{1}*H(t-t_{1}),$ (16)
$\displaystyle c_{i}(-1,t)=c_{i}^{L},\quad c_{i}(1,t)=c_{i}^{R},\;(i=1,2,3).$

Interface conditions are

$\displaystyle[\phi]=0,\quad[\epsilon_{r}\partial_{x}\phi]=0,\quad
J_{i}=0,(i=1,2,3),\quad\textrm{at}~{}~{}x=s_{b},s.$ (17)

After the bubble collapses (for the case that the dipole disappears), we have

$\displaystyle[\phi]=0,\quad[\epsilon_{r}\partial_{x}\phi]=0,\quad[c_{i}]=0,\quad[J_{i}]=0,(i=1,2,3),\quad\textrm{at}~{}~{}x=s_{b}=s.$
(18)

##  3 Results for a single channel

We first compute the initial state when the bubble is in equilibrium by
solving the system of equations with a numerical method, followed by the
results of the non-equilibrium state including the motion of the bubble and
time evolution of the concentrations and electric potential. After the bubble
collapses, the ionic fluxes reaches a steady state. In addition, we also
present the results obtained with an approximate solution (and numerical
evidence) for the intermediate quasi-static states and the final steady state.

###  3.1 Initial state and strength of dipole

We examine the case that the bubble initially occupies the entire middle
region and stays at equilibrium, i.e., $s_{b}=-s$. When $V_{0}=0$, the bubble
is in equilibrium due to symmetry. If $V_{0}\neq 0$, one the other hand,
equilibrium is achieved for an appropriate dipole strength $p_{b}$.

Near the interface $x=s_{b}$, the effect of the other terms is small compared
the dipole, and equation (12) becomes

$\displaystyle-\epsilon\partial_{x}(\epsilon_{r}\partial_{x}\phi)=p_{b}\partial_{x}(\delta(x-s_{b})),$
(19)

and integrating once gives

$\displaystyle-\epsilon\epsilon_{r}\partial_{x}\phi(x)=p_{b}\delta(x-s_{b})+C.$
(20)

By integrating again and taking the limit of $x\rightarrow s_{b}$, we obtain

$\displaystyle[\phi(s_{b})]=-\frac{p_{b}}{\epsilon}\left(\frac{1}{2\epsilon_{r0}}+\frac{1}{2\epsilon_{r1}}\right).$
(21)

Therefore, for a given $V_{0}$, we find the following relationship

$\displaystyle-\frac{p_{b}}{\epsilon}\left(\frac{1}{2\epsilon_{r0}}+\frac{1}{2\epsilon_{r1}}\right)=V_{0},$
(22)

and then the bubble will be in equilibrium as in the symmetric case with
$V_{0}=0$.

For the equilibrium profile, the fluxes are 0 and one can not distinguish the
effects of the two positive ions Na+ and K+. We can group the two positive ion
species and treat them as a single species. The boundary conditions for
$c_{1}+c_{2}$ and $c_{3}$ will be the same, and hence we will have exact
symmetry for this equilibrium case. The equilibrium profiles can be determined
analytically, and we take $V_{0}=0$ and $p_{b}=0$ in the derivation. Since the
bubble is in equilibrium, inside the bubble we have (note $s_{b}=-s$)

$\displaystyle\phi(x)=B_{1}x^{2}+\phi(0),\quad
B_{1}=\frac{{q}_{b}}{4s\epsilon\epsilon_{r0}\beta}.$ (23)

Taking the derivative and together with interface conditions at $x=s$, we have

$\displaystyle\epsilon_{r1}\phi^{\prime}(s+)=\epsilon_{r0}\phi^{\prime}(s-)=\epsilon_{r0}2B_{1}s=\frac{q_{b}}{2\epsilon\beta}.$
(24)

Due to symmetry, we only consider the right chamber $s<x<1$. It is easy to
verify that the PNP system (12) in equilibrium reduces to

$\displaystyle\epsilon\epsilon_{r1}\phi^{\prime\prime}=c_{3}-(c_{1}+c_{2})=e^{\phi}-e^{-\phi},$
(25)

where $c_{3}^{R}=1$ has been used. Integrating once gives

$\displaystyle\frac{1}{2}\epsilon\epsilon_{r1}[(\phi^{\prime}(x))^{2}-(\phi^{\prime}(s+))^{2}]=e^{\phi}+e^{-\phi}-(e^{\phi_{s}}+e^{-\phi_{s}}),$
(26)

where $\phi_{s}=\phi(s)$. Then, by combining with (23), we get

$\displaystyle(\phi^{\prime}(x))^{2}=G(\phi)$
$\displaystyle=\left(\frac{q_{b}}{2\epsilon\beta\epsilon_{r1}}\right)^{2}+\frac{2}{\epsilon\epsilon_{r1}}\left(e^{\phi}+e^{-\phi}-(e^{\phi_{s}}+e^{-\phi_{s}})\right),$
(27)

which leads to the solution

$\displaystyle x=\int_{\phi_{s}}^{\phi}\frac{1}{\sqrt{G(\phi)}}d\phi+s.$ (28)

The unknown constant $\phi_{s}$ in the solution can be determined by the
condition

$\displaystyle 1=\int_{\phi_{s}}^{0}\frac{1}{\sqrt{G(\phi)}}d\phi+s.$ (29)

Remark 2. Because of symmetry, we can estimate $\phi_{s}$ from the above
derivation as

$\displaystyle\phi_{s}\approx-\ln\left(\frac{q_{b}^{2}}{8\epsilon\epsilon_{r1}\beta^{2}}\right)$
(30)

for $q_{b}$ in a certain range.

Figure 2: Equilibrium electric potential $\phi$ and concentrations
$c_{1},c_{2},c_{3}$ when $q_{b}=2$.

For the numerical results obtained in this paper, we vary the potentials at
the two end of the domain while fixing the magnitude of permanent charge as
$q_{b}=2$. Most of the other parameter values used for the computation are
also fixed and given in Appendix A.

In Figure 2(a), the electric potential $\phi$ is plotted for the case of
$V_{0}=0$. When $V_{0}=-3.18$ (i.e., -80 mV), we obtain $p_{b}\approx 0.044$
using (22) and the electric potential is plotted in Figure 2(b), where the
jump at the interface $x=s_{b}$ is due to the presence of the dipole. The
concentrations $c_{1}$, $c_{2}$ and $c_{3}$ (which can be computed from
solution of $\phi$) are shown in Figures 2(c-e). It can be seen that
$c_{1}+c_{2}$ is symmetric (Figure 2(d)) as expected. The initial membrane
potential $V_{0}$ is balanced by the jump of $\phi$ due to the presence of the
dipole. In the non-equilibrium case (before the bubble collapses), we ignore
both the initial membrane potential $V_{0}$ and the dipole, so that the value
of $\phi$ is continuous at the interfaces. The solutions in Figure 2 will be
verified by numerical simulations in the subsequent subsections.

###  3.2 The dynamics of the bubble motion and channel currents

In this part, we present numerical solutions of the PNP system and bubble
motion. Inside the bubble, there exist no ions and their concentrations
$c_{i}$ ($i=1,2,3$) are zero. For convenience, the PNP system is solved inside
the bubble by assigning small diffusion coefficients ($D_{1}=D_{2}=10^{-15}$).

Figure 3: Equilibrium electric potential $\phi$, concentrations $c_{1}$,
$c_{2}$, $c_{1}+c_{2}$, $c_{3}$, and ionic fluxes $J_{i}$ ($i=1,2,3$) with
$h=0.0025$.

The finite difference method is used to solve the system, with a uniform mesh
$h=x_{k}-x_{k-1}$. A temporal semi-implicit discrete scheme is used with
$t_{n}=n\Delta t$ and $x_{k}=x_{0}+kh$, given by

$\displaystyle-\epsilon\frac{\epsilon_{r,k-1/2}}{h^{2}}\phi_{k-1}^{n+1}+\epsilon\frac{\epsilon_{r,k-1/2}+\epsilon_{r,k+1/2}}{h^{2}}\phi_{k}^{n+1}-\epsilon\frac{\epsilon_{r,k+1/2}}{h^{2}}\phi_{k+1}^{n+1}$
(31) $\displaystyle\hskip 142.26378pt-
c_{1,k}^{n+1}-c_{2,k}^{n+1}+c_{3,k}^{n+1}=-q_{k}^{n+1},$
$\displaystyle\frac{c_{i,k}^{n+1}-c_{i,k}^{n}}{\Delta
t}=-\frac{J_{i,k+1/2}^{n+1}-J_{i,k-1/2}^{n+1}}{h},\quad i=1,2,3,$
$\displaystyle
J_{i,k+1/2}^{n+1}=-D_{i,k+1/2}\frac{c_{i,k+1}^{n+1}-c_{i,k}^{n+1}}{h}-D_{i,k+1/2}z_{i}c_{k+1/2}^{n}\frac{\phi_{k+1}^{n+1}-\phi_{k}^{n+1}}{h},$

where harmonic average is used for the diffusion coefficient

$\displaystyle
D_{i,k+1/2}=\frac{2}{\frac{1}{D_{i,k}}+\frac{1}{D_{i,k+1}}},\quad i=1,2,3.$
(32)

In this way, we ensure that the ionic fluxes are small near the interface as
approximations of $J_{i}=0$ ($i=1,2,3$). When the bubble collapses, the
diffusion coefficient is guaranteed to be the same as that outside of the
bubble, and the continuity conditions are recovered. The quantities
$\epsilon_{r,k+1/2}$ and $q_{k}^{n}$ in (31) are defined in Appendix B. The
discrete scheme also preserves the continuity of the total current, as in the
original continuous model.

Figure 4: Electric potential $\phi$ and ionic concentrations
$c_{1},c_{2},c_{3}$ at three different times with $V_{1}=6.36$.

For $q_{b}=2$, we first compute the initial equilibrium when the bubble
occupies the entire middle region. The initial condition at $t=0$ is set as

$\displaystyle\phi(x,0)=0,\quad-1<x<1,$ (33) $\displaystyle
c_{i}(x,0)=\begin{cases}&c_{i}^{L}\quad-1<x<s_{b}=-s,\\\ &0,\quad
s_{b}<x<s,\\\ &c_{i}^{R},\quad s<x<1,\end{cases}$

where $i=1,2,3$. We also set $V_{0}=0$ and $p_{b}=0$ in the computation so
that $\phi$ is continuous. The computation is carried out until the system
reaches a steady state. For a given mesh size $h=0.0025$, Figure 3 shows the
numerical solution of electric potential $\phi$, concentrations $c_{i}$ and
ionic fluxes $J_{i}$ ($i=1,2,3$), which are in good agreement with the
analytical results in the previous subsection.

Figure 5: The fluxes $J_{i}$ ($i=1,2,3$) at three different times with
$V_{1}=6.36$.

Figure 6: The dynamics of ionic fluxes $J_{i}$ ($i=1,2,3$) with $V_{1}=6.36$.

Figure 7: The steady state with $V_{1}=6.36$.

Next, we present the results on the bubble motion and dynamic behaviour of the
PNP system. We start from the equilibrium state at $t=0$ and increase the
electric potential from zero to $\phi(1,t)=V_{1}=6.36$ (i.e., 160 mV) for
$0<t<t^{\ast}$ (we already set $V_{0}=0$ and $p_{b}=0$), where $t^{\ast}$ is
the unknown time when the bubble collapses. Figure 4 shows $\phi,c_{i}$
($i=1,2,3$) at three different times. The minimum value for $\phi$ inside the
bubble gradually increases in Figure 4(a), and the interface $s_{b}$ moves to
the right as indicated by Figure 4. Figure 5 shows the ionic fluxes at three
different times, which are small. After the bubble collapses and the dipole
disappears, we reset $\phi(1,t)=V_{0}+V_{1}=3.18$. Figure 6 shows the three
ionic fluxes at $x=\pm 1$. It can be seen that they remain small until the
bubble collapses (i.e., $s_{b}=s$) at $t^{\ast}=3.13\times 10^{6}$, which is
17.6 ms in dimensional unit. The ionic fluxes and $\phi,c_{i}$ ($i=1,2,3$)
reach a steady state soon after the bubble collapses, as shown in Figure 7. At
steady state, the dimensionless ionic flux $J_{1}$ and the dimensional current
$I$ are found to be

$\displaystyle J_{1}\approx-2.834,\quad I=|J_{1}|e_{0}AJ_{0}\approx
10\,\mathrm{pA}.$ (34)

Remark 3. The value of the steady state current $I$ obtained above is close to
that given in Figure 2(a) of [30]. When the voltage jump $V_{1}$ is reset to
zero after system reaches a steady state, the ionic fluxes reduce to zero
immediately, indicating the closure of the ion channel. In this sense, our
proposed model provides a plausible gating mechanism once the bubble is
generated. However, the mechanism of the bubble generation is not considered
here and will be the subject of a future study.

###  3.3 Quasi-static equilibrium

Since the motion of the bubble is extremely slow compared with the diffusive
timescale of the ions, ionic fluxes are essentially zero (Figure 5) before the
bubble collapses. Therefore, we can use quasi-static solution with zero ionic
fluxes as an approximation of the intermediate states. A hybrid method can be
used to determine the solution of intermediate states by first obtaining an
analytical solution (in terms of integrals), where integration constants
involved can be determined easily using a numerical method afterwards.

Given boundary condition $V_{1}$ and interface position $s_{b}$, solving the
quasi-static equilibrium is similar to that for solving the initial state. We
set $V_{0}=0$ and $p_{b}=0$ so that the continuity condition of $\phi$ can be
used at interface $s_{b}$. Inside the bubble, we have

$\displaystyle\phi(x)=B_{1}(x-s)^{2}+\tilde{\phi}_{s}(x-s)+\phi_{s},\quad
B_{1}=\frac{{q}_{b}}{2(s-s_{b})\epsilon\epsilon_{r0}\beta},$ (35)

where $\phi_{s},\tilde{\phi}_{s}$ are to be determined. The solutions of
$\phi$ outside of the bubble can be written as

$\displaystyle
x=\int_{\phi_{s}}^{\phi}\frac{1}{\sqrt{G_{1}(\phi;\phi_{s},\tilde{\phi}_{s})}}d\phi+s,\quad
s<x<1,$ (36) $\displaystyle
x=-\int_{\phi_{s_{b}}}^{\phi}\frac{1}{\sqrt{G_{2}(\phi;\phi_{s},\tilde{\phi}_{s})}}d\phi+s_{b},\quad-1<x<s_{b},$

where $\phi_{s_{b}}=\phi(s_{b})$ can be expressed by $\phi_{s}$ and
$\tilde{\phi}_{s}$. The derivation for $G_{1}(\phi),G_{2}(\phi)$ are given in
Appendix C. For given parameter values including $V_{1}$ and $s_{b}$, the two
unknowns $\phi_{s}$ and $\tilde{\phi}_{s}$ can be determined by the two
boundary conditions $\phi(1)=V_{1}$ and $\phi(-1)=0$, i.e.,

$\displaystyle
1=\int_{\phi_{s}}^{V_{1}}\frac{1}{\sqrt{G_{1}(\phi;\phi_{s},\tilde{\phi}_{s})}}d\phi+s,\quad-1=-\int_{\phi_{s_{b}}}^{0}\frac{1}{\sqrt{G_{2}(\phi;\phi_{s},\tilde{\phi}_{s})}}d\phi+s_{b},$
(37)

and $\phi(x),c_{i}(x)$ ($i=1,2,3$) can be obtained afterwards. As an example,
for $q_{b}=2$, $V_{1}=6.36$ and $s_{b}=0$, solutions of $\phi$ and $c_{i}$
($i=1,2,3$) can be computed using the procedure outlined above and plotted in
Figure 8.

For fixed $V_{1}$, we can treat $\phi(s),\phi(s_{b})$ as functions of the
parameter $s_{b}$, which can be determined by

$\displaystyle\frac{ds_{b}}{dt}=2D_{b}q_{b}\frac{\phi(s)-\phi_{(}s_{b})}{s-s_{b}}=2D_{b}q_{b}f(s_{b}).$
(38)

Integrating in time, we obtain $t^{\ast}$, the time delay after the voltage
jump and before the bubble collapses,

$\displaystyle t^{\ast}=\int_{-s}^{s}\frac{1}{2D_{b}q_{b}f(x)}dx.$ (39)

Figure 9(a) shows the dependence of quantities $\phi(s)$ and $\phi(s_{b})$ on
$s_{b}$ and Figure 9(b) shows the function $f(s_{b})$. From (39), we find that
$t^{\ast}\approx 3.26\times 10^{6}$, which is $18.3$ ms in dimensional unit,
which is slightly longer than that obtained using the finite difference method
(17.6 ms) previously.

Figure 8: The quasi-static solution of $\phi$, $c_{1}$, $c_{2}$, $c_{3}$ and
$c_{1}+c_{2}$ for $s_{b}=0$ and $V_{1}=6.36$.

Figure 9: The dependence of $\phi(s)$ and $\phi(s_{b})$ on $s_{b}$ and the
function $f(s_{b})$ for $V_{1}=6.36$.

###  3.4 The steady state after the collapse of the bubble

After the bubble collapses, interface conditions $J_{i}=0$ ($i=1,2,3$) are
replaced by continuity conditions $[J_{i}]=0$ and $q_{b}$ becomes a point
charge (a delta function). Due to the presence of $q_{b}$, the concentration
of $c_{3}$ is approximately zero near $x=s$, and we assume $J_{3}/D_{3}\approx
0$. The system at the steady state can be approximated by

$\displaystyle-\epsilon\epsilon_{r1}\phi^{\prime\prime}(x)$
$\displaystyle=c_{1}+c_{2}-c_{3}-\frac{q_{b}}{\beta}\delta(x-s),$ (40)
$\displaystyle-J_{1}$ $\displaystyle=c_{1}^{\prime}+c_{1}\phi^{\prime},$
$\displaystyle-\frac{J_{2}}{D_{2}}$
$\displaystyle=c_{2}^{\prime}+c_{2}\phi^{\prime},$ $\displaystyle 0$
$\displaystyle=c_{3}^{\prime}-c_{3}\phi^{\prime},$

where $D_{1}=1$ has been used. If we combine the effects of $c_{1},c_{2}$ and
define $J_{p}=J_{1}+\frac{J_{2}}{D_{2}}$, then the system can be reduced to a
single equation of $\phi$ (see the derivation in Appendix C)

$\displaystyle-\epsilon\epsilon_{r1}\phi^{\prime\prime}(x)=\frac{1}{2}\epsilon\epsilon_{r1}\left[(\phi^{\prime}(x))^{2}-(\phi^{\prime}(1))^{2}\right]-J_{p}(x-1)-2(e^{\phi-V}-1),$
(41) $\displaystyle\hskip 227.62204pt\mathrm{for}\quad s<x<1,$
$\displaystyle-\epsilon\epsilon_{r1}\phi^{\prime\prime}(x)=\frac{1}{2}\epsilon\epsilon_{r1}\left[(\phi^{\prime}(x))^{2}-(\phi^{\prime}(-1))^{2}\right]-J_{p}(x+1)-2(e^{\phi}-1),$
$\displaystyle\hskip 227.62204pt\mathrm{for}\quad-1<x<s,$

where $V=V_{0}+V_{1}$ and $c_{3}^{L}=c_{3}^{R}=1$ have been used. The point
charge (delta function) at $x=s$ gives the jump condition

$\displaystyle\,[\phi^{\prime}]=\frac{q_{b}}{\beta\epsilon\epsilon_{r1}}.$
(42)

Figure 10: Comparison of $\phi$ and $c_{i}$ ($i=1,2,3$) at steady state with
$V_{1}=6.36$ between semi-analytical approximation (solid lines) and the
finite-difference solution (dashed lines).

Given $\phi^{\prime}(1)$, $\phi^{\prime}(-1)$, and $J_{p}$, the solutions can
be easily determined numerically, in the two regions $s<x<1$ and $-1<x<s$. The
three constants $\phi^{\prime}(1),\phi^{\prime}(-1),J_{p}$ can be determined
by condition (42), $[\phi]=0$ and $[c_{1}+c_{2}]=0$ at $x=s$ (in practice the
numerical procedure is more stable if the ratio
$(c_{1}(s+)+c_{2}(s+))/(c_{1}(s-)+c_{2}(s-))=1$ is used instead of
$[c_{1}+c_{2}]=0$). Once $\phi(x)$ is obtained, $c_{1}$ and flux $J_{1}$ can
be computed by equation $(\ref{eq40})_{2}$ and the continuity condition
$[c_{1}]=0$ at $x=s$. Similarly, $c_{2}$ and flux $J_{2}$ can be computed by
equation $(\ref{eq40})_{3}$ and the continuity condition $[c_{2}]=0$ at $x=s$.

Figure 11: Steady state $\phi$ and $c_{i}$ ($i=1,2,3$) for $V_{1}=1.59$.

Figure 12: The dependence of the dimensionless flux $J_{1}$ and the
dimensional current $I$ on the voltage $V=V_{0}+V_{1}$ at steady state.

Figure 10 shows the semi-analytical approximation $\phi$ and $c_{i}$
($i=1,2,3$) at steady state for $V_{1}=6.36$, which agrees with those in
Figure 7 except for $c_{3}$. For $V_{1}=6.36$, ionic flux $J_{1}$ is found to
be $J_{1}\approx-2.855$ (also $\approx$ 10 pA in units), which is close to
$-2.834$ in (34) obtained by the finite difference method. The advantage of
the method in this subsection is that the computation is extremely fast
compared with the full finite difference method. It is much more efficient to
use the semi-analytical approximation to compute the steady states
(particularly the currents) with various different voltage jump $V_{1}$.

Figure 11 shows the results with $V_{1}=1.59$ (i.e., 40 mV in physical units),
and the flux is $J_{1}\approx-0.264$ (i.e., 0.933 pA in units). In above
computations, $J_{2}$ is very small since a small $D_{2}$ is used, therefore
the dimensionless total current is almost the same as $J_{1}$. Figure 12 shows
the dependence of dimensionless flux $-J_{1}$ and the dimensional current on
the voltage $V=V_{0}+V_{1}$ at steady state.

Remark 4. In the present simple model, the permanent charge is evenly
distributed in the bubble only and the size effect of different ions are not
considered. We do not expect our model to capture the current-voltage relation
for large $V_{1}$, including the saturation phenomenon observed experimentally
in the literature. To make our model more realistic, we need to know the
distribution of permanent charge (i.e., acid base side chains) along the
system. When that information is available, it can be incorporated into our
model by adding permanent charge to the channel wall in the region $-s$ to
$+s$ (see Fig. 1) as in a practical implementation [60]. The studies [54, 55]
by Weishi Liu and his group have illustrated the effects of permanent charge
on current-voltage relation. With ionic size effect and the permanent charge,
saturation phenomenon of current-voltage curves can be modelled as shown in
[51].

Remark 5. For the case that the dipole does not disappear after bubble
collapses, Figure 13 shows the results for steady state flux, which are quite
similar to those in Figure 12 for the above case when dipole disappears after
bubble collapses.

Figure 13: The dependence of the dimensional flux $J_{1}$ and the dimensional
current $I$ on the voltage $V=V_{0}+V_{1}$ at steady state for the case that
the dipole does not disappear after the bubble collapses.

##  4 Ensemble properties

In this section, we extend our model by including stochastic effect in two
aspects. We assume that the initial position of the bubble and the cross
sectional area of the channel are both random and compute the ensemble
averages of the macroscopic currents through the channel and estimate the
Cole-Moore delay based on certain statistical distributions.

First, we assume that the initial interface position $s_{b}$ is random, which
could be due to the tiny fluctuations of strength of dipoles in different
channels or due to the mechanism of bubble formation (which is not considered
in the present work). For illustration, we consider that $s_{b}$ follows a
normal distribution

$\displaystyle s_{b}\sim N(\mu,\sigma^{2}),\quad\mu=0,\quad\sigma=0.05,$ (43)

where the choice of $\sigma$ ensures that $s_{b}\in[-s,s]$ with $s=0.2$ for
almost all the generated data. We can use the previous function $f(s_{b})$ to
compute the ensemble properties of the channel, since there is negligible
effect on the curves of $f(s_{b})$ with different starting value of $s_{b}$.
With each different initial position $s_{b}$, the dynamics of the fluxes
(particularly the time delay $t^{\ast}$ for opening of the channel) will be
different. By taking the average of these fluxes, we get the ensemble curve
for the dynamics of the current through the channel (i.e., fluxes of K+).
Figure 14 (b,c) show the ensemble curves for the current $I$ and the ratio
$I/V_{1}$ with 50 channels and with 4 different voltage jumps, which are given
in Figure 14(a). Figure 14 (b) shows similar trend and scale with experiments
in Figure 2(a,c) in [30].

Figure 14: The ensemble curves for the current $I$ and the ratio $I/V_{1}$
averaged by 50 random $s_{b}$ and with 4 different $V_{1}$. Figure 15: The
function $f(s_{b})$ with 3 values of the parameter $A_{1}$.

Figure 16: The dependence of the flux $-J_{1}$ and the current $I$ on the
parameter $A_{1}$. Figure 17: The ensemble curves for the current with 50
channels, where $s_{b}\sim N(0,\sigma^{2})$ and $A_{1}(t_{i})\sim
N(1,\sigma_{A}^{2})$, with $\sigma=0.05,\sigma_{A}=0.03$.

The case is more complicated when cross sectional area $A$ is random. We set
$A=A_{0}A_{1}$, where $A_{1}\sim N(1,\sigma_{A}^{2})$ with $\sigma_{A}=0.03$
and $A_{0}=(0.7\mathrm{nm})^{2}$ (the same as the value in Appendix A). The
area $A$ will affect the dimensionless parameter $\beta$, and hence influences
the effective permanent charge $q_{b}/\beta$.

We start by examining the effect of $A_{1}$ on $t^{\ast}$. Figure 15 shows
$f(s_{b})$ with 3 different values of $A_{1}$, indicating that the effect of
$A_{1}$ on $f(s_{b})$ and hence on $t^{\ast}$ is very small. Therefore, the
previous curve $f(s_{b})$ can be used to compute $t^{\ast}$ as an
approximation. We study the effect of $A_{1}$ on the flux $J_{1}$ or the
current $I$ at steady state. Figure 16(a) shows the dependence of $-J_{1}$ on
$A_{1}$, indicating that the magnitude of $J_{1}$ will slightly decrease with
increase of $A_{1}$. Since the final dimensional current also depends on the
scaling factor which contains $A_{1}$, Figure 16(b) shows the dependence of
the current $I$ on the parameter $A_{1}$, indicating that the current
increases with $A_{1}$. Figure 16(b) also shows the approximate current
$A_{1}I(1)$ where $I(1)$ is taken from previous computation with $A_{1}=1$,
which is close to the exact curve. Therefore, the main effect of $A_{1}$ on
the current is due to the scaling factor. We conclude that $A_{1}I(1)$ can be
used as an approximation for the current in the following figures.

We fix $A_{1}$ for each channel during the evolution of the bubble, while
allowing it (together with $s_{b}$) to vary randomly among 50 channels. The
ensemble curves for the current with 50 channels are very similar to those in
Figure 14. We also consider the case that $A_{1}$ fluctuates randomly when the
bubble evolves. We take $A_{1}(t_{i})\sim N(1,\sigma_{A}^{2})$ with
$\sigma_{A}=0.03$ for each discrete time $t=t_{i}$ and for each channel.
Figure 17 shows the ensemble curve for the current with 50 channels, where
$s_{b}\sim N(0,\sigma^{2})$ and $A_{1}(t_{i})\sim N(1,\sigma_{A}^{2})$, with
$\sigma=0.05,\sigma_{A}=0.03$, and 400 discrete $t_{i}$ are used for the time
interval of 30 ms. Figure 17 shows similar trend and fluctuations with those
in experiments, see Figure 2(a,c) in [30] and Figure 3.17 in [4].

To model the Cole-Moore delay [38], we can treat the mean value $\mu$ and
standard variation $\sigma$ in (43) as a function of the holding potential
$V_{0}$. For illustration, we take

$\displaystyle s_{b}\sim N(\mu,\sigma^{2}),$ (44)
$\displaystyle\mu(V_{0})=s\tanh(k(V_{0}-V_{0}^{\ast})),\quad
V_{0}^{\ast}=-80\mathrm{mV},\quad k=0.002/\mathrm{mV},$
$\displaystyle\sigma(V_{0})=\sigma_{0}=0.05$

where $V_{0}$ is the initial holding potential, and $V_{0}^{\ast}$ is a
reference value. Figure 18(a) shows the ensemble curves for the current with
100 channels and $V_{1}=160$ mV, for 7 different holding potential $V_{0}$
which are [-52, -72, -93, -113, -133, -162, -212]mV, corresponding to curves
from left to right. The ensemble curves show similar features as experimental
curves in Figure 18(b), which is reproduced from Figure 5(a) in [38]. Figure
19 shows the ensemble curves for the current with 100 channels and 600
channels and with $V_{1}=160$ mV, for 2 different holding potential
$V_{0}=-52,-212$ mV. It can be observed from the figures that the delay is
longer when holding potential $V_{0}$ is smaller.

Remark 6. We note that the Cole-Moore effect may also arise in the hydrophobic
gasket of the voltage sensor region of the channel and show itself as a delay
in gating current [41]. The bubble in the voltage sensor itself would not
collapse, and the gating current would be given by our equations (50-52). We
speculate that some of the gating current could flow in the adjacent
conduction pore, and open it, perhaps by collapsing a bubble in the conduction
pore.

Figure 18: (a) The ensemble curves for the current with 100 channels and
$V_{1}=160$ mV, for 7 different holding potentials $V_{0}$, (b) the ensemble
curve reproduced based on the experiments in Figure 5(a) of [38].

Figure 19: The ensemble curves for the current for $V_{1}=160$ mV and 2
different holding potential $V_{0}$, (a) with 100 channels, and (b) with 600
channels.

##  5 Conclusion

In this paper, we present a macroscopic bubble model for the gating of Kv
Channels. The time delay in the opening of a single channel is determined by
the motion of the bubble before it collapses. The bubble motion is coupled
with a Poisson-Nernst-Planck system, which is solved by a full numerical
computation as well as a quasi-static approximation method. We also present a
stochastic model for the bubble and channel sizes and the ensemble properties
of the Kv channel are consistent with experimental observations. Furthermore,
the Cole-Moore delay is explored by assuming the dependence of bubble
properties on the holding potential.

Although the present simple model captures some key features in the ensemble
properties, some parts are oversimplified and there is room for improvement.
The permanent charges in the channel are lumped together in the model, and the
distinction and effects of charges on bubble and on the channel wall could be
examined in the future. The selectivity of channel is not considered in detail
here, which depends on the ion sizes (which makes the PNP system very
complicated). This is circumvented by assuming small diffusion constants of
other ions except K+ in the present work. The generalization to high-
dimensional case is also interesting and nontrivial, since the bubble
interface will have a curved shape and specific forces (e.g., some force due
to the maxwell stress) can act on the interface.

## References

  * [1] M. Biel, C. Wahl-Schott, S. Michalakis, X. Zong, Hyperpolarization-activated cation channels: from genes to function, Physiological reviews 89 (3) (2009) 847–885.
  * [2] M. Ø. Jensen, V. Jogini, D. W. Borhani, A. E. Leffler, R. O. Dror, D. E. Shaw, Mechanism of voltage gating in potassium channels, Science 336 (6078) (2012) 229–233.
  * [3] D. A. Jacobson, F. Mendez, M. Thompson, J. Torres, O. Cochet, L. H. Philipson, Calcium-activated and voltage-gated potassium channels of the pancreatic islet impart distinct and complementary roles during secretagogue induced electrical responses, The Journal of physiology 588 (18) (2010) 3525–3537.
  * [4] B. Hille, Ion channels of excitable membranes, Sinauer Associates, Inc., 2001.
  * [5] A. Huxley, The quantitative analysis of excitation and conduction in nerve. from nobel lectures, physiology or medicine 1963–1970 (1972).
  * [6] A. F. Huxley, Kenneth Stewart Cole 1900-1984. A biographical memoir by Sir Andrew Huxley, National Academies Press, Washington DC, 1996.
  * [7] A. Huxley, From overshoot to voltage clamp, Trends in Neurosciences 25 (11) (2002) 553–558.
  * [8] A. Hodgkin, Chance and Design, Cambridge University Press, New York, 1992.
  * [9] S. J. Gould, The Flamingo’s Smile: Reflections in Natural History, Norton, New York, 1987.
  * [10] A. L. Hodgkin, A. F. Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve, The Journal of physiology 117 (4) (1952) 500.
  * [11] R. MacKinnon, Nobel lecture. potassium channels and the atomic basis of selective ion conduction, Biosci Rep 24 (2) 75–100.
  * [12] E. Neher, Ion channels for communication between and within cells Nobel Lecture, December 9, 1991, World Scientific Publishing Co, Singapore, 1997, pp. 10–25.
  * [13] D. Boda, W. Nonner, M. Valisko, D. Henderson, B. Eisenberg, D. Gillespie, Steric selectivity in na channels arising from protein polarization and mobile side chains, Biophys J 93 (6) (2007) 1960–80.  doi:10.1529/biophysj.107.105478.
URL http://www.ncbi.nlm.nih.gov/pubmed/17526571

  * [14] T. Dudev, C. Lim, Factors governing the na+ vs k+ selectivity in sodium ion channels, Journal of the American Chemical Society 132 (7) (2010) 2321–2332.
  * [15] T. Dudev, C. Lim, Ion selectivity strategies of sodium channel selectivity filters, Accounts of chemical research 47 (12) (2014) 3580–3587.
  * [16] C. Lim, T. Dudev, Potassium versus sodium selectivity in monovalent ion channel selectivity filters, Springer, 2016, pp. 325–347.
  * [17] F. Bezanilla, How membrane proteins sense voltage, Nat Rev Mol Cell Biol 9 (4) 323–32.  doi:nrm2376[pii]10.1038/nrm2376.
  * [18] F. Bezanilla, Ion channels: from conductance to structure, Neuron 60 (3) 456–68.  doi:S0896-6273(08)00900-8[pii]10.1016/j.neuron.2008.10.035.
  * [19] F. Bezanilla, Gating currents, The Journal of General Physiology 150 (7) (2018) 911–932.  doi:10.1085/jgp.201812090.
URL https://doi.org/10.1085/jgp.201812090

  * [20] J. J. Lacroix, H. C. Hyde, F. V. Campos, F. Bezanilla, Moving gating charges through the gating pore in a kv channel voltage sensor, Proceedings of the National Academy of Sciences 111 (19) (2014) E1950–E1959.
  * [21] L. Catacuzzeno, F. Franciolini, Simulation of gating currents of the shaker k channel using a brownian model of the voltage sensor, Biophysical journal 117 (10) (2019) 2005–2019.
  * [22] L. Catacuzzeno, L. Sforna, F. Franciolini, Voltage-dependent gating in k channels: experimental results and quantitative models, Pflügers Archiv-European Journal of Physiology 472 (1) (2020) 27–47.
  * [23] L. Catacuzzeno, L. Sforna, F. Franciolini, R. Eisenberg, Why are voltage gated na channels faster than k channels? one multi-scale hierarchical model, bioRxiv. Cold Spring Harbor Laboratory 11.
  * [24] C. A. Bassetto, J. L. Carvalho-de Souza, F. Bezanilla, Molecular basis for functional connectivity between the voltage sensor and the selectivity filter gate in shaker k+ channels, Elife 10 (2021) e63077.
  * [25] L. Catacuzzeno, F. Franciolini, F. Bezanilla, R. S. Eisenberg, Gating current noise produced by brownian models of a voltage sensor, Biophysical Journal 120 (September 21, 2021) (2021) 1–19.  doi:10.1016/j.bpj.2021.08.015.
URL https://dx.doi.org/10.1016/j.bpj.2021.08.015

  * [26] L. Catacuzzeno, L. Sforna, F. Franciolini, R. S. Eisenberg, Multiscale modeling shows that dielectric differences make nav channels faster than kv channels, Journal of General Physiology 153 (2).
  * [27] T.-L. Horng, R. S. Eisenberg, C. Liu, F. Bezanilla, Gating current models computed with consistent interactions, Biophysical Journal 110 (3) (2016) 102a–103a.
  * [28] T.-L. Horng, R. S. Eisenberg, C. Liu, F. Bezanilla, Continuum gating current models computed with consistent interactions, Biophysical journal 116 (2) (2019) 270–282.
  * [29] I. Llano, R. J. Bookman, Ionic conductances of squid giant fiber lobe neurons., The Journal of general physiology 88 (4) (1986) 543–569.
  * [30] I. Llano, C. K. Webb, F. Bezanilla, Potassium conductance of the squid giant axon. single-channel studies., The Journal of general physiology 92 (2) (1988) 179–196.
  * [31] I. Kim, A. Warshel, Coarse-grained simulations of the gating current in the voltage-activated kv1. 2 channel, Proceedings of the National Academy of Sciences 111 (6) (2014) 2128–2133.
  * [32] F. Bezanilla, Voltage sensor movements, The Journal of general physiology 120 (4) (2002) 465–473.
  * [33] A. Hodgkin, A. Huxley, B. Katz, Ionic currents underlying activity in the giant axon of the squid, Arch. Sci. physiol. 3 (1949) 129–150.
  * [34] B. Sakmann, E. Neher, Single Channel Recording., 2nd Edition, Plenum, New York, 1995\.
  * [35] J. Zheng, M. C. Trudeau, Handbook of ion channels, CRC Press, 2015.
  * [36] O. P. Hamill, A. Marty, E. Neher, B. Sakmann, F. J. Sigworth, Improved patch-clamp techniques for high-resolution current recording from cells and cell-free membrane patches, Pflügers Archiv 391 (2) (1981) 85–100.
  * [37] D. Werry, J. Eldstrom, Z. Wang, D. Fedida, Single-channel basis for the slow activation of the repolarizing cardiac potassium current, iks, Proceedings of the National Academy of Sciences 110 (11) (2013) E996–E1005.
  * [38] K. S. Cole, J. W. Moore, Potassium ion current in the squid giant axon: dynamic characteristic, Biophysical Journal 1 (1) (1960) 1–14.
  * [39] J. W. Moore, K. S. Cole, Resting and action potentials of the squid giant axon in vivo, The journal of general physiology 43 (5) (1960) 961–970.
  * [40] R. E. Taylor, J. W. Moore, K. S. Cole, Analysis of certain errors in squid axon voltage clamp measurements, Biophysical journal 1 (2) (1960) 161–202.
  * [41] M. F. Priest, E. E. Lee, F. Bezanilla, Tracking the movement of discrete gating charges in a voltage-gated potassium channel, Elife 10.  doi:10.7554/eLife.58148.
URL https://www.ncbi.nlm.nih.gov/pubmed/34779404

  * [42] P. S. Langan, V. G. Vandavasi, W. Kopec, B. Sullivan, P. V. Afonne, K. L. Weiss, B. L. de Groot, L. Coates, The structure of a potassium-selective ion channel reveals a hydrophobic gate regulating ion permeation, IUCrJ 7 (5) (2020) 835–843.
  * [43] N. Schoppa, F. Sigworth, Activation of shaker potassium channels: I. characterization of voltage-dependent transitions, The Journal of general physiology 111 (2) (1998) 271–294.
  * [44] F. Bezanilla, E. Perozo, E. Stefani, Gating of shaker k+ channels: Ii. the components of gating currents and a model of channel activation, Biophysical journal 66 (4) (1994) 1011–1021.
  * [45] J. Tytgat, P. Hess, Evidence for cooperative interactions in potassium channel gating, Nature 359 (6394) (1992) 420–423.
  * [46] L. Delemotte, M. Tarek, M. L. Klein, C. Amaral, W. Treptow, Intermediate states of the kv1. 2 voltage sensor from atomistic molecular dynamics simulations, Proceedings of the National Academy of Sciences 108 (15) (2011) 6109–6114.
  * [47] L. Delemotte, M. A. Kasimova, D. Sigg, M. L. Klein, V. Carnevale, M. Tarek, Exploring the complex dynamics of an ion channel voltage sensor domain via computation, BioRxiv (2017) 108217.
  * [48] A. Peyser, W. Nonner, Voltage sensing in ion channels: Mesoscale simulations of biological devices, Physical Review E 86 (1) (2012) 011910.
  * [49] A. Dryga, S. Chakrabarty, S. Vicatos, A. Warshel, Coarse grained model for exploring voltage dependent ion channels, Biochimica et Biophysica Acta (BBA)-Biomembranes 1818 (2) (2012) 303–317.
  * [50] Z. Song, X. Cao, T.-L. Horng, H. Huang, Electric discharge of electrocytes: Modelling, analysis and simulation, Journal of Theoretical Biology 498 (2020) 110294\.
  * [51] Z. Song, X. Cao, T.-L. Horng, H. Huang, Selectivity of the kcsa potassium channel: Analysis and computation, Physical Review E 100 (2) (2019) 022406.
  * [52] X. Cao, Z. Song, T.-L. Horng, H. Huang, Electric potential generation of electrocytes: Modelling, analysis, and computation, Journal of Theoretical Biology 487 (2020) 110107.
  * [53] Z. Song, X. Cao, H. Huang, Electroneutral models for a multidimensional dynamic poisson-nernst-planck system, Physical Review E 98 (3) (2018) 032404.
  * [54] B. Eisenberg, W. Liu, Poisson–nernst–planck systems for ion channels with permanent charges, SIAM Journal on Mathematical Analysis 38 (6) (2007) 1932–1966.
  * [55] L. Zhang, W. Liu, Effects of large permanent charges on ionic flows via poisson–nernst–planck models, SIAM Journal on Applied Dynamical Systems 19 (3) (2020) 1993–2029.
  * [56] D. Colquhoun, A. Hawkes, On the stochastic properties of single ion channels, Proceedings of the Royal Society of London. Series B. Biological Sciences 211 (1183) (1981) 205–235.
  * [57] D. Colquhoun, A. G. Hawkes, The principles of the stochastic interpretation of ion-channel mechanisms, in: Single-channel recording, Springer, 1995, pp. 397–482.
  * [58] B. Eisenberg, X. Oriols, D. Ferry, Dynamics of current, charge and mass, Computational and Mathematical Biophysics 5 (1) (2017) 78–115.
  * [59] B. Eisenberg, N. Gold, Z. Song, H. Huang, What current flows through a resistor?, arXiv preprint arXiv:1805.04814.
  * [60] H. Miedema, M. Vrouenraets, J. Wierenga, W. Meijberg, G. Robillard, B. Eisenberg, A biological porin engineered into a molecular, nanofluidic diode, Nano letters 7 (9) (2007) 2886–2891.

##  Appendix A Parameter values

We adopt the following values for the physical parameters [29, 30]

$\displaystyle k_{B}=1.38\times 10^{-23}\mathrm{J}/\mathrm{K},\quad
e_{0}=1.602\times 10^{-19}\mathrm{C},\quad\epsilon_{0}=8.854\times
10^{-12}\mathrm{C}/(\mathrm{V}\cdot\mathrm{m}),$ (45) $\displaystyle
T=292.15\mathrm{K},\quad k_{B}T/e_{0}\approx
25.17\textrm{mV},\quad\epsilon_{r0}=2,\quad\epsilon_{r1}=40,$ $\displaystyle
L=0.75\textrm{nm},\quad s=0.15\textrm{nm},\quad A=(0.7\textrm{nm})^{2},\quad
c_{0}=560~{}\textrm{mM}\approx 3.37\times 10^{26}/\textrm{m}^{3},$
$\displaystyle c_{1}^{L}=10\textrm{mM},\quad c_{2}^{L}=550\textrm{mM},\quad
c_{3}^{L}=560\textrm{mM},$ $\displaystyle c_{1}^{R}=400\textrm{mM},\quad
c_{2}^{R}=160\textrm{mM},\quad c_{3}^{R}=560\textrm{mM},$ $\displaystyle
D_{0}=D_{1}=10^{-10}\textrm{m}^{2}/\textrm{s},\quad
D_{2}=D_{3}=10^{-12}\textrm{m}^{2}/\textrm{s},\quad,D_{b}=10^{-19}\textrm{m}^{2}/\textrm{s},$
$\displaystyle AJ_{0}=AD_{0}c_{0}/L\approx 2.2\times 10^{7}/\textrm{s},\quad
e_{0}AJ_{0}\approx 3.53~{}\textrm{pA},$ $\displaystyle
t_{0}=\frac{L^{2}}{D_{0}}=5.625\times 10^{-9}\textrm{s},\quad
V_{0}=-80\textrm{mV},\quad V_{1}=160\textrm{mV}.$

The dimensionless quantities are

$\displaystyle\epsilon=\frac{\epsilon_{0}k_{B}T}{e_{0}^{2}c_{0}L^{2}}\approx
7.3\times 10^{-3},\quad\beta=LAc_{0}\approx 0.12,$ (46) $\displaystyle
c_{1}^{L}\approx 0.018,\quad c_{2}^{L}\approx 0.982,\quad c_{3}^{L}=1,$
$\displaystyle c_{1}^{R}\approx 0.71,\quad c_{2}^{R}\approx 0.29,\quad
c_{3}^{R}=1,$ $\displaystyle D_{1}=1,\quad D_{2}=D_{3}=0.01,\quad
D_{b}=10^{-9},$ $\displaystyle V_{0}=-3.18,\quad V_{1}=6.36.$

##  Appendix B Continuity of the total current

###  B.1 The continuous system

The total current consists of three different types of current in different
regions

  * (i)

the current from the change of electric field (for the whole interval/channel)

  * (ii)

the current from the ionic fluxes (outside of the bubble)

  * (iii)

the current from the motion of the bubble charge (in the bubble)

We will illustrate the continuity of the total current by the dimensional
system in Section 2.1. For the region outside of the bubble
($-L<x<s_{b},s<x<L$), we define the total current (per unit cross-sectional
area) as $I_{total}^{pnp}$

$\displaystyle
I_{total}^{pnp}(x,t)=\epsilon_{0}\epsilon_{r}\partial_{t}E+\sum_{i=1}^{3}e_{0}z_{i}J_{i}=-\epsilon_{0}\epsilon_{r}\partial_{tx}\phi+e_{0}(J_{1}+J_{2}-J_{3}).$
(47)

Taking the time derivative of $(\ref{eq1})_{1}$ and using $(\ref{eq1})_{2}$,
we get

$\displaystyle\partial_{x}I_{total}^{pnp}=0,$ (48)

which implies the continuity of current outside of the bubble.

In the bubble, the define the total current (per unit cross-sectional area) as

$\displaystyle
I_{total}^{bubble}(x,t)=\epsilon_{0}\epsilon_{r}\partial_{t}E+\partial_{t}{Q}_{b}=-\epsilon_{0}\epsilon_{r}\partial_{tx}\phi+\partial_{t}{Q}_{b},$
(49)

where ${Q}_{b}$ is the total bubble charge (per unit area) stored in the
interval $[s_{b},x]$ (it is the magnitude of total negative charge)

$\displaystyle{Q}_{b}=\int_{s_{b}}^{x}\frac{q_{b}}{V_{b}}dx=\frac{q_{b}(x-s_{b})}{A(s-s_{b})}.$
(50)

If ${Q}_{b}$ increases, that means some positive current of the bubble charge
goes across the interface at $x$. Another interpretation is based on the
velocity of the cross sectional surface at $x$

$\displaystyle v(x)=\frac{(s-x)}{(s-s_{b})}\frac{ds_{b}}{dt},$ (51)

and one can easily verify that

$\displaystyle\partial_{t}{Q}_{b}=J_{b}=-\frac{q_{b}}{V_{b}}v(x).$ (52)

Taking the time derivative of (3), we get the continuity of the total current
in the bubble

$\displaystyle\partial_{x}I_{total}^{bubble}=0.$ (53)

###  B.2 The discrete numerical scheme

The quantities $q_{k}^{n+1}$ and $\epsilon_{r,k+1/2}$ are defined as

$\displaystyle
q_{k}^{n+1}=\int_{x_{k-1/2}}^{x_{k+1/2}}\frac{1}{\beta}\frac{{q}_{b}}{(s-s_{b})}dx,$
(54)

and

$\displaystyle\epsilon_{r,k+1/2}=\frac{h}{\frac{h_{1}}{40}+\frac{h_{2}}{2}},\quad
h_{1}=s_{b}-x_{k},\quad
h_{2}=x_{k+1}-s_{b},\quad\mathrm{if}\,s_{b}\in[x_{k},x_{k+1}].$ (55)

Next we show the continuity of the total current in the discrete scheme. The
equation of $\phi$ can be written as

$\displaystyle
J_{\phi,k+1/2}^{n+1}-J_{\phi,k-1/2}^{n+1}=hc_{1,k}^{n+1}+hc_{2,k}^{n+1}-hc_{3,k}^{n+1}-hq_{k}^{n+1},$
(56) $\displaystyle
J_{\phi,k+1/2}^{n+1}=-\epsilon\epsilon_{r,k+1/2}\frac{\phi_{k+1}^{n+1}-\phi_{k}^{n+1}}{h}.$

Summing over $k=1,..,N$ gives

$\displaystyle
J_{\phi,N+1/2}^{n+1}-J_{\phi,1/2}^{n+1}=h\sum_{k=1}^{N}c_{1,k}^{n+1}+\sum_{k=1}^{N}hc_{2,k}^{n+1}-\sum_{k=1}^{N}hc_{3,k}^{n+1}-\sum_{k=1}^{N}hq_{k}^{n+1}.$
(57)

Repeat it for $J^{n}$ with the time step $t_{n}$ , take the difference, divide
it by $\Delta t$, and then we get

$\displaystyle\frac{J_{\phi,N+1/2}^{n+1}-J_{\phi,N+1/2}^{n}}{\Delta
t}-\frac{J_{\phi,1/2}^{n+1}-J_{\phi,1/2}^{n}}{\Delta t}$ (58) $\displaystyle=$
$\displaystyle h\sum_{k=1}^{N}\frac{c_{1,k}^{n+1}-c_{1,k}^{n}}{\Delta
t}+h\sum_{k=1}^{N}\frac{c_{2,k}^{n+1}-c_{2,k}^{n}}{\Delta
t}-h\sum_{k=1}^{N}\frac{c_{3,k}^{n+1}-c_{3,k}^{n}}{\Delta t}$
$\displaystyle-\frac{h}{\Delta
t}\left(\sum_{k=1}^{N}q_{k}^{n+1}-\sum_{k=1}^{N}q_{k}^{n}\right)$
$\displaystyle=$
$\displaystyle-\left(J_{1,N+1/2}^{n+1}-J_{1,1/2}^{n+1}\right)-\left(J_{2,N+1/2}^{n+1}-J_{2,1/2}^{n+1}\right)+\left(J_{3,N+1/2}^{n+1}-J_{3,1/2}^{n+1}\right)$
$\displaystyle-\frac{h}{\Delta t}\left(Q_{b}^{n+1}-Q_{b}^{n}\right).$

Rearranging the terms leads to

$\displaystyle\frac{J_{\phi,N+1/2}^{n+1}-J_{\phi,N+1/2}^{n}}{\Delta
t}+\left(J_{1,N+1/2}^{n+1}+J_{2,N+1/2}^{n+1}-J_{3,N+1/2}^{n+1}\right)+\frac{h}{\Delta
t}\left(Q_{b}^{n+1}-Q_{b}^{n}\right)$ (59) $\displaystyle=$
$\displaystyle\frac{J_{\phi,1/2}^{n+1}-J_{\phi,1/2}^{n}}{\Delta
t}+J_{1,1/2}^{n+1}+J_{2,1/2}^{n+1}-J_{3,1/2}^{n+1},$

where the three terms on the left hand side are the discrete version of the
three types of currents defined in (i), (ii), (iii) in the previous
subsection. If the sum is over the entire interval (i.e., $N+1/2$ is the right
end), the term $Q_{b}^{n+1}-Q_{b}^{n}$ disappears since the total bubble
charge $Q_{b}^{n}$ is conserved by definition, and the total current is
conserved at the two ends.

##  Appendix C Derivation for quasi-static state and steady state

###  C.1 The quasi-static state

For the quasi-static state, we ignore the dipole and $V_{0}$ (equivalently the
$\phi$ is shifted up by a constant $V_{0}$ and continuity condition of $\phi$
will be used at interface). We first consider the right part $s<x<1$. We get

$\displaystyle c_{1}=c_{1}^{R}e^{-(\phi-V_{1})},\quad
c_{2}=c_{2}^{R}e^{-(\phi-V_{1})},\quad c_{3}=e^{\phi-V_{1}},$ (60)
$\displaystyle\epsilon\epsilon_{r1}\phi^{\prime\prime}=e^{\phi-
V_{1}}-e^{-(\phi-V_{1})},$

where $c_{3}^{R}=c_{1}^{R}+c_{2}^{R}=1$ have been used. Integrating once gives

$\displaystyle\frac{1}{2}\epsilon\epsilon_{r1}[(\phi^{\prime}(x))^{2}-(\phi^{\prime}(s+))^{2}]=e^{\phi-
V_{1}}+e^{-(\phi-V_{1})}-B_{2},$ (61)

where

$\displaystyle
B_{2}=e^{\phi_{s}-V_{1}}+e^{-(\phi_{s}-V_{1})},\quad\phi^{\prime}(s+)=\frac{\epsilon_{r0}}{\epsilon_{r1}}\phi^{\prime}(s-)=\frac{\epsilon_{r0}}{\epsilon_{r1}}\tilde{\phi}_{s}.$
(62)

Then, we obtain

$\displaystyle(\phi^{\prime}(x))^{2}=G_{1}(\phi;\phi_{s},\tilde{\phi}_{s})$
$\displaystyle=\left(\frac{\epsilon_{r0}}{\epsilon_{r1}}\tilde{\phi}_{s}\right)^{2}+\frac{2}{\epsilon\epsilon_{r1}}[e^{\phi-
V_{1}}+e^{-(\phi-V_{1})}-B_{2}]$ (63)

and

$\displaystyle
x=\int_{\phi_{s}}^{\phi}\frac{1}{\sqrt{G_{1}(\phi;\phi_{s},\tilde{\phi}_{s})}}d\phi+s.$
(64)

For the left part $-1<x<s_{b}$, we get

$\displaystyle\epsilon\epsilon_{r1}\phi^{\prime\prime}=e^{\phi}-e^{-\phi},$
(65)
$\displaystyle\frac{1}{2}\epsilon\epsilon_{r1}[(\phi^{\prime}(x))^{2}-(\phi^{\prime}(s_{b}-))^{2}]=e^{\phi}+e^{-\phi}-B_{3},$

with

$\displaystyle
B_{3}=e^{\phi_{s_{b}}}+e^{-\phi_{s_{b}}},\quad\phi_{s_{b}}=B_{1}(s_{b}-s)^{2}+\tilde{\phi}_{s}(s_{b}-s)+\phi_{s}$
(66)
$\displaystyle\phi^{\prime}(s_{b}-)=\frac{\epsilon_{r0}}{\epsilon_{r1}}\phi^{\prime}(s_{b}+)=\frac{\epsilon_{r0}}{\epsilon_{r1}}[\tilde{\phi}_{s}+2B_{1}(s_{b}-s)],$

where $B_{1}$ is given in (35). Then, we get

$\displaystyle(\phi^{\prime}(x))^{2}=G_{2}(\phi;\phi_{s},\tilde{\phi}_{s})$
$\displaystyle=(\phi^{\prime}(s_{b}-))^{2}+\frac{2}{\epsilon\epsilon_{r1}}[e^{\phi}+e^{-\phi}-B_{3}]$
(67)

and

$\displaystyle
x=-\int_{\phi_{s_{b}}}^{\phi}\frac{1}{\sqrt{G_{2}(\phi;\phi_{s},\tilde{\phi}_{s})}}d\phi+s_{b}.$
(68)

###  C.2 The steady state

Now we consider the steady state. Define

$\displaystyle p=c_{1}+c_{2},\quad J_{p}=J_{1}+\frac{J_{2}}{D_{2}},$ (69)

then the two equations for $c_{1}$ and $c_{2}$ lead to

$\displaystyle-J_{p}=p^{\prime}(x)+p\phi^{\prime}(x).$ (70)

Let $V=V_{0}+V_{1}$. For the right part $x>s$, we get

$\displaystyle c_{3}=c_{3}^{R}e^{\phi-V}=e^{\phi-V}.$ (71)

Multiplying $\phi^{\prime}$ in the equation of $\phi$ (i.e., equation
$(\ref{eq40})_{1}$, and the delta function is put into the jump conditions)
gives

$\displaystyle-\epsilon\epsilon_{r1}\phi^{\prime\prime}(x)\phi^{\prime}(x)=p\phi^{\prime}-c_{3}\phi^{\prime}=-J_{p}-p^{\prime}-c_{3}^{\prime},$
(72)

and integrating gives

$\displaystyle
p=\frac{1}{2}\epsilon\epsilon_{r1}\left[(\phi^{\prime}(x))^{2}-(\phi^{\prime}(1))^{2}\right]-J_{p}(x-1)-c_{3}+2,$
(73)

where the boundary conditions at $x=1$ have been used. Substituting into
$(\ref{eq40})_{1}$, we get

$\displaystyle-\epsilon\epsilon_{r1}\phi^{\prime\prime}(x)=\frac{1}{2}\epsilon\epsilon_{r1}\left[(\phi^{\prime}(x))^{2}-(\phi^{\prime}(1))^{2}\right]-J_{p}(x-1)-2(e^{\phi-V}-1),$
(74)

for $x>s$. Similarly for the left part $x<s$, we have

$\displaystyle c_{3}=c_{3}^{L}e^{\phi},\quad
c_{1}=\frac{1}{2}\epsilon\epsilon_{r1}\left[(\phi^{\prime}(x))^{2}-(\phi^{\prime}(-1))^{2}\right]-J_{p}(x+1)-c_{3}+2,$
(75)
$\displaystyle-\epsilon\epsilon_{r1}\phi^{\prime\prime}(x)=\frac{1}{2}\epsilon\epsilon_{r1}\left[(\phi^{\prime}(x))^{2}-(\phi^{\prime}(-1))^{2}\right]-J_{p}(x+1)-2(e^{\phi}-1).$