MATHEMATICAL MODELING AND SIMULATION OF MECHANOELECTRICAL
TRANSDUCERS AND NANOFLUIDIC CHANNELS
By
Jin Kyoung Park

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Mathematics - Doctor of Philosophy
2014

ABSTRACT
MATHEMATICAL MODELING AND SIMULATION OF
MECHANOELECTRICAL TRANSDUCERS AND NANOFLUIDIC
CHANNELS
By
Jin Kyoung Park
Remarkable advances in nanotechnology and computational approaches enable researchers
to investigate physical and biological phenomena in an atomic or molecular scale. Smallerscale approaches are important to study the transport of ions and/or molecules through ion
channels in living organisms as well as exquisitely fabricated nanofluidic channels. Both
subjects have similar physical properties and hence they have common mathematical interests
and challenges in modeling and simulating the transport phenomena. In this work, we first
propose and validate a molecular level prototype for mechanoelectrical transducer (MET)
channel in mammalian hair cells. Next, we design three ionic diffusive nanofluidic channels
with different types of atomic surface charge distribution, and explore the current properties
of each channel.
We construct the molecular level prototype which consists of a charged blocker, a realistic
ion channel and its surrounding membrane. The Gramicidin A channel is employed to
demonstrate the realistic channel structure, and the blocker is a positively charged atom of
radius 1.5˚
A which is placed at the mouth region of the channel. Relocating this blocker
along one direction just outside the channel mouth imitates the opening and closing behavior
of the MET channel. In our atomic scale design for an ionic diffusive nanofluidic channel, the
atomic surface charge distribution is easy to modify by varying quantities and signs of atomic
charges which are equally placed slightly above the channel surface. Our proposed nanofluidic

systems constitutes a geometrically well-defined cylindrical channel and two reservoirs of KCl
solution.
For both the mammalian MET channel and the ion diffusive nanofluidic channel, we employ
a well-established ion channel continuum theory, Poisson-Nernst-Planck theory, for three
dimensional numerical simulations. In particular, for the nano-scaled channel descriptions, the
generalized PNP equations are derived by using a variational formulation and by incorporating
non-electrostatic interactions. We utilize several useful mathematical algorithms, such as
Dirichlet to Neumann mapping and the matched interface and boundary method, in order to
validate the proposed models with charge singularities and complex geometry. Moreover, the
second-order accuracy of the proposed numerical methods are confirmed with our nanofluidic
system affected by a single atomic charge and eight atomic charges, and further study the
channels with a unipolar charge distribution of negative ions and a bipolar charge distribution.
Finally, we analyze electrostatic potential and ion conductance through each channel
model under the influence of diverse physical conditions, including external applied voltage,
bulk ion concentration and atomic charge. Our MET channel prototype shows an outstanding
agreement with experimental observation of rat cochlear outer hair cells in terms of open
probability. This result also suggests that the tip link, a connector between adjacent stereocilia,
gates the MET channel. Similarly, numerical findings, such as ion selectivity, ion depletion
and accumulation, and potential wells, of our proposed ion diffusive realistic nanochannels
are in remarkable accordance with those from experimental measurements and numerical
simulations in the literature. In addition, simulation results support the controllability of the
current within a nanofluidic channel.

To my parents, Jongyoul Park and Sunkyoung Choi, and my sister, Minkyoung Park

iv

ACKNOWLEDGMENTS

Over the past eight years, it has been a challenging and rewarding journey for me to complete
the PhD study. I would never have been able to finish my dissertation without the guidance
of my committee members, support from my family, and help from my friends and colleagues.
Within the limited space, I wish to acknowledge those who contribute to my achievement.
I would like to thank all my comprehensive exam and dissertation defense committee
members, Dr. Guowei Wei, Dr. Chichia Chiu, Dr. Jongeun Choi, Dr. Andrew Christlieb,
Dr. Russell Schwab and Dr. Moxun Tang for their valuable time and support. Especially, I
would like to express my deepest and sincerest gratitude to my advisor, Dr. Guowei Wei,
for his thoughtful guidance, encouragement, caring and patience. He has provided me a
systematic training in computational skills and a broad understanding of several research
fields, including Mathematics, Biology, Chemistry, Physics and Engineering. He also led me
to build an interdisciplinary perspective. Five years’ study with him taught me to have great
passion in research.
I would like to express immeasurable gratitude to my lovely family. My father is always
my first and everlasting mentor, my mother supports me in every way and my wonderful
sister is my lifetime companion.
Especially, I deeply appreciate my precious friends, Ms. Yeoreum Lee and Mr. Daewoo
Pak for helping with the preparation for my dissertation defense presentation. I also would
like to thank Dr. Sanghyup Jeong, who has been a mentor in that he gave me considerate
advice on overcoming difficulties as a PhD student.
Finally, I would like to thank my other friends as a colleague. Firstly, I appreciate Dr.
Kelin Xia for his valuable help in programming tools and research tips whenever I met
v

difficulties in my work. Many thanks also go to the rest of former and current members in
Dr. Wei’s group, Dr. Duan Chen, Dr. Zhan Chen, Dr. Langhua Hu, Mr. Yin Cao and Ms.
Weijuan Zhou, for their cooperation. Further, I would like to express gratitude to my other
MSU friends, Dr. Kyungbae Park, Ms. Hana Cho, Ms. Bosu Choi, Mr. Seonghak Kim, Mr.
Sugil Lee and Mr. Woong Bae Park, for their encouragement.

vi

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

x

Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Introduction to mechanoelectrical transducer channel in auditory system
1.2 Introduction to nanofluidic channels . . . . . . . . . . . . . . . . . . . . .
1.3 Review of current computational models . . . . . . . . . . . . . . . . . .
1.3.1 Three representative computational approaches . . . . . . . . . .
1.3.1.1 Molecular dynamics . . . . . . . . . . . . . . . . . . . .
1.3.1.2 Stochastic dynamics . . . . . . . . . . . . . . . . . . . .
1.3.1.3 Continuum models . . . . . . . . . . . . . . . . . . . . .
1.3.2 Modeling and simulation of biological ion channels . . . . . . . . .
1.3.3 Modeling and simulation of synthetic nanofluidic channels . . . .
1.4 Existing challenges . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4.1 Existing challenges in mechanoelectrical tranducer channel . . . .
1.4.2 Existing challenges in nanofluidics . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

1
1
4
10
11
11
12
14
18
21
23
23
24

Chapter 2 Theoretical models and mathematical algorithms
2.1 Generalized Poisson-Nernst-Planck (PNP) model . . . . . .
2.1.1 Computational domain . . . . . . . . . . . . . . . . .
2.1.2 Energy functional . . . . . . . . . . . . . . . . . . . .
2.1.3 Governing equations . . . . . . . . . . . . . . . . . .
2.2 Mathematical algorithms . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

26
26
27
29
34
38

Chapter 3 Modeling and simulating a mechanoelectrical transducer (MET)
channel in mammalian hair cells in molecular level . . . . . . . .
3.1 A molecular level prototype of a MET channel . . . . . . . . . . . . . . . . .
3.2 Computational setup for the PNP systems . . . . . . . . . . . . . . . . . . .
3.3 MET channel prototype simulations . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Experimental data . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Computational results . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41
42
44
46
46
47
53

Chapter 4 Modeling and simulating three types of ionic diffusive nanofluidic channels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1 Atomic scale design of an ionic diffusive nanofluidic channel . . . . . . . . .

56
57

vii

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

4.2

4.3

4.4

Numerical validation . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Analytical solution system . . . . . . . . . . . . . . . .
4.2.2 A cylindrical nanochannel with a single atomic charge .
4.2.3 A cylindrical nanochannel with eight atomic charges . .
4.2.4 A negatively charged ionic diffusive nanofluidic channel
4.2.5 A bipolar ionic diffusive nanofluidic channel . . . . . .
Nanofluidic simulations . . . . . . . . . . . . . . . . . . . . . .
4.3.1 A negatively charged ionic diffusive nanofluidic channel
4.3.2 A bipolar ionic diffusive nanofluidic channel . . . . . .
4.3.3 A double-well ionic diffusive nanofluidic channel . . . .
Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

Chapter 5 Thesis achievements and future work . . . . . . . . . . . . . . . .
5.1 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.1 On the modeling and simulating a MET channel in mammalian hair cells
5.1.2 On the modeling and simulating charged nanofluidic channels . . . .
5.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
APPENDICES . . . . . . . . . . . . . . . . . . . . . .
A
Mathematical algorithms . . . . . . . . . . . . . . . . . .
A.1
Dirichlet to Neumann mapping . . . . . . . . . .
A.2
Matched interface and boundary (MIB) method .
A.3
Iterations of Poisson and Nernst-Planck equations
B
Unit system . . . . . . . . . . . . . . . . . . . . . . . . .
C
Atomic charge distribution . . . . . . . . . . . . . . . . .
C.1
A negatively charged channel . . . . . . . . . . .
C.2
A bipolar channel . . . . . . . . . . . . . . . . . .
C.3
A double-well channel . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

59
60
61
62
64
66
68
69
76
82
85
88
89
89
91
92
94
95
95
96
99
102
105
105
106
107

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

viii

LIST OF TABLES

Table 4.1

Numerical errors and orders for the cylindrical channel with a single
atomic charge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63

Numerical errors and orders for the cylindrical channel with eight
atomic charges. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63

Table 4.3

Numerical errors and orders for the negatively charged channel. . . .

67

Table 4.4

Numerical errors and orders for the bipolar channel. . . . . . . . . .

68

Table 1

Centimeter-Gram-Second system of units . . . . . . . . . . . . . . . 102

Table 2

Gaussian units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

Table 3

Positions and charges of all atomic charges in the negatively charged
channel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

Table 4

Positions and charges of all atomic charges in the bipolar channel. . 106

Table 5

Positions and charges of all atomic charges in the double well channel. 107

Table 4.2

ix

LIST OF FIGURES

Figure 2.1

Figure 3.1

Figure 3.2

Figure 3.3

Figure 3.4

An illustration of the computational domain Ω consisting of Ωm and
Ωs . (a) For the MET channel prototype, the domain is divided into
four regions: channel protein, membrane layers, bulk regions and
channel pore. Here, the red circle represents the blocker, a charged
atom, for channel gating effect and the arrow indicates the direction
of its movement. (b) For the ionic diffusive nanofluidic channel, the
ion inclusion region Ωs is composed of two reservoirs and the inside
of the cylindrical channel with a static boundary. Here, red circles
indicate the charged atoms around the channel to generate channel
surface charge effect. . . . . . . . . . . . . . . . . . . . . . . . . . .

27

An illustration of the MET channel model consisting of GA channel
protein, the artificial membrane, and a positive ion which presents
the gating effect. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

Electrostatic potential profiles at different locations xb of the blocker
with charge Q = 2ec when C 0 = 0.1M and Φ0 = 0.1V. The region
between two dashed lines indicates the channel pore region. As the
blocker position gets closer to 0.9˚
A along the x-direction, the blocker
produces a stronger barrier near the channel mouth as shown in the
electrostatic potential profile. . . . . . . . . . . . . . . . . . . . . . .

48

Electrostatic potential profiles along the channel direction for the
MET prototype (a) at different bulk ion concentrations C 0 = 0.1M
(circles) and C 0 = 0.2M (diamonds) under the fixed applied voltage
Φ0 = 0.1V and (b) at different applied external voltages Φ0 = 0.1V
(circles) and Φ0 = 0.2V (diamonds) under the fixed bulk concentration
C 0 = 0.1M. In both investigations, we locate the blocker at xb = 0.9˚
A
and set the blocker charge as Q = 2ec . Increase either in the bulk
concentration or in the external voltage creates higher potential within
the channel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

The total current versus the relative displacement dx of the charged
blocker. (a) Total current behavior under two bulk ion concentrations
C 0 = 0.1M (circles) and C 0 = 0.2M (diamonds) when Q = 2ec and
Φ0 = 0.1V are fixed. (b) Total current behavior under two external
voltages Φ0 = 0.1V (circles) and Φ0 = 0.2V (diamonds) when Q = 2ec
and C 0 = 0.1M are fixed. Consequently, increase either in bulk ion
concentration or in external voltage induces increase in channel current. 50

x

Figure 3.5

(a) Electrostatic potential profiles and (b) total current profiles for
four different charges Q = 0.5ec (triangles), Q = 1ec (diamonds),
Q = 1.5ec (squares) and Q = 2ec (circles) when C 0 = 0.1M and
Φ0 = 0.1V are fixed. The region between two dashed lines indicates
the channel pore region in (a). As the blocker charge gets increased,
the barrier at the gate of the channel gets higher. In conclusion, the
amplitude of the current curve also gets enlarged. . . . . . . . . . .

51

Open probability PO is plotted against the relative displacement dx
under four different atomic charges Q = 0.5ec (triangles), Q = 1ec
(diamonds), Q = 1.5ec (squares) and Q = 2ec (circles) where C 0 =
0.1M and Φ0 = 0.1V are set. The solid dots are experimental data
of the normalized MET current versus normalized rat hair bundle
displacement obtained from the literature [123]. At each blocker
charge, the open probability in response to the relative displacement
forms a sigmoidal shape and, especially, the charge Q = 2ec gives a
remarkable agreement with the experimental result. . . . . . . . . .

53

Sensitivity test of the molecular level MET model for the prediction of
open-closed probability according to the relative blocker displacement
when (a) the bulk ion concentration is doubled; (b) the applied
external voltage is doubled. Both cases are fairly consistent with the
experimental finding [123]. . . . . . . . . . . . . . . . . . . . . . . .

54

Illustration of a nanofluidic system geometry. (a) A 3D view of
a schematic cylindrical nanochannel whose ends are connected to
two reservoirs of KCl solution. (b) A 2D cross-section view of the
cylindrical channel whose diameter is 10˚
A and length is 49˚
A in the xz0
0
plane. Here, ΦL and ΦR , respectively, represent the external voltage
applied at the left and right electrodes, and C 0 represents the bulk
ion concentration of both K+ and Cl− in two reservoirs. . . . . . . .

57

Illustration of atomic charge distribution. (a) A 2D cross-section view
and (b) a 3D view of schematic diagram consisting of the cylindrical
channel and four atomic charges which are equally placed around the
nanochannel at z = 0˚
A. . . . . . . . . . . . . . . . . . . . . . . . . .

58

Figure 4.3

Illustration of the geometries of two simple numerical test cases. . .

62

Figure 4.4

Surface plot of electrostatic potential profiles on (a) a negatively
charged channel and (b) a bipolar channel. As represented in the color
bar, blue colors represent negative values, while red colors represent
positive values. The negatively charged channel surface is mostly blue,
but the color of the bipolar channel surface is changed from red to blue. 64

Figure 3.6

Figure 3.7

Figure 4.1

Figure 4.2

xi

Figure 4.5

Three contour plots of electrostatic potential through the negatively
charged channel at z = −10, z = 0 and z = 10 on the xy-plane when
the mesh size is h = 0.2. The blue colors just outside the channel
imply the negative surface charge. . . . . . . . . . . . . . . . . . . .

66

Figure 4.6

Three contour plots of electrostatic potential through the bipolar
channel at z = −10, z = 0 and z = 10 on the xy-plane when the mesh
size is h = 0.2. As indicated in the color bar, blue colors represent
negative values, while red colors represent some small positive values.
At the cross section of z = −10 the inside of the channel is blue (i.e.,
negatively charged solution), whereas at the cross section of z = 10
the inside of the channel is red (i.e., slightly positively charged solution). 67

Figure 4.7

Effect of applied voltage on electrostatics and dynamics through the
negatively charged nanofluidic channel. (a) Electrostatic potential
profiles and (b) ion concentration distributions along the channel
length (z−axis). Here ∆Φ is varied from 0V (squares), 0.2V (triangles), 0.4V (asterisks), 0.6V (diamonds), 0.8V (circles) to 1V (plus
signs), where ∆Φ is the difference of the applied voltage between
two electrodes. Each surface atomic charge is Qk = −0.08ec and
bulk concentration is C 0 = 0.05M. Two dashed vertical lines indicate
the ends of the nanochannel. As the applied voltage difference gets
larger, the potential at the right part of the inner channel is increased.
Consequently, more K+ ions (dashed line) are accumulated at the
left part of the inner channel. However, there is little change in the
concentration of Cl− ion (solid line). . . . . . . . . . . . . . . . . . .

70

Ionic current for each ion species versus the applied potential difference
∆Φ through (a) the uncharged nanochannel (Qk = 0ec ) and (b) the
negatively charged nanochannel (Qk = −0.08ec ). Here, the bulk ion
concentrations at both reservoirs are fixed at C 0 = 0.05M. When
Qk = 0ec , both current-voltage graphs are linear and the positive
current (squares) is roughly double of the negative current (triangles).
When Qk = −0.08ec , the positive current-voltage graph (squares) is
nonlinear and the negative current-voltage graph (triangles) is almost
always zero. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

Figure 4.8

xii

Figure 4.9

Figure 4.10

Figure 4.11

Figure 4.12

Figure 4.13

Effect of atomic surface charge on ionic current through the negatively
charged nanochannel. The ionic currents in response to the change of
the atomic charge Qk are depicted. Six different charges Qk = 0ec ,
Qk = −0.02ec , Qk = −0.04ec , Qk = −0.06ec , Qk = −0.08ec and
Qk = −0.1ec are simulated under the fixed bulk concentration C 0 =
0.05M and applied voltage difference ∆Φ = 1V. Herein, |Qk | is the
magnitude of the charge of the each atom. Stronger atomic charges
amplify the K+ current (squares) sharply, but reduce the Cl− current
(triangles) to near zero. . . . . . . . . . . . . . . . . . . . . . . . . .

73

Effect of bulk ion concentration on total current through the negatively
charged nanochannel. The graphs of total channel current against
the external voltage difference (I-V) are shown at five different bulk
concentrations C 0 = 0.01M (squares), C 0 = 0.05M (triangles), C 0 =
0.1M (diamonds), C 0 = 0.2M (circles) and C 0 = 0.4M (plus signs)
when Qk = −0.08ec . As the bulk ion concentration is increased, the
total current gets higher and, moreover, the I-V characteristic becomes
near linear. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

74

The effect of bulk ion concentration on the normalized current through
the negatively charged nanochannel. (a) The normalized ionic current
and (b) the normalized channel current with respect to the increase in
bulk ion concentration when the atomic charge Qk = −0.08ec and the
applied voltage difference ∆Φ = 1V are fixed. The negatively charged
channel with a sufficiently larger bulk ion concentration behaves like
a uncharged channel because the normalized values becomes one. . .

75

(a) Electrostatic and (b) ion concentration profiles of the bipolar
nanofluidic channel along the channel axis at (i) no bias ∆Φ = 0V,
(ii) forward bias ∆Φ = 1V and (iii) reverse bias ∆Φ = −1V. Fix the
bulk ion concentration C 0 = 0.1M and the atomic surface charges
|Qk | = 0.08ec . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77

Ionic current versus applied voltage difference in the bipolar nanofluidic channel. The atomic charge |Qk | = 0.08ec and the bulk ion
concentration C 0 = 0.1M are assumed. Under reverse bias, both ion
species cannot pass through the channel. However, under forward
bias, the currents of both ion species get enlarged. Especially, the
cation current (squares) increases faster. . . . . . . . . . . . . . . . .

79

xiii

Figure 4.14

Figure 4.15

Figure 4.16

Figure 4.17

Effect of atomic surface charge on the current through the bipolar
nanofluidic channel. Three sets of atomic charges, that is, |Qk | =
0.04ec (squares), |Qk | = 0.08ec (triangles) and |Qk | = 0.12ec (diamonds), are studied. Here, the bulk ion concentration C 0 is fixed
at 0.1M. All I-V curves increase when ∆Φ varies from −1V to 1V.
Greater magnitude of the atomic charge results in a higher channel
current with forward bias, but insufficient atomic charge may generate
a current leakage because it weakens the depletion zone with reverse
bias. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

Effect of bulk ion concentration on the channel current through the
bipolar nanofluidic channel. Three sets of bulk ion concentrations,
such as C 0 = 0.05M (squares), C 0 = 0.1M (triangles) and C 0 = 0.2M
(diamonds), are explored in terms of the current-voltage relation.
The atomic charges are given by |Qk | = 0.08ec . As the bulk ion
concentration gets higher, the amplitude and gradient of the currentvoltage relationship gets maximized. . . . . . . . . . . . . . . . . . .

81

The schematic diagram of a three-dimensional double-well nanofluidic
channel. The channel length is divided into three parts. The first and
last parts are negatively charged, but the middle part is positively
charged. The red dots indicate atoms with negative charge Qk =
−0.12ec and the green dots indicate atoms with positive charge Qk =
0.04ec . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

Effect of applied voltage on electrostatics and dynamics through the
double-well nanofluidic channel. (a) Electrostatic potential profiles
and (b) ion concentration distributions along the channel length
(z−axis). Here ∆Φ is varied from 0V (squares), −0.2V (triangles),
−0.4V (asterisks), −0.6V (diamonds), −0.8V (circles) to −1V (plus
signs), where ∆Φ is the difference of the applied voltage between two
ends of the system. The bulk concentration is C 0 = 0.05M for both
ions K+ and Cl− . Two dashed vertical lines indicate the ends of
the cylindrical channel. Each electrostatic potential graph shows two
potential wells, which brings about two piles of K+ ions along the
channel axis. Moreover, higher applied voltage at the left end of the
system, Φ0L , breaks the symmetry of the potential wells. The left well
becomes weaker and the right well becomes stronger. . . . . . . . . .

83

xiv

Figure 4.18

Figure 1

Figure 2

Effect of bulk ion concentration on the channel current through the
double-well nanofluidic channel. The graphs of total channel current
against the external voltage difference (I-V) are studied at four different
bulk concentrations: C 0 = 0.01M (squares), C 0 = 0.05M (triangles),
C 0 = 0.1M (diamonds), and C 0 = 0.2M (circles). Higher bulk ion
concentration encourages the total current to increase and the I-V
characteristic to become linear. . . . . . . . . . . . . . . . . . . . . .

84

An illustration of the MIB method for an irregular geometry. The
kth grid line intersects the interface at (xi , yo , zk ) (star), which
involves two irregular points (circle) xi , yj , zk and xi , yj+1 , zk .
Here (xi , yo , zk+2 ) and (xi , yo , zk+1 ) (square) are two auxiliary points
needed to discretize ( Φz )z and (xi−1 , yo , zk ) and (xi−2 , yo , zk ) (square)
are two auxiliary points needed to discretize ( Φx )x . The other points
(triangle) indicate the points selected to interpolate Φ along the y-grid
lines at (xi , yo , zk ), (xi , yo , zk+1 ) and (xi , yo , zk+2 ). . . . . . . . . . .

97

Flowchart of the numerical implementation to solve the standard PNP
system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

xv

Chapter 1
Introduction

1.1

Introduction to mechanoelectrical transducer channel in auditory system

Ion channels are pore-forming membrane proteins which are very permeable and highly
selective and their opening and closing are suitably managed [112]. Ion channels are involved in
diverse physiological functions including altering membrane potentials, controlling electrolyte
movements for cell volume regulation and polarized transport of salt, and generating electrical
signals which are utilized to regulate hormone secretion, neurotransmitter release and muscle
contraction [94]. Traditionally, ion channels are classified and named by their two crucial
properties such as ion selectivity and gating mechanism [94, 98]. Most channels allow only a
few specific ions to penetrate through them although there are exceptions like non-selective
cation channels and, in particular, potassium, sodium, calcium and chloride are four principal
permeable ions. The gating mechanism of an ion channel is classified according to a particular
stimulus which the channel gate react strongly. Typical kinds of stimuli are the voltage
gradient across the transmembrane (voltage-gated), the binding and interaction of ligands
with an ion channel (ligand-gated) and a mechanical force (mechanically-gated). Besides,
other ways to gate ion channels can be triggered by photonic and thermal stimuli.
The mechanoelectrical transducer (MET) channel in hair cells is one of the remarkable ion

1

channel research subjects. Hair cells in the inner ear play a pivotal role in the mechanotransduction for the sense of hearing. Mechanotransduction is the conversion of a mechanical
stimulus into an electrical signal, which is fundamental for the senses of hearing, touch and
balance [88]. The peripheral auditory system, the sensory system for hearing, is mainly
divided into three parts: the outer, middle and inner ear [21]. A wide range of airborne
sound waves enters the external canal of the outer ear and through the middle ear they are
converted to pressure waves which vibrate the basilar membrane in the cochlea of the inner
ear [159]. In the cochlea, basilar membrane vibrations deflect the stereociliary bundles on
the hair cells and hence the sensory cells are depolarized by allowing cations, predominantly
Ca2+ , to penetrate into the mechanoelectrical transducer (MET) channels [88].
The human cochlea contains three rows of outer hair cells (OHCs) and one row of inner
hair cell (IHC) with different functions and shapes [175]. OHCs amplify the amplitude
as well as the frequency selectivity of basilar membrane vibrations. IHCs convey acoustic
information to afferent neurons [78, 121]. Each hair cell possess a mechanically sensitive
bundle consisting of a number of actin-filled microvilli, also named stereocilia, on its apical
surface [159]. The stereocilia in a bundle are arranged in rows of increasing height and,
moreover, they are tightly connected to each other by extracellular linkages including tip
links and top-connectors [175]. In particular, it is speculated that the tip link, which extends
from the tip of a stereocilium to the side of an adjacent taller stereocilium, directly opens
the mechanically gated ion channels [121, 175]. Additionally, it is postulated that an elusive
elastic element, a so-called gating spring, contributes to unfasten the channel inlet [121].
The biophysical principles underlying the mechanotransduction process in hair cells have
been intensively investigated in the past few decades. The recently proposed model of
activation and adaptation of the hair cell tranducer channel is composed of three stages
2

[121]. It is postulated that deflection of a hair bundle toward the tallest stereocilia, namely
positive deflection, increases the open probability of the MET channel at the lower end of
the tip link. Then Ca2+ ions travel through the channel pore and then contribute to fast
adaptation by binding to a molecule inside or near the channel. Finally, slow adaptation
including channel closure and tension restoration is accomplished by a myosin motor at the
upper end of the tip link. Despite these advances in MET channel study, the molecular detail
of the mechanotransducer machinery has been proved to be elusive. Particularly, there is no
direct structural confirmation for the molecular building blocks of the MET channel.
There have been several studies to elucidate the features of the MET channel such as gating
mechanism, ion conductance, molecular-level structure and location because charactering the
channel is beneficial to comprehend not only hair cell transducer process but also hearing
mechanism [77]. First of all, the MET channel is a non-selective cation channel with a
considerably high permeability of Ca2+ [156]. Then many empirical experiments have been
focused on the localization of this mechanosensitive channel, which is essential to validate
that the channel is gated by mechanical stimuli [62, 109]. Especially, Beurg et al. established
that the MET channels were located only at the bottom of the tip links as indicated by
fast confocal calcium-imaging [16]. In spite of these findings, the MET channel is a notable
example of ion channels whose molecular nature is still obscure, partly due to few channels
per hair cell and few hair cells per organ [159]. Although various potential candidates
including any known mechanosensitive channels and transient receptor potential channels
were compared with the MET channel from biophysical perspectives in the literature, none of
them was perfectly matched with the MET channel [77, 159]. However, Ricci et al. showed
that the conductance of hair cells in the turtle cochlea was tonotopically changed with their
characteristic frequency, which implies that the MET channel might be composed of several
3

subunits [165]. Farris and his colleagues measured the dimensions of the MET channel using
various antagonists of the candidate channel classes and simple amine compounds [74]. Then
they found that the minimal pore size was about 12˚
A and channel length was approximately
31˚
A [74]. More intrinsic characteristics of the MET channel are yet to be discovered in order
to construct its molecular components.

1.2

Introduction to nanofluidic channels

Nanofluidics, one of the most vibrant scientific research fields, refers to the study and
application of the transport of molecules and/or ions dissolved in a solution as well as fluid
behavior through or past structures with one or more nanometer dimensions [55, 172]. The
nanofluidic studies have been highly motivated by the efforts to design a solid-state DNA
separation system and, moreover, extensive research has been driven by the dramatic advance
in both nanofabrication techniques and theoretical tools to describe fluid motion on the
nanoscale [188]. Another impetus is the discovery of new mechanical and electrochemical
phenomena which are non-existent or less influential in macrofluidic or microfluidic systems
[188]. In a nanofluidic channel, the combination of the remarkably large surface-to-volume
ratio, the electrostatic interactions between the fluid and the charged wall, and the channel’s
characteristic dimension comparable with the size of biomolecules generate unique transport
patterns such as double-layer overlap, ion-current rectification, ion permittivity and diffusion
[68, 223].
The extraordinary transport features of nanofluidics have been received a great deal of
attention in chemistry, physics, biology, material science, medicine and several engineering
fields [172]. In particular, the design and fabrication of nanofluidic devices for molecular

4

biology applications is a new interdisciplinary field that takes advantage of precise control and
manipulation of fluids at submicrometer and nanometer scales in order to study the behavior of
molecular and biological systems. Nanofabrication techniques are generally classified into two
groups: top-down and bottom-up methods [151, 160]. The top-down method creates patterns
on a large scale with nanometer lateral dimensions, while the bottom-up method places
the atoms and molecules in nanostructures by using chemical reactions. Sophisticatedly
synthesized nanochannels are applied to biosensing of high-throughput, regulating and
separating ions and molecules and energy harvesting [93]. As the characteristic length
scale of the fluid is comparable with the length scale of the biomolecule and/or the Debye
length, nanofluidics can be applied to a variety of interesting powerful tools for genomics
or proteomics [1, 154]. For example, Han et al. introduced an entropic trapping sieving
mechanism for long DNA molecules [95] and Fu et al. designed and investigated a patterned
anisotropic nanofluidic sieving structure for size-based separation of DNA and proteins as
well as electrostatic separation of proteins [82]. For a protein analysis, nanofilter is used to
maximize protein concentration in a sample [204]. Schoch and his colleagues demonstrated
a pH-controlled diffusional separation of proteins in a nano-scaled channel [171]. Besides,
nanofluidic techniques have been instrumented for macromolecule accumulator [44, 210],
electronic circuits [116, 118, 215], local charge inversion [96], photonic crystal circuits [73]
and nanofluidic dynamic array [202]. Despite rapid development in nanotechnology, it is
still an intriguing challenge for engineers to generate inexpensive nanostructures which are
more feasible in diverse areas, for example biomedical research [154]. Since nanofluidic
device prototyping and fabrication are technically challenging and financially expensive, it is
desirable to further advance the field by mathematical/theoretical modeling and simulation.
One major factor to characterize a nanofluidic system is its structure. The novel nanofab5

rication skills enable the production of diverse nanofluidic devices from 0, 1 to 2 dimensions,
where the dimension indicates the number of non-nanometer length component of the device
[68]. Generally, 0-D, 1-D, and 2-D nanochannels are also called as nanopores, nanotubes, and
nanoslits, respectively. Another parameter to represent the structure is the ratio of channel
height to width, namely the aspect ratio (AR), and, especially, low AR and near-unity AR
channels are referred to as planar and square channels [1]. Nanopores are usually formed
perpendicularly through diverse substrate materials, and the most well-known nanopore
sensors are pore-forming proteins such as α-hemolysin and silicon nitride membranes with
solid-state pores [196, 223]. Nanopore-based sensing is a cost-effective label-free approach
without amplification at a relatively high speed and it has been promoted to meet the needs
of researchers interested in DNA sequencing [22, 196]. Kasianowicz and his group tried
to detect single-stranded RNA and DNA molecules using α-hemolysin [120] and Li et al.
designed a set of solid-state nanopores in a silicon nitride to observe the folding behavior
of individual double-stranded DNA [134]. Nanochannels whose depth and/or width are
reduced to the nanoscale are eligible to combine with other sophisticated devices and to
demonstrate the transport inside the channel [223]. For example, Perry and his colleagues
designed funnel-shape nanofluidic channels with different taper angles and proved that the
taper angle of the funnel influenced the ion-rectifying effect [161]. Typically, nanochannels
with comparatively longer lateral length have either a cylindrical or a conical geometry [223].
While surface charge and applied external potential predominantly control the counterion
current in a cylindrical channel, the flow direction is also critical to determine the ion conductance patterns through a conical channel. In the work presented here, we will concentrate
more on nanofluidic channels.
The other important property that distinguishes a nanofluidic system from a higher
6

dimensional fluid system is its unprecedented transport phenomena. The fluid transport
through a nanoscale channel is explained by three factors: the external driving force generated
by either an electrical potential gradient or a pressure gradient, a variety of molecular-level
interactions, and the friction forces caused by solvent-wall interactions and hydrodynamic
solute-wall interactions [188]. Due to the extraordinary surface-to-volume ratio in nanofluidic
systems, fluid interactions with the wall dominantly induce notable transport behaviors
including ion permselectivity, ion enrichment and depletion, fast mixing and rapid injection
of a small amount of reagent [162]. Especially, steric/hydration interactions, van der Waals
interactions and electrostatic interactions play a central role in nanofluidic systems and
these intermolecular interactions determine the characteristic length scale of nanofluidics [55].
The interactions within a nanoscale fluid system are mostly influenced by its physical and
chemical properties such as geometrical confinement and charge, and flow condition such
as ion composition and concentration. Consequently, microscopic interactions dominate the
solute and solvent transport in nanofluidic structures, but continuum fluid mechanics regulate
it in macrofluidic and microfluidic structures.
Electrostatic distribution in a solution and electrokinetic flow of the solute molecules are
fundamental concepts to understand the transport in nanofluidic devices [172]. Specifically, the
surface charge of the channel wall plays a pivotal role to derive electrostatic interactions and
electrokinetic effects within nanoscale fluidic channels when charged particles are sufficiently
adjacent to the wall [56, 173, 190, 218]. The fixed charge at the channel wall induces
electrostatic interactions with ions dissolved in a solution; hence the wall attracts oppositely
charged ions, but repels ions with the same charge [188]. The interplay between the long-range
electrostatic forces and short-range van der Waals forces creates a screening region, so-called
an electrical double layer (EDL), in an electrolyte solution confined to the nanostructure so
7

that this oppositely charged region is able to attain electroneutrality at the interface between
the solid and the liquid [172]. The combination effect between two forces at the charged
solid-liquid interface is well-established by the Derjaguin-Landau-Verwey-Overbeek (DLVO)
theory [63, 197]. The structure of the EDL is described well in the Gouy-Chapman-Stern
model, which suggests three layers of the EDL: the inner Helmholtz plane, outer Helmholtz
plane and diffuse layer [172]. While the inner Helmholtz plane consists of non-hydrated coions
and counterions that are attached to the channel surface, the outer Helmholtz plane contains
bound, hydrated and partially hydrated counterions. Moreover, the part between the inner
and outer Helmholtz planes is called the Stern layer, where the charge and potential are
linearly varied. The diffuse layer, the farthest layer, contains mobile coions and counterions.
In particular, a thorough understanding of the potential variance in accordance with the
structure of the EDL is indispensable because the principal mechanism of the transport is
electrokinetic behavior [162]. The readers can refer to the literature [6, 91, 157] for a more
comprehensive understanding of the EDL.
Characteristic length scale, such as Reynolds number, Biot number and Nusselt number,
is a reference quantity used to explain particular characteristics in a physical system. For
example, Reynolds number, the ratio of inertial forces to viscous ones, is a dimensionless
number to determine flow patterns in fluidic mechanics [105]. One of the most representative
characteristic length scales in nanofluidic structures is the Debye length λD =
where ε is the dielectric constant of the solvent, ε0 is the permittivity of vacuum,

εε0 kB T
,
Cα0 qα2
α
kB is the

Boltzmann constant, T is the absolute temperature, and Cα0 and qα are, respectively, the bulk
ion concentration and the charge of ion species α [55]. The Debye length only depends on the
flow conditions, including the solute composition and concentration, and also describes the

8

thickness (or, precisely, 1e th reduction) of the EDL where the potential decays exponentially
[172]. Essentially, ionic fluid behaves like a microscopic flow within the EDL region, but it
acts as a macroscopic flow far beyond the Debye length.
When the internal diameter of a nanochannel is comparable to or smaller than the Debye
length, the electrolyte solution becomes a unipolar solution with a charge of sign opposite to
the channel surface charge [55]. However, this unipolar transport never occurs in microfluidic
systems because the channel pore is much larger than the Debye length. If the EDLs are
overlapped across the channel pore dimension, only the control of counterion flow is possible,
but otherwise the flow of both counterions and coions can be governed [56]. Ion selectivity is
another important feature which enables nano-sized channels to work as an ionic filter and it
is defined as the ratio of the difference between currents of cations and anions to the total
current delivered by both ions [199]. Vlassiouk and his group examined the ion selectivity
of single nanometer channels under various conditions including channel dimension, buffer
concentration and applied voltage.
The rectification of ionic current in nanofluidic devices has been of great interest because
it is possible to control the ionic flow by simple electrical switches [41]. This asymmetry in ion
conductance can be derived by applying the external voltage when the symmetry of surface
charge distribution, bulk concentration, channel geometry or a combination of these is broken
along the longitudinal axis [41]. In order to interpret the rectification phenomena through
nanofluidic systems, several experiments and theoretical modeling have been performed
in various types of structures. For example, Pu et al. conducted experiments to observe
ion-enrichment and ion-depletion effects in 60-nm-deep nanochannels, which usually resulted
in the ion-current rectification [163]. In their design, an applied field gave rise to accumulation
of both cations and anions at the cathode and absence of all ions at the anode of the channels.
9

Siwy and her colleagues showed the rectification phenomena in a single conical nanopore and
in a single conical nanotube when the surface was adequately charged [184, 183, 180, 181, 182].
Karnik et al. designed a cylindrical nanofluidic diode with uneven surface charge distribution
[117] and Daiguji et al. theoretically described the ionic transport through such channels
[56].

1.3

Review of current computational models

Channels with nanoscale pores are ubiquitous in biological systems, while recently astonishing
advance in synthesis techniques enables manufacturing such small-sized pore channels for
various applications[153]. Modeling and analyzing the ionic transport and fluidic dynamics
through both natural and artificial nanofluidic channels have been a great challenge and
a crucial task over the past few decades. Membrane channels are obviously different from
solid-state channels. For instance, the ionic selectivity of membrane channels is mainly based
on the sizes of ions and channel pore; in contrast, that of solid-state channels depends on
electrostatic effects [55]. In spite of such differences, ion transport through both types of
nanochannels can be simulated using similar theoretical and computational tools due to
common physical features.
An enormous number of theoretical approaches, from phenomenological to fundamental,
have been proposed and developed in order to elucidate the transport phenomena within
cell membrane channels as well as man-made channels [46, 55, 112, 133, 146, 153, 168]. In
this section, we account for three well-known representative theoretical approaches: all-atom
molecular dynamics (MD), stochastic dynamics and continuum models such as PoissonBoltzamann (PB) and Poisson-Nernst-Planck (PNP) equations. Then we discuss the history

10

of modeling biological ion channels and introduce several theoretical models and simulations
of ion channels using each model in the literature. Finally, we summarize a list of modeling
and simulation of synthetic nanofluidic channels.

1.3.1

Three representative computational approaches

1.3.1.1

Molecular dynamics

In the MD simulation, the dynamical motions of all the atoms, including ions, water, protein
and lipid in the simulated system, are described via Newton’s second law of motion, and
the empirical force fields express the potential energy under the influence of interatomic
interactions [153].
mi

d2
ri (t) = Fi (t) = −∇V (ri (t), . . . , rN (t)),
dt2

where ri and mi denote, respectively, the position and the mass of the atom i, N the total
number of atoms in the system and Fi the force acting on the atom i that is the gradient of
the potential energy V (ri (t), . . . , rN (t)).
The potential energy is composed of bonded potential and unbonded potential [146].
While the former is related to bond lengths, bond angles, improper dihedral angles and
torsional angles, the latter contains the Coulomb potential for electrostatic interactions and
the Lennard-Jones potential for van der Waals interactions. The widely used molecular
dynamics simulation packages with a set of force fields are CHARMM [23, 143], AMBER
[51, 209], GROMOS [97], and OPLS-AA [113].
The greatest advantage of the MD simulation is to give the most accurate time-dependent
atomic detail of the system of interest which is also comparable with experimental data
[146]. In this respect, it is a powerful theoretical approach for biophysical and biochemical
11

analysis. Especially, the MD calculates a valuable parameter, the potential of mean force
(PMF). It is very beneficial to discover the permeation through nanometer pores or channels
and, moreover, it can provide the selectivity sequences of monovalent cations using the
free-energy profiles [48]. Despite such great benefits, the MD has several limitations: the
high computational cost due to microscopic detailedness of description, the wide range of
time scales from femtoseconds over miliseconds to seconds, the difficulty to account for
polarization and pH value, and the difficulty to apply an electrostatic potential across the
system [153]. The abrupt advance of computer powers raises the time scales of simulation up
to 1 microsecond [125]. Ab initio MD framework has progressed to analyze the polarization
effects in biomolecular systems and, for instance, these simulations successfully demonstrate
dissociation of NaCl in water [148, 192]. Many studies have purported to surmount the other
puzzling issue of the availability of the external electrostatic potential in MD simulations
[153]. Since Aksimentiev and Schulten simulated the ion conductance through α-Hemolysin
under the applied external electrostatic field [4], a similar investigation method has been
utilized to study ion transport in a variety of systems [20, 126, 145, 152, 179, 174].

1.3.1.2

Stochastic dynamics

The stochastic dynamics scheme is proposed as a great compromise between two computational
approaches in that it can simulate only a reasonable number of ions explicitly by considering
the solvent molecules implicitly [55]. The simplest and most commonly used form of stochastic
dynamics to explore complex many-body systems is the Brownian dynamics (BD) [48].
Specifically, for the BD simulations of nanofluidic channels, the dynamics depict the behavior
of each solute ion while considering the channel as a rigid structure and the solvent molecules
implicitly under frictional and stochastic forces [146]. By reducing the number of degrees of
12

freedom in calculation, this statistical approach provides very useful computational data such
as current-voltage and conductance-concentration curves at a relatively lower cost [48, 153].
Herein, the Langevin equation is repeatedly solved to detect the random motions of the
individual ions of interest under the influence of the dissipative and fluctuating forces [48]:

mi

d2
ri (t) = Fi (t, r) − mi γi r˙ i (t) − ξi (t),
dt2

where Fi denotes the systematic force including electrostatic forces derived from interactions
between the charged atoms, mi γi r˙ i (t) the frictional force denoting mi the mass, γi the
reciprocal of the relaxation time and r˙ i the velocity, and ξi the stochastic force. Two
necessary parameters required to conduct a BD simulation are the diffusion coefficient for
each ion and the force applied to each ion which corresponds to the multi-ion PMF [146].
These two factors are usually obtained from all-atom MD simulations.
In particular, the BD simulations are very valuable for ion channel studies in that
they empower to calculate the current flow and ionic concentrations in the system using a
trajectory analysis and to determine the valence selectivity [48]. Although the BD somewhat
compensates the defects of both the MD and the continuum approaches, this framework
also has critical shortcomings and hence many scholars have extensively researched solutions
to such problems [48, 146, 153]. The stochastic dynamics of the reduced system involves
two major simplifications assuming water as a continuum and protein as a fixed structure,
which are closely involved in the shortcomings of the simulations [133]. Another weakness of
the BD method is using the position-independent dielectric constant for water, but the MD
demonstrates that the polarizability of water is decreased due to confinement in a real system
[48]. Thus, some BD simulations use position-dependent dielectric constants validated by

13

MD simulations [107, 131]. For ion channel analysis, the conformational changes of proteins
may have influence on ion permeation, which needs to be discovered using experiments and
MD simulations and, further, to be applied in BD simulations [48].
The dynamic Monte Carlo (DMC) is an alternative statistical theory to simulate nonequilibrium or relaxation systems [169]. The fundamental principle of the DMC is that the
dynamics of a system can be described by a sequence of states generated by random particle
displacements in time [54, 169]. In this simulation, a randomly-chosen ion of an arbitrary
species is transferred to a randomly-chosen new location within a maximum displacement of
its former location [54]. Cs´cnyi et al. satisfactorily computed ionic current through a sodium
channel model under nearly physiological bulk ion concentration [54]. A theoretical approach
coupling grand canonical Monte Carlo (GCMC) and Brownian dynamics (BD) is proposed by
Roux and his colleagues in order to incorporate bulk ion concentration and transmembrane
potential in ion channel microscopic simulations [107, 108]. They validated the proposed
algorithm to study ion permeation and selectivity of the OmpF porin of Escherichia coli.

1.3.1.3

Continuum models

In a continuum based simulation, the mathematical theory treats all the components of
a system of interest including ions, water and proteins as a continuous element [48]. The
most well-established continuum models are Poisson-Boltzmann (PB) for electrostatics and
Poisson-Nernst-Planck (PNP) for ion transport [146]. Due to such simplifying approximation,
the continuum approach has a number of limitations; however, it plays a crucial role to
predict structural and physical features from the current-voltage-concentration relation at
high computational efficiency [146].
In the early 1900s, the PB equation was first proposed by Gouy (1910) and Chapman (1913)
14

in order to elucidate electrostatic patterns around a biomolecule [29, 90]. It is generalized
as a continuum approach to compute electrostatic free energy of small spherical ions in an
ionic solution by Debye and H¨
uckel in solution chemistry [61]. Since such worthwhile usage
of the PB equation in electrochemistry, this model has been studied and applied in various
research fields as well [139]. For example, it is the foundation of DLVO theory which describes
electrostatic effects in colloid systems [63, 197] and, moreover, it is also utilized in biophysics
[59, 103]. The non-linear PB equation is usually written in the following form:

qα Cα0 exp

−∇ · ( (r)∇Φ(r)) = 4πρm (r) + 4πλ(r)
α

−qα Φ(r)
kB T

,

where (r) denotes the spatial-dependent dielectric constant, Φ(r) the electrostatic potential
and ρm (r) the fixed charge density of the biomolecule. For each mobile ion species α within
an aqueous solution, Cα0 and qα represent, respectively, the concentration at bulk region and
the charge. Moreover, kB is the Boltzmann constant, T is the absolute temperature, and λ(r)
is the characteristic function whose value is 1 in the solvent region which ions can penetrate
through, but whose value is 0 in the biomolecule region impermeable to ions.
The fundamental assumption of the PB model is that the electrostatic potential in an
ionic solution is determined by the Boltzmann distribution of solute charges [81, 146]. The
PB equation is largely solved to compute the electrostatic potential at the solvent-accessible
molecular surface, the reaction rates between molecules in a solution, the free energy of
biomolecular association and its salt dependence and pKa shifts in proteins [81]. Furthermore,
it can be conveniently integrated into molecular mechanics and dynamics to maximize
computational efficiency [81]. Since the PB is involved with a mean-field approximation,
the PB modeling produces several weaknesses [26]. First, the important properties of an

15

ion, such as the size and the discrete surface charge, are not considered while the ion is
simply considered as a point charge. Second, the non-electrostatic ion-ion interactions and
ion-solid interactions are ignored. Third, the dielectric constant of the solvent is assumed
to be the same through the system, but in reality the permittivity of the medium can be
altered due to charge distribution near the surface. Such problems notwithstanding, the
PB framework is definitely a useful tool in ion channel studies in that it facilitates the
calculation of the free energy required to move an ion from the bulk region to the channel
inside, the characterization of ion-channel interactions and the description of transmembrane
electrostatic potential distribution [146].
The utmost remarkable continuum theory for ion transport through a nanofluidic channel
is an electrodiffusion model, namely Poisson-Nernst-Planck (PNP) theory, which treats
both solutes and solvent as continuous entities, but describes membrane protein in atomic
detail. Specifically, the PNP model describes the solvent water molecule as a dielectric
continuum, treats ion species by continuum density distributions and, in principle, retains
the discrete atomic detail and/or charge distribution of the channel or pore constitution
[14, 72, 129, 221, 220]. The Nernst-Planck equation (NP) combines two fundamental physical
laws, that is, Ohm’s law and Fick’s law, based on the fact that the flux of ions through a
nanochannel is driven by the potential and concentration gradients [48]. For each solvent ion
species α,
Jα (r) = −Dα (r) ∇Cα (r) +

Cα (r)
∇ (qα Φ(r) + Uα (r)) ,
kB T

with the Boltzmann constant kB and the absolute temperature T . Herein, Jα (r), Dα (r),
Cα (r), qα and Uα (r) are, respectively, the flux density, the diffusion coefficient, the concentration, the charge valence and the sum of potentials induced from interactions with the

16

channel for ion species α. Additionally, the electrostatic potential Φ(r) is obtained from the
Poisson equation:
−∇ · ( (r)∇Φ(r)) = 4πρm (r) + 4π

qα Cα (r),
α

where (r) and ρm (r) are the position-dependent dielectric constant and the partial charge density of the nanopore, respectively. By solving the NP and Poisson equations self-consistently
and simultaneously, one can calculate the electrostatic potential, the ionic concentration
and the ionic currents [153]. To simulate the PNP model with a nanofluidic channel, the
required parameters need to be gained empirically or, sometimes, by all-atom simulations
[48, 153]. The parameters include the channel geometry and its surface charge composition,
the dimension and bulk ion concentration at each reservoir, the dielectric constant at each
region, the diffusion coefficients of an ion species and the applied external potential [48].
Although the PNP theory is proved to analyze and predict ion transport at a comparably
small computational cost, it is similar with the PB theory in its limitations. In particular, the
PNP theory is inadequate to describe the transport phenomena through narrow nanochannels
owing to the finite size effects of the ions [114] and short-range non-electrostatic interactions
[52].
The Poisson-Nernst-Planck theory, combined with the density functional theory, is a
useful extension of the PNP model [86, 87]. In this framework, ions are considered as charged
rigid spheres and the density functional theory is employed to calculate the chemical potential
of the ions. Other potent theoretical tools are Poisson-Boltzmann-Nernst-Planck (PBNP),
which simplifies the calculation for multiple ion species with Boltzmann distributions of
ion concentrations [140, 221] and variational multiscale models for charge transport, which
employ diverse variational formulations based on differential geometry theory [208].

17

1.3.2

Modeling and simulation of biological ion channels

In 1952, Hodgkin and Huxley designed a voltage clamp experiment technique which is used
to determine the membrane permeability [100]. Using this brilliant idea, they elucidated the
action potential in squid giant axons by analyzing electrical movements of Na+ and K+ through
the ion channels [99]. Similarly, early studies on ion channels were biophysical experiments
to study conductance by measuring voltage and current [19]. Then the understanding
of ion channels at molecular level have been expanded by structural studies [19]. For
example, electron microscopy, nuclear magnetic resonance (NMR) spectroscopy and X-ray
crystallography have been established to reveal atomic-level structural information of ion
channels [112]. The discovery of the crystallized structure of the Streptomyces lividans KcsA
potassium channel by MacKinnon and his collaborators was a remarkable achievement, which
also helps to understand the ion conductance mechanism at atomic basis [67, 144]. As another
good example of the well-studied ion channel structure, the pore and channel molecular
structures of a Gramicidin A (GA) channel have been explored by X-ray crystallographic
and/or NMR spectroscopic methods [201]. The molecular structures of other potassium
channels [135], mechanosensitive channels [9, 28] and a chloride channel [69] have been
discovered. Therefore, the abrupt advance in crystallographic analysis gives rise to atomicresolution structures of many ion channels [48]. The detected structures are used to identify
their functions from the essential laws of physics in electrolyte solutions [48].
Parallel to progress in experimental methods, theoretical tools have been so sophisticated that researchers could investigate structure, function, dynamics and transport of ion
channels in electrophysiology, biochemistry, molecular biology, computational chemistry
and bioinformatics from the whole cell studies to single channel studies [112]. Theoretical

18

modeling is beneficial to obtain an overall physical description of an ion channel, which
gives several advantages [46, 168]. First of all, it is very effective not only to organize and
visualize experimental data but also to predict physical phenomena. Moreover, it possibly
provides a connection between the structure of an ion channel and its functions. A number of
theoretical frameworks have been proposed and developed to answer the following puzzling
issues: the detailed process of ion permeation, the positions of binding sites in a channel, the
rate-limiting steps in conduction, and the relationship between the molecular composition of
an ion channel and its ion selectivity [153]. Powerful computational approaches including
all-atom molecular dynamics (MD), Brownian dynamics (BD), Poisson-Boltzmann (PB)
model, Poisson-Nernst-Planck (PNP) model, Poisson-Boltzmann-Nernst-Planck (PBNP)
model and variational multiscale models have been established over years to explore ion
channels [112, 146, 153]. Each method has strengths and weaknesses, so a researcher should
choose an appropriate method by considering the computational complexity and cost [146].
MD, a dynamic description of the detailed motions of individual atoms in a system, has
been developed to deal with more complex biological or chemical systems [5, 119, 149]. Since
the first simulation of a folded globular protein, its usage has been extended to account
for dynamic conformational changes and longer time scales [5, 119, 149]. Mackay and
his coworkers applied MD method to observe the transport of four ions (lithium, sodium,
potassium, and cesium) through the GA channel, which was the first MD simulation of
an ion channel [142]. Especially, this theoretical approach elucidates the mechanism of
ion selectivity, for instance, water permeation through human aquaporin-1 [60] and Ca2+
permeation through an L-type calcium channel [7]. BD describes the ions of interest explicitly,
but it treats the solvent molecules implicitly [146]. Chung et al. preformed BD simulations
to study currents of sodium and chloride ions through the designed membrane channel
19

whose structure is similar to ACh channel [47]. PB is an equilibrium approach beneficial
to describe the electrostatics around charged biomolecules in a solution, while PNP is a
nonequilibrium approach useful to explore the ion fluxes depending on concentration gradient
and electrostatic potential gradient [146]. Due to high efficiency of computational cost, PNP
has been applied to study diverse ion channels including GA channel [31, 49, 101, 129] and
OmpF porin [3, 107]. Zheng and her colleagues suggested PBNP model in order to simplify
the calculations of Nernst-Planck equations and tested that model with two protein molecules
such as cytochrome c551 and GA [221].
Recently, a variety of differential geometry based multiscale models were introduced
for charge transport by our group [207, 208, 206]. The differential geometry theory of the
molecular surface provides a natural means to separate the microscopic domain of biomolecules
from the macroscopic domain of solvent so that appropriate physical laws are applied to
appropriate domains. Our variational formulation is able to efficiently bridge macro-micro
scales and synergically couple macro-micro domains [206]. One type of our multiscale models
is the combination of Laplace-Beltrami (LB) equation, PB equation and Kohn-Sham equations
to study the molecular mechanism of proton transport [30, 33]. Another type incorporates LB
equation into the generalized PNP equations for the dynamics and transport of ion channels
and transmembrane transportors [208, 206]. The other type alternates MD simulations
and continuum elasticity (CE) descriptions of the solute molecule as well as continuum
fluid mechanics formulation of the solvent [207, 208, 206, 213]. We proposed the theory
of continuum elasticity with atomic rigidity (CEWAR) to treat the shear modulus as a
continuous function of atomic rigidity so that the dynamics complexity of a macromolecular
system is separated from its static complexity [213]. As a consequence, the time-consuming
dynamics is approximated by the CE theory, while the less time-consuming static analysis is
20

carried out in atomic level. Efficient geometric modeling strategies associated with differential
geometry based multiscale models have been carried out with both Lagrangian Eulerian
[75, 76] and Eulerian representations [212].

1.3.3

Modeling and simulation of synthetic nanofluidic channels

Nanofluidics has been extensively studied in chemistry, physics, biology, material science,
and many areas of engineering [172]. The modeling and simulation of nanofluidic systems
have been of enormous importance and a fast growing research field in the past decade.
Basically, modeling a fluidic system with one or more dimensions at nanoscale is different
from modeling a microfluidic system in that the liquid confined to such a narrow system
contains fewer molecules and the wall considerably influences on the transport [189]. Thus
the simulation models utilized in microfluidic systems need to be carefully tested if each
model is applicable in nanofluidic analysis [189]. Most microfluidic devices involve a fluid
flow through a microfluidic pore or channel, which has been an interesting concern in the
theoretical modeling. Gad-al-Hak classified the fluidics models as two categories: continuum
models including Euler, Navier-Stokes (NS) and Burnett equations, and molecular models
including molecular dynamics (MD), direct simulation Monte Carlo (DSMC) and latticeBoltzmann method (LBM) [83]. The characteristic length scales such as Renolds number
and Debye Screening length are key factors to explore the transport features in nanofluidics
as well as microfluidics [55, 105]. Especially, for a nanofluidic device, if the characteristic
length scale of its pore is greater than 5nm, the PNP equations and Stokes equations well
describe electrokinetics and liquid flow within the device, respectively [55, 189]. Otherwise,
the transport analysis requires the discreteness of substances; thus, in particular, MD is an
efficient tool in this respect [55].
21

In every fluidic system, the transport behavior is determined by the forces between the
individual atoms, which implies that considering all the atomic interactions is the best way
to depict fluid mechanics [70]. In spite of several critical limitations such as time scale, MD is
the foremost modeling tool to study the structural and dynamic characteristics of nanoscale
devices [55]. For nanofluidic analysis, MD simulations have been employed to discover
significant dynamics phenomena such as ion current rectification and charge inversion in
various types of solid-state nanochannels [153]. Lyden-Bell et al. designed a simple cylindrical
channel model to investigate the mobility and solvation of ions and uncharged molecules in
an aqueous solution [141]. This study found the solvation shell of water molecules by means
of the MD simulations. All-atom molecular dynamics simulations have been performed to
delineate the ionic transport and several remarkable phenomena within silica devices [153].
Lorenz et al. investigated ζ-potential, Stern layer conductance, charge inversion and ionic
mobilities using a detailed description of ionic transport, electro-osmotic flow and streaming
currents [136, 137]. Curz-Chu et al. demonstrated the effect of ion-binding spots at the
silica surface on the ionic current rectification [53]. Moreover, Shirono and his colleagues
successfully provided a atomic-scale description of ionic transport within ultra-narrow silica
nanopores with 2nm diameter pore [178]. Although MD is a powerful theoretical approach,
it is inadequate to examine a dilute fluid because the intermolecular interactions have little
influence and hence alternative theoretical approaches are mandatory [105]. For example, a
lattice-Boltzmann method (LBM) based on mesoscopic kinetic equations with microscopic
models has been advanced to simulate not only microfluidic systems but also nanofluidic
systems [37, 105].
The Navier-Stokes (NS) equations have been a famous continuum dynamics model to
elucidate microfluidic phenomena in lap-on-a-chip systems at comparatively less computational
22

cost and its simpler version is the Stokes equations [105]. The Navier-Stokes/Stokes equations
have been coupled with various other equations such as PB equations for electrokinetic flow
dynamics and PNP equations for ion transport analysis [105]. Especially, the PNP is the most
renowned model for charge transport in nanofluidics [14, 27, 57, 71]. Daiguji and his coworkers
intensively characterized a variety of transport phenomena through nanofluidic channels
including unipolar solution, energy conversion and ionic rectifying effect by coupling the NS
equations with the PNP equations [57, 58, 56]. Choi et al. investigated the electrokinetic
flow-induced currents through silica nanofluidic channels using the performances of the PB
and PNP models with appropriate boundary conditions [43]. They also validated the proposed
models by comparing with experimental data The BD simulation of ions in a nanopore channel
was combined with the continuum PNP model for regions away from the nanopore channel
[2]. A further simplified model is the Lippmann-Young equation, which is able to predict the
liquid-solid interface contact angle and interface morphology under an external electric field
[186].

1.4
1.4.1

Existing challenges
Existing challenges in mechanoelectrical tranducer channel

Mechanosensitive channel is one of the most famous targets of computational studies, but
the molecular mechanism of mechanosensation is still elusive [146]. For only few simple
mechanosensitive channels, their crystal structures were discovered and they have been studied
to elucidate ion conductance and the mechanism of mechanical gating using simulations
[9, 28, 187, 191]. Since the transduction process of hair cells attracted a great deal of

23

attention, many principal physiological properties of a hair cell mechanotranducer have been
demonstrated [79]. These features contain exceptional sensitivity to a mechanical force,
kinetic principles such as activation and adaptation, ion selectivity, and conductance of a
single MET channel. Especially, the channel permeation properties are critical foundations to
build its molecular building blocks [158]. However, the molecular identity of the MET channel
is not confirmed yet because it is fairly formidable to connect properly its physiological
functions to intrinsic cytoarchitecture and accessory proteins [158]. Parallel to this difficulty,
there are two obscure, but essential, questions - whether the transducer channel is activated
by mechanical stimuli through accessory proteins or by the deformation change in the lipid
bilayer and which part of the channel protein calcium ions bind to for fast adaptation [79].
It is one of the most challenging as well as indispensable works to complete the molecular
structure of the MET channel, which will provide more substantial understanding of the
transduction process in hearing.

1.4.2

Existing challenges in nanofluidics

Nanofluidics has been intensively and extensively developed not only in terms of fabrication
skills but also in terms of theoretical studies over the past decades. In nanofluidic modeling,
computation and analysis, there are many standing theoretical and technical problems.
For example, nanofluidic processes may induce structural modifications and even chemical
reactions [115, 195], which are not described in the present nanofluidic simulations. The
PNP model, one of the most widely used theories in nanofluidics, can incorporate atomic
charge details in its pore or channel description, which are vital to channel gating and fluid
behavior. However, atomic charge details beyond the coarse description of surface charge
are usually neglected in most nanofluidic simulations. Although Stern layer and ion steric
24

effect are significant for the EDL as discussed earlier, they are not appropriately described
in the conventional PNP model. Furthermore, nanofluidic simulations have been hardly
performed in three dimensional realistic settings with physical parameters. Consequently,
simulation results can only be used for qualitative (i.e., phenomenological) comparison and
not for quantitative prediction. Finally, the material interface induced jump conditions in
the Poisson equation are seldom enforced in nanofluidic simulations with realistic geometries.
Therefore, it is imperative to address these issues in the current nanofluidic modeling and
simulation.

25

Chapter 2
Theoretical models and mathematical
algorithms
This chapter contributes to derive a generalized Poisson-Nernst-Planck (PNP) model using a
variational formulation in order to describe the ion transport within not only a biological
ion channel but also a synthetic nanopore channel. In this model, we consider van der
Waals interactions as part of our free energy functional to describe solvent-solute and ion-ion
non-electrostatic interactions. Then we develop Dirichlet to Neumann mapping (DNM) for
dealing with charge singularities and the matched interface and boundary (MIB) method
for material interfaces in order to compute the PNP equations with 3D irregular channel
geometries and singular charges.

2.1

Generalized Poisson-Nernst-Planck (PNP) model

Although the PNP theory is quite standard to describe the dynamics within nanochannels
[14, 35, 36, 45, 50, 110, 194, 199, 200, 203, 206, 222, 225], it does not include non-electrostatic
interactions. In this section, we present a generalized PNP theory by incorporating nonelectrostatic interactions between the solution and the nanoscale channel pore, and between
solvent molecules such as waters and ions. We utilize a variational formulation to derive
generalized PNP equations.
26

2.1.1

Computational domain

(a) The molecular-level MET channel prototype

(b) Ionic diffusive nanofluidic channel

Figure 2.1: An illustration of the computational domain Ω consisting of Ωm and Ωs . (a) For
the MET channel prototype, the domain is divided into four regions: channel protein, membrane
layers, bulk regions and channel pore. Here, the red circle represents the blocker, a charged atom,
for channel gating effect and the arrow indicates the direction of its movement. (b) For the ionic
diffusive nanofluidic channel, the ion inclusion region Ωs is composed of two reservoirs and the
inside of the cylindrical channel with a static boundary. Here, red circles indicate the charged atoms
around the channel to generate channel surface charge effect.

Let us consider a total computational domain Ω ⊂ R3 . We denote Ωm and Ωs respectively
the microscopic channel region (or the ion exclusion region) and the solution region (or the
ion inclusion region). Essentially, we seek continuum descriptions of ions in the solution
domain Ωs , while we consider a discrete atomistic description in the channel structure domain
Ωm . Interface Γ separates two regions Ωm and Ωs so that Ωm

Γ

Ωs = Ω. Since biological

channels and artificial channels differ in the specific composition of the domain Ω, the domain
setting is a careful process. Moreover, it is critical, not simple, to determine a well-defined
interface because it is related to the structural features of the channel.
In this work, we simulate two nano-scaled systems: the molecular prototype for the
MET channel and the ionic diffusive nanofluidic channels with atomic surface charge. In the
27

following, we will discuss in detail about the computational domain for each channel type.

Biological ion channel domain setting
As depicted in Fig. 2.1(a), the molecule region Ωm consists of channel proteins and membrane
layers, but the solvent region Ωs is composed of channel pore and bulk regions which ions
and molecules travel through. Since in reality, the shape of a protein channel continuously
varies, for example by its interactions with water molecules, identifying a static structure
has been a puzzling issue in a theoretical modeling of a complex macromolecule [166]. The
molecular surface (MS) is one of the standard methods to visualize the molecule-solvent
interface [166, 170]. The minimal molecular surface (MMS) is proposed in order to minimize
the surface free energy [10, 11]. To describe the boundary of our MET prototype, we employ
the hypersurface function S : R3 → R3 from our earlier differential geometry based multiscale
models [207, 208, 206]. Basically, the function S(r) is a solvent-solute characteristic function
whose value is 1 in Ωs and 0 in Ωm , whereas it takes proper values between 0 and 1 near
the boundary of the biomolecule to become a smooth function [206]. Thus, the hypersurface
function determines movable biomolecule-solvent interface whose dynamics is driven by
mechanical and electrostatic forces, for example using the Laplace-Beltrami equation [208].
With this hypersurface function S(r), we can derive a generalized PNP model using the
differential geometry formalism of solvent-solute interface as in [206].

Synthetic nanofluidic channel domain setting
In solid-state channel analysis, the region Ωs is divided into the channel inside and two
reservoirs as in Fig. 2.1(b). Usually, an external field leads ions and molecules to pass
through the channel inside from one reservoir to another. Herein, we introduce a domain

28

characteristic function S : R3 → R3 such that Ωm = SΩ and Ωs = (1 − S)Ω. Obviously, S
and (1 − S) are the indicators for the channel domain and the solution domain, respectively.
Unlike the charge and material transport in biomolecular systems, the charge and material
transport in nanofluidic systems induces a negligible reconstruction of the solid-fluid interface
compared to the system scale. Therefore, the fixed solid-fluid interface is predetermined in
the present model as opposed to using the hypersurface function S(r) in ion channel model.
With this domain characteristic function S(r), we present a new variational derivation of the
PNP model for the solid-state nanofluidic channel in this thesis work.

2.1.2

Energy functional

Generally, a nanochannel, which can be either biological or artificial, system is mainly composed of an ionic solution and a microscopic channel composition. Inter-particle correlations
occur between components of the solution or between the solution component and the channel
component at the channel surface. These interactions are mainly classified to be either
electrostatic or non-electrostatic.

Electrostatic energy functional
Electrostatic interactions are ubiquitous in all biomolecular systems and processes, and
are also the dominating effects for nanofluidic behaviors [55, 81]. Especially, electrostatic
interactions in biomolecules are salient to elucidate the several behaviors such as folding,
binding, conformational stability and enzyme activity [155]. The electrostatic interactions are
typically modeled by a number of theoretical approaches: the Poisson-Boltzmann (PB) theory
[59, 81, 130, 177], the polarizable continuum solvation theory [150, 193] and the generalized
Born model, an approximate continuum solvent model, [8, 66]. Among these methods, the PB
29

theory is the most popular in structural and functional analysis of biomolecules. A variation
formulation of this theory was originally introduced by Sharp and Honig in 1990 [176], and
was extended to a multiscale formalism [207, 206] as well as an electrostatic force derivation
[89].
In the present work, we consider the following electrostatic energy functional

Gelectr =






χ −

m

2

|∇Φ|2 + Φρm + (1 − χ) −

s

2

Ns

|∇Φ|2 + Φ
α=1



qα Cα  dr,


(2.1)

where Φ is the electrostatic potential, and s and m are the dielectric constants of the ionic
Nm

Qk δ(r − rk ) with Qk denoting the

solution and channel regions, respectively. ρm =
k=1

partial atomic charge at position rk and Nm the total number of fixed charges on the channel
surface. For each ion species α, Cα and qα , respectively, represent the concentration and the
charge valence which is zero for an uncharged solution component, and Ns represents the
number of mobile ion species in the solution domain. Moreover, the function χ(r) in Eq.(2.1)
is the characteristic function to determine the channel-solution interface. This characteristic
function means the hypersurface function S(r) for membrane channel modeling, but the
domain characteristic funciton S(r) for solid-state channel modeling.

Non-electrostatic interactions
The examples of non-electrostatic interactions are van der Waals interactions, dispersion
interactions, ion-water dipolar interactions, ion-water cluster formation or dissociation, steric
effects, and so forth. In particular, these interactions near the solution-channel interface
have a crucial influence on the transport features due to confinement or dilution of solution;
however, continuum approaches ignore such important correlations [12, 26, 52, 114]. This

30

simplification may result in problematic discrepancy between simulation predictions and
experimental data and hence the extensions of continuum models including non-electrostatic
correlations have been explored. For example, the PB and PNP, two representative continuum
based theories, have been studied in terms of ionic finite size effects in the past [13, 12, 24,
30, 92, 106, 127, 132, 198].
To account for solvent-channel interactions as well as solvent-solvent interactions, the
non-electrostatic interaction energy functional takes the form

Gnon−electr =

(1 − χ) U dr.

U dr =
Ωs

(2.2)

Ω

In this modeling, we assume that the aqueous environment has multiple ion species labeled
by α. The non-electrostatic interaction potential U can be obtained by summing all the
pairwise inter-particle interactions within the solution or near the solution-channel interface:

U (r) =

Cα (r)Uα (r) =
α

Cα (r)Uαk (r) +
k

Cα (r)Uαβ (r),

(2.3)

β

where Cα (r) is the density of αth ionic component of the solution, which may be either
charged or uncharged, and Uα is a non-electrostatic interaction potential with αth species.
For each α, Uαk (r) is an interaction potential with the kth atom of the channel structure
and Uαβ (r) is a potential for ion-ion or ion-water non-electrostatic interactions with βth
component of the solution.
The solvent-solute interactions in solvation analysis have been represented by the LennardJones (LJ) potential [40, 38, 39]. The Weeks-Chandler-Andersen (WCA) decomposition [205]

31

was utilized to split the LJ potential into attractive and repulsive parts [39, 208]



 − αk ,
0 < r − rk < σk + σα ,
att,WCA
Uαk
(r) =


 V LJ ,
r − rk
σk + σα , and
αk


LJ +

 Vαk
αk , 0 < r − rk < σk + σα ,
rep,WCA
Uαk
(r) =


 0,
r − rk
σk + σα ,

(2.4)

where αk is the depth of the potential well, σk and σα are respectively the radii of the
kth atom in the channel region and the αth ion species in the solvent region, and rk
represents the location of the kth atom in the channel. For the MET channel analysis,
rk represents the position of the partial charge of the channel protein; however, for the
ionic diffusive channel analysis, rk represents that of the atomic surface charge. Here,
LJ = 4
Vαk
αk

σ 12
−
rd

σ 6
rd

represents the Lennard-Jones potential for the interaction

between the kth atom in the channel region and the αth ion species in the solvent region
with denoting σ the distance at which the interaction potential is zero and rd the distance
between two particles.
The inter-particle interaction term Uαβ (r) within the solution does not affect the derivation
and the form of other expressions. More detailed description of Uαβ (r) for ion channel
transport can be found in our earlier work [33, 208].

Chemical potential related free energy
Chemical potential related free energy combines a reference homogeneous term and an
entropic term of mobile charges [80], which is essential for the description of mobile charges

32

in the nanofluidic system.
Ns

µ0α − µα0 Cα + kB T Cα ln

Gchem =
α=1

Cα
Cα0

− kB T (Cα − Cα0 ) dr,

(2.5)

where kB is the Boltzmann constant and T is the temperature. For each αth species, µα0 is
a relative reference chemical potential which implies the difference between the equilibrium
concentrations of different solvent species. µ0α and Cα0 , respectively, represent a reference
chemical potential and its associated ion concentration, given that Φ = Uα = µα0 = 0.
Especially, the term kB T Cα ln CCα is the entropy of mixing and the term −kB T (Cα − Cα0 )
α0
is a relative osmotic term [147].
For each species α, it is standard that the variation of the energy Gchem with respect to
Cα , namely,

δGchem
, generates the chemical potential as follows [208].
δCα
µchem
= µ0α − µα0 + kB T ln
α

Cα
Cα0

.

(2.6)

Note that, at equilibrium, µchem
is nontrivial, and Cα and Cα0 are not equal due to possible
α
external electrical potentials, solvent-solute interactions, and charged particles.

33

Total free energy functional
The total free energy functional for the nanofluidic system is

GPNP
total [χ, Φ, {Cα }]



Ns

χ − m |∇Φ|2 + Φρm + (1 − χ) − s |∇Φ|2 + Φ
=
qα Cα 

2
2
α=1

(2.7)

Ns

+(1 − χ)

Cα Uα
α=1
Ns

µ0α − µα0 Cα + kB T Cα ln

+(1 − χ)
α=1

Cα
Cα0

− kB T (Cα − Cα0 ) + λα Cα




dr.



This shows three contributions to the total free energy: the electrostatic free energy including
both fixed and mobile charges, the free energy of non-electrostatic interactions, and the
chemical potential related energy. Note that the electrostatic free energy functional is the
same as the polar solvation free energy functional in the literature [38, 39, 206]. Here λα
is a Lagrange multiplier, which is included to ensure appropriate physical properties at
equilibrium [80].
A powerful feature of the present total free energy function formulation (2.7) is to employ
the free energy functional of the non-electrostatic interactions U for nanofluidic dynamics.
Therefore, the present theory is able to deal with a variety of non-electrostatic interactions
and the solvent microstructure near the channel can be predicted correctly as well.

2.1.3

Governing equations

The total free energy functional (2.7) is a function of electrostatic potential Φ and ion
concentration Cα . The governing equations for this thesis work are derived by using the
34

variational principles [207, 208, 206].

Generalized Poisson equation
The variation of the total free energy functional (2.7) with respect to the electrostatic potential
Φ results in the classical Poisson equation
Nm

−∇ · ( (r)∇Φ(r)) = χ(r)

Ns

Qk δ(r − rk ) + (1 − χ(r))

qα Cα (r),

r ∈ Ω,

(2.8)

α=1

k=1

where (r) = χ(r) m + (1 − χ(r)) s is the dielectric profile whose value is m in the channel
domain Ωm and s in the solution domain Ωs . On the right-hand side, the first term accounts
for the fixed charge density of the channel and the second term represents the mobile charge
density within the solution. Due to the characteristic function χ(r), the Poisson equation
(2.8) can be split into two equations
Nm

−∇ · ( m ∇Φ(r)) =

Qk δ(r − rk ),

r ∈ Ωm

(2.9)

qα Cα (r),

r ∈ Ωs .

(2.10)

k=1
Ns

−∇ · ( s ∇Φ(r)) =
α=1

Since the electrostatic potential Φ(r) is defined on the whole computational domain Ω, the
following jump conditions at the solution-channel interface Γ are to be implemented to ensure
the well-posedness of the generalized Poisson equation [85, 102, 217].

[Φ(r)] = Φm (r) − Φs (r) = 0

(2.11)

[ (r)∇Φ(r)] · n = m ∇Φm (r) · n − s ∇Φs (r) · n = 0,

(2.12)

35

where n is the outward unit norm of the interface Γ. In Eq. (2.8), the densities of ions Cα
should be determined by the variational principle for the generalized Nernst-Planck equation
in the following subsection.
The boundary conditions of Eq. (2.8) are typically mixed boundary conditions depending
on experimental settings. In our numerical experiments, we use the Dirichlet boundary
condition Φ(r) = Φ0 on ∂Ω and the Neumann boundary condition in other types of boundaries,
where Φ0 is the external voltage applied at electrodes. However, the boundary conditions
except the applied voltages are somewhat irrelevant as long as boundaries are sufficiently far
from the channel pore [220, 221].

Generalized Nernst-Planck equation
gen

For each ion species α, we can derive the relative generalized potential µα

by the variation

of the total free energy functional GPNP
total with respect to the ion density Cα , that is,

gen

µα

= µ0α − µα0 + kB T ln

Cα
Cα0

δGPNP
total ,
δCα

+ qα Φ + Uα + λα
(2.13)

= µchem
+ qα Φ + Uα + λα .
α
gen
At equilibrium, we require not µchem
but µα to vanish, which results in
α

λα = −µ0α

and Cα = Cα0 exp −

gen

Consequently, the relative generalized potential µα

gen

µα

= kB T ln

Cα
Cα0

36

qα Φ + Uα − µα0
kB T

.

(2.14)

is simplified as

+ qα Φ + Uα − µα0 .

(2.15)

We derived a similar quantity from a slightly different perspective in our earlier work [221]. It
gen

is interesting to note that the relative generalized potential µα

consists of contributions from

the entropy of mixing, the electrostatic potential, the non-electrostatic interaction with αth
ion species and the position-independent relative reference chemical potential. In practice,
the biological or synthetic nanopore system is out of equilibrium due to applied external field
gen

and/or inhomogeneous concentration across the nanofluidic channel, which implies that µα
does not vanish.
The Fick’s first law generates the ion diffusive flux
gen

Jα (r) = −Dα (r)Cα (r)∇

µα (r)
kB T

(2.16)

with Dα denoting the diffusion coefficient of species α and the general conservation law gives
the diffusion equation
∂Cα (r)
= −∇ · Jα (r)
∂t

(2.17)

[128]. Therefore, the generalized Nernst-Planck equation is

∂Cα (r)
Cα (r)
= ∇ · Dα (r) ∇Cα (r) +
∇(qα Φ(r) + Uα (r))
∂t
kB T

,

α = 1, · · · , Ns , (2.18)

where qα Φ + Uα can be regarded as the potential of the mean field. At the absence of
non-electrostatic interactions, Eq. (2.18) reduces to the standard Nernst-Planck equation. At
the steady state, one has

∇ · Dα (r) ∇Cα (r) +

Cα (r)
∇(qα Φ(r) + Uα (r))
kB T

37

= 0,

α = 1, · · · , Ns .

(2.19)

Herein, Eq. (2.18) does not involve the characteristic function χ because it has already been
restricted to the solution domain Ωs even though the generalized Poisson equation (2.8) is
defined on the whole domain Ω.
Mixed boundary conditions are also applied to the generalized NP equations (2.18) for all
ion species. Since the Nernst-Planck equation is only defined in the ion inclusion region Ωs ,
it is enough to consider the boundary ∂Ωs . On ∂Ω ∩ ∂Ωs , Cα = Cα0 , where Cα0 represents the
bulk ion concentration for species α. However, at the interface Γ = Ωm ∩ Ωs , the zero-flux
condition − Dα ∇Cα +

Cα
∇(qα Φ + Uα )
kB T

= 0 is required, where 0 = (0, 0, 0) is a null

vector.

2.2

Mathematical algorithms

In the present study, we only solve the standard Poisson-Nernst-Planck equations in a selfconsistent manner with appropriate initial/boundary conditions. To emphasize the primary
effects of atomic charges in both molecular-level prototype for MET channel and ionic diffusive
nanofluidic channel design, we neglect the non-electrostatic interactions U . However, since
non-electrostatic interactions are important for nanoscopic scale systems [13, 12], the situation
that U = 0 will be investigated in our future work.
Moreover, it is an indispensable task to establish rigorous numerical schemes and computational algorithms in the study of realistic biological, physical, chemical and engineering
problems [208]. Therefore, we employ several novel theoretical schemes to obtain more
effective and reasonable calculation results from the PNP model. They are the Dirichlet-toNeumann mapping (DNM) technique, the matched interface and boundary (MIB) method,
and the successive over relaxation (SOR)-like iterative procedure

38

First of all, the DNM is a computational technique using the Green’s function in order to
remove charge singularities due to the Dirac delta function δ (r − rk ) in the Poisson equation
(2.8) [85, 217, 221, 220]. It decomposes the electrostatic potential Φ(r) into the regular part
Φ(r) and the singular part Φ(r), which results in the following equations:
Ns

−∇ ·

(r)∇Φ(r) = (1 − χ(r))

qα Cα (r),

r ∈ Ω,

(2.20)

α=1

and




 ∇2 Φ(r) = 0, r ∈ Ω

m

(2.21)



 Φ(r) = −Φ∗ (r), r ∈ Γ.
This process generates a system of equations which can be numerically well-interpolated
without losing the order of accuracy. For more information, readers may refer to Section A.1.
The MIB scheme has been extensively developed to deal with discontinuous coefficients
and irregular complex geometries [208, 33, 30, 84, 85, 214, 217, 216, 219, 221, 220, 224, 226].
Basically, this method uses fictitious values and auxiliary points to extend the computational
domain so that the desired order of convergence is achieved [217, 216, 226, 219]. Appendix
A.2 outlines the principal ideas of the MIB method and introduces the required discretization
equations. The advantages of this numerical scheme are using standard finite difference
method on a simple rectangular coordinate system, carrying out low order of physical jump
conditions and locally reducing three dimensional interface problem to one dimensional
problem [221]. The standard centered finite difference scheme is implemented to find Φ0
from Eq (2.21) and Cα from Eq. (2.19) without non-electrostatic interaction term Uα away
from the complex boundary. However, careful interpolation is required near the boundary for
both equations [217, 220]. Especially, the MIB method is utilized to solve Eq. (2.20) for Φ

39

with two jump conditions (2.11) and (2.12). Since these equations are coupled, an iterative
procedure is required to obtain convergent results as explained in detail in Section A.3.

40

Chapter 3
Modeling and simulating a
mechanoelectrical transducer (MET)
channel in mammalian hair cells in
molecular level
In this chapter, we introduce a molecular level prototype of a MET channel and computational
studies to explore the mechanoelectrical gating mechanism of the MET channel in mammalian
hair cells. Our channel model consists of a realistic ion channel, a membrane and an additional
charged lip. Our theoretical analysis has been carried out with the standard PNP model
which is able to provide channel current under various experimental conditions, such as bulk
ion concentration, applied voltage, lip charge and lip position. In particular, we examine
the response of total current according to lip displacement. Our results are compared with
experimental data of rat hair bundle displacement and current relation given by Kennedy et
al. [123].
The rest of this chapter is organized as follows. In Section 3.1, the structure of our
molecular level prototype for the MET channel is presented. Since the molecular structure of
the MET channel is unknown, an alternative channel is mandatory to design the channel

41

model in microscopic detail. We use the Gramicidin A (GA) ion channel (PDB ID: 1MAG)
for the following reasons. First, the structure of GA is simple so that a blocker can be
added with 100% confidence in the blocker’s geometry. In practice, many other channels
have additional domains and complex composition, which may cause unwanted complication
in the result interpretation. Second, the GA is one of the most studied realistic channels
and its behaviors are well-known [220, 33, 208]. Consequently, the change in fluidic behavior
through this channel can be easily attributed to the blocker effect. A particle-like lip with an
adjustable charge, also called a blocker, is mounted at the mouth region of the GA channel to
enable the mechanoelectrical gating of the channel. Section 3.2 introduces the preparation of
the GA atomic structural data and explains a few important simulation parameters including
the dielectric constants and the diffusion coefficients. Section 3.3 is devoted to the numerical
study of our prototype for the MET channel. We validate our model by analyzing the
electrostatic potential and current under a variety of channel conditions via the second-order
PNP solver recently developed in our lab [220]. Among our computational results, open
probability distributions are compared with experimental data in the literature. This chapter
ends with a conclusion in Section 3.4.

3.1

A molecular level prototype of a MET channel

Howard and Hudspeth [104] proposed the original gating spring model, in which the elastic
component, the so-called gating spring, regulates opening and closing of the transducer
channel. Deflection of a hair bundle toward the tallest row of stereocilia increases the tension
of the spring which is attached to the hypothetical channel gate. As a result, the channel
attains the open state. It was assumed that the tip link acted as a gating spring in the

42

beginning. However, recent research findings suggest that the tip link may be too stiff to
express a spring effect and thus an elastic filament may be located at the bottom of the
channel [88, 122, 175]. It is still mysterious how the deflection of a hair bundle gives rise to
opening the MET channels. Several possibilities have been suggested, but none of them is
conclusive owing to the lack of molecular level structural information for the MET channel
and insufficient experimental evidence. Moreover, the channel opening may involve not only
a conformational change but also a complicated molecular interplay [159, 164].

Figure 3.1: An illustration of the MET channel model consisting of GA channel protein, the
artificial membrane, and a positive ion which presents the gating effect.

In this work, we design a molecular level prototype to simulate the MET channel gating
process in a realistic setting. The GA channel (PDB code: 1MAG), one of the most wellstudied ion channels, is employed to represent the channel structure in molecular level. A
positively charged atom of radius 1.5˚
A, called as a “blocker”, is placed at the mouth region
of the channel so that it behaves as a lip, see Fig. 3.1. Since the GA channel is a cation
channel, it can be effectively blocked by an additional cation with a pertinent charge. The
43

computational domain of the MET channel prototype incorporates four different regions such
as channel pore, channel protein, bulk regions and membrane as in Fig. 2.1(a).
The channel pore region is placed along the z-axis. The membrane part is added to the
geometry from z = −15˚
A to z = 9˚
A (See Fig. 2.1(a)). To demonstrate the blocker’s gating
behavior, we allow it to move along one direction just outside the channel entry. Then the
displacement of the blocker imitates the deflection of a hair bundle [123]. In order to find
out the optimal position of the blocker, several locations along the channel axis (z-axis)
were tested; as a result, zb = −15.95˚
A is selected to be the most suitable z-coordinate of
the blocker. We fix the y-coordinate of the blocker at 0.1˚
A and allow the blocker to move
only along the x-direction from xb = −0.5˚
A to xb = 0.9˚
A, which involves the transition
between the open and closed states of the channel model. Additionally, in order to intensify
the efficiency of the blocker role, we vary its charge over a certain range in our numerical
experiments.

3.2

Computational setup for the PNP systems

The standard PNP model, a relatively simple and computationally inexpensive continuum
model, is employed to examine our MET channel prototype, in order to optimize the
computational efficiency and to focus on the gating behavior of the prototype. Our earlier
work already validated that the standard PNP model simulated the GA channel at the second
order accuracy using a number of advanced mathematical techniques, and the simulation
prediction showed a good agreement with experimental result from the literature [220].
Further, the PNP model can incorporate with the structural detail of ion channel proteins. In
the present simulation study, the molecular structural data of the GA channel with a blocker

44

is prepared in the procedure as described in our earlier work [220]. Specifically, for each
atom in the biomolecule, its van der Waals radius is assigned with the CHARMM22 force
field [143] in order to achieve a full all-atom structure. Then partial charges for the channel
protein are calculated by the PDB2PQR software [64, 65], which are accounted in the fixed
charge density in the Poisson equation (2.8) in our computation. However, partial charges in
the membrane layers are neglected to attain a good approximation for ion transport. The
molecular surface of the GA channel is generated by the MSMS program with density 10
and probe radius 1.4˚
A [170] to identify the channel-solvent interface Γ which separates two
domains Ωm and Ωs .
We need to treat carefully parameters involved in the PNP equations. First, it is necessary
to assign pertinent values for the dielectric constant function (r) because the Poisson equation
(2.8) requires the spatial-dependent dielectric constant. In our computation, we set m = 2
for the biomolecular region impermeable to water or ions and s = 80 for the aqueous solution
region permeable to water and simple ions. Another significant parameter is the diffusion
coefficient for each ion species of interest in the solvent region. The diffusion coefficient in
bulk region is different from, usually larger than, that in channel pore region. Practically, the
overall diffusion coefficient distribution is intricate to be determined, but the bulk values for
relatively simple biomolecules are possibly available from experiments [221]. Furthermore,
the diffusion coefficient function Dα (r) should be deliberately defined to be differentiable
using the buffering zone between the pore and bulk regions as discussed in our earlier work
[220].


pore
pore


Dα ,
r ∈ Ωs




pore
pore
bulk f (r), r ∈ Ωbuffering
Dα (r) = Dα
+ Dα − Dα
s






 Dbulk ,
r ∈ Ωbulk
α
s .
45

(3.1)

Here, the function f (r) is defined as

f (r) = f (z) = n

n+1
z − z pore
− (n + 1)
z bulk − z pore

n
z − z pore
, z ∈ z bulk , z pore , (3.2)
z bulk − z pore

where z bulk and z pore , respectively, indicate the bulk region and the channel pore region.
In our computation, we use Dbulk
= 1.96 × 10−5 cm2 /s and Dbulk
= 2.06 × 10−5 cm2 /s.
K+
Cl−
pore

Additionally, the diffusion coefficient for the pore region is Dα

bulk /21 for each species
= Dα

α [220].

3.3

MET channel prototype simulations

In the present simulation work, we mainly study the electrostatic distribution and transport
features of potassium (K+ ) chloride (Cl− ) in our channel prototype. In order to verify our
prototype, the computational results are compared with the experimental findings from the
reference [123].

3.3.1

Experimental data

Experiments were conducted on outer hair cells of Sprague-Dawley rats. At first, the apical
and middle turns of the organ of Corti were excised from rats between 6 and 14 days postnatal
and then chemically treated as previously reported in the literature [17, 18, 123, 124]. The
separated turns were fixed in the experimental chamber and the axial motion of a glass pipette
deflects hair bundles, which activates the MET channels. Currents and bundle motions were
recorded mostly at the beginning of the apical turn and averages of 10 stimuli were calculated
to obtain a standard error of ±1 [123].

46

We mainly focus on the MET channel current which reveals the state of the channel.
The time course of the MET channel current can be divided into three steps [123]. The
first step represents that the channels respond to the bundle motions within microseconds,
which proves the direct relationship between bundle motion and channel activation. Then
the second part describes fast adaptation, which indicates the process of the channel closure.
The third part refers to slow adaptation of the MET channel mechanism, which is related to
the restoration of the sensitivity of the hair bundle.
The next significant feature of the MET channel is the normalized current (I/IMax ),
which can be regarded as the channel open probability. It is also presented in the literature
[123]. The relation between the normalized current and the displacement of the bundle shows
a nonlinear sigmoidal behavior, which demonstrates the following. The positive deflection
(toward the longest stereocilium) increases the open probability; on the contrary, the negative
deflection (away from the longest sterecilium) decreases the open probability.

3.3.2

Computational results

We explore the behavior of our MET prototype under various conditions including the applied
voltage, bulk ion concentration, and the charge and position of the blocker. One of the
fundamental issues is the electrostatic potential profile projected to the channel direction
(z-direction) at the middle of the channel pore. Except specified, we set the charge of the
blocker as Q = 2ec with an elementary charge ec , the bulk concentrations of both ion species
as C 0 = C 0 + = C 0 − = 0.1M (molar) and the applied external voltage at the electrode near
K
Cl
the channel gate as Φ0 = 0.1V. Note that the concentration for both ion species are always
identical at the bulk regions.
We first examine the impact of moving the charged blocker along the x-direction. The
47

Figure 3.2: Electrostatic potential profiles at different locations xb of the blocker with charge
Q = 2ec when C 0 = 0.1M and Φ0 = 0.1V. The region between two dashed lines indicates the channel
pore region. As the blocker position gets closer to 0.9˚
A along the x-direction, the blocker produces
a stronger barrier near the channel mouth as shown in the electrostatic potential profile.

goal of this study is to verify whether the proposed MET prototype suitably represents the
MET channel in terms of its mechanically sensitive behavior. We transfer the charged blocker
from xb = −0.5˚
A to 0.9˚
A with fixing its y- and z-coordinates. As illustrated in Fig. 3.2, at
locations xb = 0.6˚
A and xb = 0.9˚
A, the charged blocker creates a large electrostatic barrier
at the channel mouth region, which hinders the inward flow of the cation K+ . Consequently,
the electrostatic potential level inside the channel becomes very lower due to the fact that
the positive ions cannot effectively enter through the channel. Therefore, the proposed MET
prototype shows a desirable gating effect.
Next, we examine the behavior of our channel model in response to the change in bulk
ion concentration C 0 . We place the blocker with the charge of Q = 2ec at xb = 0.9˚
A. We
set Φ0 = 0.1V for the applied voltage. Two different bulk ion concentrations C 0 = 0.1M
and C 0 = 0.2M are tested in our numerical experiments. Fig. 3.3(a) illustrates the response

48

(a) Effect of bulk ion concentration

(b) Effect of applied external voltage

Figure 3.3: Electrostatic potential profiles along the channel direction for the MET prototype (a)
at different bulk ion concentrations C 0 = 0.1M (circles) and C 0 = 0.2M (diamonds) under the fixed
applied voltage Φ0 = 0.1V and (b) at different applied external voltages Φ0 = 0.1V (circles) and
Φ0 = 0.2V (diamonds) under the fixed bulk concentration C 0 = 0.1M. In both investigations, we
locate the blocker at xb = 0.9˚
A and set the blocker charge as Q = 2ec . Increase either in the bulk
concentration or in the external voltage creates higher potential within the channel.

of the present model at these two bulk ion concentrations by comparing the electrostatic
potential curve along the channel pore direction. The variance in the bulk ion concentration
has little influence on the electrostatic profile outside of the channel pore region. Especially,
the electrostatic property of the outside part is dominated by the applied voltage and also
reflects the electrical neutralization of cations and anions. However, the increase in the bulk
ion concentration leads a higher concentration of the positive ion to pass through the channel
pore region, which enhances the electrostatic potential within the channel as well.
We further investigate the behavior of our model under two different applied voltages. To
this end, we set xb = 0.9˚
A, Q = 2ec and C 0 = 0.1M. The external voltage applied at the
electrode near the channel mouth is doubled from Φ0 = 0.1V to Φ0 = 0.2V. The electrostatic
potential profiles in response to this change are depicted at Fig. 3.3(b). We see an obvious

49

change in the electrostatic potential in the left bulk region near the channel mouth and in
the left part of the channel pore region.

(a) Effect of bulk ion concentration

(b) Effect of applied external voltage

Figure 3.4: The total current versus the relative displacement dx of the charged blocker. (a)
Total current behavior under two bulk ion concentrations C 0 = 0.1M (circles) and C 0 = 0.2M
(diamonds) when Q = 2ec and Φ0 = 0.1V are fixed. (b) Total current behavior under two external
voltages Φ0 = 0.1V (circles) and Φ0 = 0.2V (diamonds) when Q = 2ec and C 0 = 0.1M are fixed.
Consequently, increase either in bulk ion concentration or in external voltage induces increase in
channel current.

To show the effectiveness of the blocker, we study the impact of the displacement of the
blocker along the x-direction under the fixed blocker charge Q = 2ec as shown in Figure 3.4.
We analyze the channel current pattern under two bulk ion concentrations C 0 = 0.1M and
C 0 = 0.2M as shown in Fig. 3.4(a) and under two applied voltages Φ0 = 0.1V and Φ0 = 0.2V
as shown in Fig. 3.4(b). In this numerical experiment, we define a relative displacement dx as
the scaled distance of the blocker from its optimal position. When the relative displacement
is zero, i.e., dx = 0, the blocker locates right on top of the channel mouth, which essentially
blocks the channel and creates the closed state. When dx = 1, there is no blocking effect
and hence the system reaches its normal open state. Clearly, the current across the channel

50

gets higher as the relative displacement increases. Meanwhile, the channel current increases
as the bulk ion concentration increases because of permeating more cations through the
channel, which is consistent with the increase in the electrostatic potential within the channel
pore region shown in Fig. 3.3(a). Similarly, the channel current increases as the applied
external voltage increases due to the increase in the electrostatic potential as shown in Fig.
3.3(b). Moreover, these results agree qualitatively with those of hair bundle deflections in the
literature [123, 124].

(a) Electrostatic potential

(b) Total current

Figure 3.5: (a) Electrostatic potential profiles and (b) total current profiles for four different
charges Q = 0.5ec (triangles), Q = 1ec (diamonds), Q = 1.5ec (squares) and Q = 2ec (circles) when
C 0 = 0.1M and Φ0 = 0.1V are fixed. The region between two dashed lines indicates the channel
pore region in (a). As the blocker charge gets increased, the barrier at the gate of the channel gets
higher. In conclusion, the amplitude of the current curve also gets enlarged.

Having established the importance of the blocker position to the channel open-closed
transition, we are interested in other properties of the blocker that may contribute to its
gating mechanism. Two most relevant properties are the blocker size and charge. In the
present system, the change of the molecular size of the blocker is ineffective in the channel
behavior because the standard PNP model does not consider the finite volume effect [114, 138].
51

However, it can be studied how the blocker charge affects the ion conductance through the
MET prototype. To this end, we test our model with a number of different atomic charges,
namely, Q = 0.5ec , Q = 1ec , Q = 1.5ec and Q = 2ec , while set C 0 = 0.1M and Φ0 = 0.1V.
Our calculation results are presented in Fig. 3.5. In Fig. 3.5(a), the electrostatic potential
at the entrance of the channel is dramatically influenced by the magnitude of the blocker
charge. A larger charge generates a higher potential barrier, which results in less density of
the cation and a lower electrostatic potential within the channel pore region. In Fig. 3.5(b),
we illustrate the impact of the blocker charge on the relation between the channel current and
the relative blocker displacements dx . It is interesting to notice that only when the atomic
charge reaches an appropriate threshold, the blocker can effectively close the channel gate.
Therefore, an ionic current leakage exists when the blocker does not carry sufficient charge or
is not located at its optimized position.
It is attractive to compare our model prediction with the experimental finding of the
relation between the open probability and the rat hair bundle displacement given in the
literature [123]. To this end, we compute the channel open probability PO at a number
of relative displacement dx with four different atomic charges when we set C 0 = 0.1M
and Φ0 = 0.1V. For an explicit comparison, our calculation data are presented in Fig. 3.6
together with experimental data [123]. Amazingly, there is an excellent agreement between
our theoretical prediction and experimental measurement when the blocker charge is Q = 2ec .
Thus we can conclude that Q = 2ec is the optimal charge in our molecular level prototype of
the MET channel. Our simulation result implies that our molecular level prototype is able to
reveal the mechanoelectrical gating mechanism of the MET channel in hair cells.
Finally, our optimal model predictions are obtained with fixed applied voltage and bulk
ion concentration. Obviously, it remains to investigate whether our model prediction is
52

Figure 3.6: Open probability PO is plotted against the relative displacement dx under four different
atomic charges Q = 0.5ec (triangles), Q = 1ec (diamonds), Q = 1.5ec (squares) and Q = 2ec (circles)
where C 0 = 0.1M and Φ0 = 0.1V are set. The solid dots are experimental data of the normalized
MET current versus normalized rat hair bundle displacement obtained from the literature [123]. At
each blocker charge, the open probability in response to the relative displacement forms a sigmoidal
shape and, especially, the charge Q = 2ec gives a remarkable agreement with the experimental
result.

sensitive to the variance of these two experimental conditions. As shown in Fig. 3.7, although
bulk ion concentration and applied external voltage are doubled in our numerical simulations,
the computational data does not change very much in terms of the open probability. There
is still an outstanding consistency between our MET model prediction and experimental
measurement [123]. This result further validates the robustness of our molecular level
prototype for the rat MET channel.

3.4

Concluding remarks

The auditory system is one of the most significant sensory systems for mammals. Mechanoelectrical transducer (MET) is a principal device in mammalian auditory system for the brain

53

(a) Effect of bulk ion concentration

(b) Effect of applied external voltage

Figure 3.7: Sensitivity test of the molecular level MET model for the prediction of open-closed
probability according to the relative blocker displacement when (a) the bulk ion concentration is
doubled; (b) the applied external voltage is doubled. Both cases are fairly consistent with the
experimental finding [123].

to perceive sounds. It is generally speculated that hair cell deflections regulate the MET
channels located at the stereocilia on each hair cell top. However, the molecular building
blocks of the MET channel are not yet available to date and hence its gating mechanism is
still evasive. In this work, we construct a molecular level prototype for the mammalian MET
channel in order to elucidate the mechanoelectrical gating mechanism of mechanotransduction
in mammals. Our MET prototype consists of a realistic ion channel, namely, the Gramicidin
A (GA) channel, an additional charged lip, called a blocker, which is positioned at the GA
channel mouth, and membrane layers.
To explore the physical properties of the proposed MET prototype, we employ a welltested theoretical model, the Poisson-Nernst-Planck (PNP) model, for three dimensional (3D)
numerical simulations of the MET channel transport. Advanced computational techniques
which have been well-established in our earlier work, such as Dirichlet-to-Neumann mapping

54

(DNM) and matched interface and boundary (MIB) method, are utilized for the present
numerical simulations. We design extensive numerical experiments to analyze in detail the
electrostatic potential and channel current in response to variations of blocker charge, blocker
relative displacement, bulk ion concentration and external voltage applied at the electrode in
the channel entrance region. Finally, we compare our prototype prediction with experimental
measurement from rat outer hair cells [123] in regard to the relation between channel open
probability and relative displacement. Remarkable consistency between two data is observed
under pertinent physical conditions. We further demonstrate that our molecular based model
is insensitive to the change of both ion concentration in bulk region and applied external
voltage. Our numerical findings indicate that hair-cell tip links, which connect hair bundles,
efficiently convey mechanical force to mechanosensitive transduction channels and manage
their opening and closing.

55

Chapter 4
Modeling and simulating three types
of ionic diffusive nanofluidic channels
In this chapter, we introduce a nanofluidic system model to analyze a realistic ionic diffusive
nanofluidic channel with atomic charge details and propose a second-order convergent numerical method for nanofluidic problems. First of all, we validate the second-order accuracy
of our proposed numerical solver for Poisson-Nernst-Planck (PNP) equations with ionic
diffusive nanochannels in three dimensional (3D) realistic settings. Next, we investigate
the impact of atomic charge distribution on the fluidic behavior of a few 3D nanoscale
channels. We demonstrate that atomic charges give rise to specific and efficient control of
nanochannel flows. Furthermore, the change of the distribution in atomic surface charges
is orchestrated with the variation of applied external voltage and bulk ion concentration in
order to understand nanofluidic currents. Therefore, we are able to elucidate quantitatively
the transport phenomena of three types of ionic diffusive nano-scaled channels, including a
negatively charged channel, a bipolar channel and a double-well channel. These flow patterns
are examined in terms of electrostatic potential profiles, ion concentration distributions and
current-voltage characteristics.
The rest of this chapter is organized as follows. In Section 4.1, we introduce an atomic
scale nanofluidic structure including a cylindrical channel with atomic charges and two bulk

56

reservoirs. Section 4.2 is devoted to validate the present PNP calculation with synthetic
nanoscale channels. We first test a cylindrical nanochannel with one charged atom at the
middle of the channel length and then examine the channel with eight atomic charges that
are equally placed around the channel. Since PNP equations admit no analytical solution
in general, we design analytical solutions for a modified PNP system, which has the same
mathematical characteristics as the standard PNP system. In Section 4.3, we investigate the
atomic scale control and regulation of cylindrical nanofluidic systems. Three ionic diffusive
nanofluidic channels, such as a negatively charged channel, a bipolar channel and a doublewell channel, are inspected in terms of electrostatic potential, ion concentration and current.
Finally, this paper ends with concluding remarks in Section 4.4.

4.1

Atomic scale design of an ionic diffusive nanofluidic
channel

(a) A 3D view

(b) A 2D cross-section view

Figure 4.1: Illustration of a nanofluidic system geometry. (a) A 3D view of a schematic cylindrical
nanochannel whose ends are connected to two reservoirs of KCl solution. (b) A 2D cross-section
view of the cylindrical channel whose diameter is 10˚
A and length is 49˚
A in the xz-plane. Here, Φ0L
and Φ0R , respectively, represent the external voltage applied at the left and right electrodes, and C 0
represents the bulk ion concentration of both K+ and Cl− in two reservoirs.

In the present work, we construct a 16 × 16 × 56 cuboid nanofluidic system. It contains an

57

ionic diffusive cylindrical nanochannel which is placed at the center of the system and whose
both ends are connected to reservoirs of potassium chloride (KCl) solution as illustrated
in Fig. 4.1(a). The radius of the channel pore is 5˚
A and the length of the channel is 49˚
A
as depicted in Fig. 4.1(b). In our simulation study, the computational domain Ω is the
nanofluidic system which is mainly divided into ion inclusion region Ωs and ion exclusion
region Ωm as in Fig. 2.1(b). The ion inclusion region contains the region inside the channel
and two reservoirs where ions may penetrate and travel through, while the ion exclusion
region is the rest where there is no mobile ion, but there are fixed charged particles. In
contrast to our differential geometry based multiscale models [207, 208, 206], the interface
Γ between two regions Ωs and Ωm is predetermined by the channel structure and does not
change during our simulation.

(a) A 2D cross-section view

(b) A 3D view

Figure 4.2: Illustration of atomic charge distribution. (a) A 2D cross-section view and (b) a 3D
view of schematic diagram consisting of the cylindrical channel and four atomic charges which are
equally placed around the nanochannel at z = 0˚
A.

A number of properly charged atoms are positioned around the channel at an interval
of about 1.8˚
A so that the channel flow can be regulated by electrical charges. In reality,
these charged atoms can be realized by appropriate dopants. The z-coordinates of the atoms
along the channel length are determined first. Then at each cross section perpendicular
58

to the channel axis, the atomic charges are aligned along a concentric circle whose size is
sufficiently bigger than that of the channel pore. The locations of the atoms are equally
spaced according to the circumference of the circle. Fig. 4.2 shows an example of placing four
negative atomic charges around the cylindrical channel at z = 0˚
A. In the cross section on the
xy-plane, we divide a concentric circle with radius 6.5˚
A into four parts and then locate each
anion as described in Fig. 4.2(a). Managing the number, magnitudes and signs of atomic
charges enables generating various types of surface charge distribution on the ionic diffusive
nanofluidic channel.

4.2

Numerical validation

In our earlier work [220], we develop and verify a second order of PNP solver for proteins
and biological ion channels. In this work, we construct analytically solvable solid-state
nanofluidic systems to validate the proposed numerical methods. The analytic solution of
the PNP equations is unknown for realistic geometries. However, it is a standard procedure
to design an analytical solution for slightly modified PNP equations, which share the same
mathematical features with the original PNP equations [220]. Consequently, the numerical
convergence of the designed solution algorithms can be validated.
First, we consider two simple nanochannel examples, one with a single atomic charge,
and the other with eight atomic charges, in order to verify the second order convergence of
our numerical methods. Next, both a negatively charged nanofluidic channel and a bipolar
nanofluidic channel are utilized to further validate the proposed numerical schemes.
In our computation, we match the nanofluidic system with a cuboid which has the size
of [−8, 8] × [−8, 8] × [−28, 28]. We set m = 1 and s = 80 in all the numerical tests in

59

this section. Therefore, there is a sharp discontinuity in the dielectric coefficient across the
channel-solvent interface.

4.2.1

Analytical solution system

We consider a set of Nm charged atoms, each of which is placed at rk = (xk , yk , zk ) in the
microscopic channel region Ωm and possesses a fixed charge Qk , where k = 1, 2, · · · , Nm .
The geometry of the analytically solvable system can be arbitrary in principle. However, one
can refer to the cylindrical geometry described in Fig. 4.1 and the computational domain
illustrated in Fig. 2.1(b). We define the electrostatic potential as follows:

Φ (r) =


Nm




 cos x cos y cos z +

k=1 m

Qk
(x − xk )2 + (y − yk )2 + (z − zk )2

,




0.4π


cos x cos y cos z,
3s

r ∈ Ωm
(4.1)
r ∈ Ωs .

For two mobile ion species, the concentrations are given by

C1 (r) = 0.2 cos x cos y cos z + 0.1 and C2 (r) = 0.1 cos x cos y cos z + 0.1

(4.2)

for r ∈ Ωs , but C1 (r) = C2 (r) = 0 for r ∈ Ωm because ions are able to move only in the
solution region Ωs .
This set of solutions (4.1) and (4.2) satisfies the following PNP-like equations

Nm




− ∇ · ( (r)∇Φ(r)) = 4π
Qk δ(r − rk ) + 4π [C1 (r) − C2 (r)] + R(r),




k=1


∇ · [∇C1 (r) + C1 (r)∇Φ(r)] = R1 (r),







 ∇ · [∇C2 (r) − C2 (r)∇Φ(r)] = R2 (r),

r∈Ω
r ∈ Ωs
r ∈ Ωs .

60

(4.3)

Herein, we set, respectively, the charge of two mobile ion species as q1 = 1 and q2 = −1, and
their diffusion coefficients as D1 (r) = D2 (r) = 1 for all r ∈ Ωs . Moreover, for every r ∈ Ωs ,




R(r) = −3 cos x cos y cos z





R1 (r) = −0.6 cos x cos y cos z + 0.4π
3 s ∇ · [(0.2 cos x cos y cos z + 0.1)∇(cos x cos y cos z)]






 R2 (r) = −0.3 cos x cos y cos z − 0.4π ∇ · [(0.1 cos x cos y cos z + 0.1)∇(cos x cos y cos z)] .
3 s
However, R(r) = R1 (r) = R2 (r) = 0 for all r ∈ Ωm . The jump conditions at the interface Γ
can be specifically given in the following:



 [Φ(r)] = Φ∗ (r) + 1 − 0.4π
3 s cos x cos y cos z


 [ (r)∇Φ(r)] · n = m ∇ [cos x cos y cos z + Φ∗ (r)] · n − s ∇ 0.4π cos x cos y cos z · n,
3 s
Nm

where n is the outward unit normal vector and Φ∗ (r)

Qk

=

(x − xk )2 + (y − yk )2 + (z − zk )2
In order to investigate the convergence order, we apply two error measurements
k=1 m

num − F exact | and L =
L∞ = max | Fi,j,k
2
i,j,k
i,j,k

1
N

num − F exact
Fi,j,k
i,j,k

2

,

i,j,k

num and F exact , respectively, represent the numerical and exact values of a function
where Fi,j,k
i,j,k

F at (xi , yj , zk ), and N is the total number of computational nodes.

4.2.2

A cylindrical nanochannel with a single atomic charge

We first test the cylindrical channel with a single atomic charge. The atomic charge is placed
at (6.5, 0, 0) and its charge is −0.08 as shown in Fig. 4.3(a). The set of analytical solutions
(4.1) and (4.2) of the PNP equations introduced in Subsection 4.2.1 is used to compare with
61

.

(a) With a single negative atomic charge

(b) With eight negative atomic charges

Figure 4.3: Illustration of the geometries of two simple numerical test cases.

numerical results by solving Eq. (4.3). We calculate the electrostatic potential Φ(r), the
positive ion concentration C1 (r) and the negative ion concentration C2 (r) at four different
mesh sizes: h = 0.4, h = 0.32, h = 0.2, and h = 0.16. Table 4.1 demonstrates numerical errors
and convergence orders for different number of computational nodes. The computational
result attains a good second order convergence. Moreover, the errors and orders for both
concentrations show little difference.

4.2.3

A cylindrical nanochannel with eight atomic charges

Next, we explore the same dimensional cylindrical nanochannel with eight fixed atoms with
a charge of −0.08 over the channel surface as illustrated in Fig. 4.3(b). Consider two
cross-sections perpendicular to the channel length at z = −11, and z = 11. In order to
maintain the same distance of each atom from the channel surface and the same angle
difference between two adjacent atoms, we employ the polar coordinate system. For each

62

Table 4.1: Numerical errors and orders for the cylindrical channel with a single atomic charge.

Φ

C1

C2

Mesh size
h = 0.4
h = 0.32
h = 0.2
h = 0.16
h = 0.4
h = 0.32
h = 0.2
h = 0.16
h = 0.4
h = 0.32
h = 0.2
h = 0.16

L∞
Error
1.3742E-2
8.7238E-3
3.5614E-3
2.2282E-3
3.8661E-3
2.5648E-3
8.4896E-4
5.4286E-4
2.0352E-3
1.3130E-3
4.7824E-4
2.7990E-4

Order
−
2.0362
1.9062
2.1016
−
1.8390
2.3524
2.0039
−
1.9642
2.1488
2.4327

L2
Error
3.1647E-3
2.0215E-3
8.7926E-4
5.3548E-4
1.0406E-3
6.4972E-4
2.5014E-4
1.5920E-4
5.3493E-4
3.2772E-4
1.3661E-4
8,0089E-5

Order
−
2.0088
1.7712
2.2224
−
2.1106
2.0309
2.0250
−
2.1958
1.8617
2.3930

Table 4.2: Numerical errors and orders for the cylindrical channel with eight atomic charges.

Φ

C1

C2

Mesh size
h = 0.4
h = 0.32
h = 0.2
h = 0.16
h = 0.4
h = 0.32
h = 0.2
h = 0.16
h = 0.4
h = 0.32
h = 0.2
h = 0.16

L∞
Error
1.3745E-2
8.7249E-3
3.5571E-3
2.2300E-3
3.8661E-3
2.5648E-3
8.4892E-4
5.4286E-4
2.0353E-3
1.3130E-3
4.7824E-4
2.7790E-4

63

L2
Order
−
2.0369
1.9090
2.0925
−
1.8390
2.3525
2.0037
−
1.9643
2.1488
2.4326

Error
3.1649E-3
2.0214E-3
8.7894E-4
5.3564E-4
1.0406E-3
6.4972E-4
2.5015E-4
1.5920E-4
5.3493E-4
3.2772E-4
1.3659E-4
8.0091E-5

Order
−
2.0092
1.7720
2.2194
−
2.1106
2.0308
2.0252
−
2.1958
1.8620
2.3923

atomic charge, the distance from the center of the channel pore is always 6.5˚
A and the
inter-atomic angle is a right angle. Therefore, the central coordinates of these atoms are
6.5 cos

π
π
π
π
i , 6.5 sin
i , −11 and 6.5 cos
i , 6.5 sin
i , 11 for i = 0, 1, 2, and
2
2
2
2

3. As a result, we obtain four equally spaced atomic charges at each cross section. We test
this nanofluidic channel system at four different mesh sizes: h = 0.4, h = 0.32, h = 0.2, and
h = 0.16. As shown in Table 4.2, the errors and orders in solving the PNP equations with this
atomic charge setting generate little difference from those with a single atomic charge. This
validation experiment also indicates the second order convergence of our proposed algorithms.

4.2.4

A negatively charged ionic diffusive nanofluidic channel

(a) A negatively charged channel

(b) A bipolar channel

Figure 4.4: Surface plot of electrostatic potential profiles on (a) a negatively charged channel and
(b) a bipolar channel. As represented in the color bar, blue colors represent negative values, while
red colors represent positive values. The negatively charged channel surface is mostly blue, but the
color of the bipolar channel surface is changed from red to blue.

Now, we perform another numerical test to verify the convergence and accuracy of the
proposed PNP calculation on a more realistic nanofluidic channel. At first, we design a
64

negatively charged nanochannel, or called a unipolar nanochannel. The channel length on
the z-axis is divided into 27 subdivisions and then there are 28 nodes between z = −23, and
z = 23. Thus the distance between two adjacent atoms along the z-direction is about 1.8. On
a circular cross section at each node along the z-axis, eight atoms with a charge of −0.08, all
of which are 1.5 apart from the cylinder surface They are equally spaced in a circular pattern
and so the interatomic angular difference is half a right angle. The x- and y-coordinates of
each atomic charge are calculated by using the polar coordinate system; consequently, The
detailed coordinates of the atomic charge distribution are given in Table 3.
In order to ensure that the channel surface is negatively charged, we first examine colored
surface plot and contour plots of the electrostatic potential Φ(r) by numerically solving Eq.
(4.3). The computational results are demonstrated in Fig. 4.4(a) and Fig. 4.5. The colored
surface plot is useful to understand the effect of the atomic charge on the nano-scaled channel.
Fig. 4.4(a) shows the electrostatic potential profile over the surface of this nanochannel. The
greater part of the channel surface has blue colors with a little different darkness. Since the
blue colors indicate the negative electrostatic potential values, it is obvious that the channel
surface possesses negative charge.
From the boundary line of the channel clearly shown in contour plots in Fig. 4.5, it
is interesting to discern that our proposed PNP solver works very well with ion diffusive
nanofluidic channels. Three contour plots of the electrostatic potential Φ(r) at z = −10,
z = 0 and z = 10 are described in Fig. 4.5. In every picture, we can notice that the right
outside of the channel is dark blue, which implies that the channel surface is negatively
charged.
Then we use the same analytical solutions (4.1) and (4.2) to calculate the numerical
errors and orders at four different mesh sizes: h = 0.4, h = 0.32, h = 0.2 and h = 0.16. The
65

(a) Cross section: z = −10

(b) Cross section: z = 0

(c) Cross section: z = 10

Figure 4.5: Three contour plots of electrostatic potential through the negatively charged channel
at z = −10, z = 0 and z = 10 on the xy-plane when the mesh size is h = 0.2. The blue colors just
outside the channel imply the negative surface charge.

calculation results are given in Table 4.3. Through this numerical observation, we confirm
that our proposed PNP numerical schemes achieve the second order accuracy in computing
the potential and ion concentrations for the nano-sized channel with negative atomic charges
on its surface.

4.2.5

A bipolar ionic diffusive nanofluidic channel

We also consider a bipolar nanofluidic channel which functions as a nanofluidic diode. It is a
nano-sized channel whose atomic coordinates are the same as those of the aforementioned
negatively charged nanochannel. However, their charges are altered from positive to negative
or vice versa at the middle of the channel axis [41, 56]. In our ionic diffusive bipolar channel,
the first half cylinder (from z = −23 to z = −0.8519) is affected by atoms with charge of
0.08 and the atomic charges on the other half (from z = 0.8519 to z = 23) are −0.08. Table
4 presents all of the detailed central positions and charges of the atomic charge distribution
of this bipolar nanofluidic channel.
The computational results of the electrostatic potential through the bipolar channel are

66

Table 4.3: Numerical errors and orders for the negatively charged channel.

Φ

C1

C2

Mesh size
h = 0.4
h = 0.32
h = 0.2
h = 0.16
h = 0.4
h = 0.32
h = 0.2
h = 0.16
h = 0.4
h = 0.32
h = 0.2
h = 0.16

(a) Cross section: z = −10

L∞
Error
1.3803E-2
8.7573E-3
3.5184E-3
2.3190E-3
3.8725E-3
2.5664E-3
8.6571E-4
5.4302E-4
2.0402E-3
1.3126E-3
4.8312E-4
2.7784E-4

L2
Order
−
2.0389
1.9401
1.8682
−
1.8437
2.3121
2.0902
−
1.9765
2.1265
2.4793

Error
3.1672E-3
2.0218E-3
8.7433E-4
5.3538E-4
1.0406E-3
6.4976E-4
2.5037E-4
1.5936E-4
5.3508E-4
3.2767E-4
1.3660E-4
8.0291E-5

(b) Cross section: z = 0

Order
−
2.0116
1.7836
2.1980
−
2.1103
2.0291
2.0244
−
2.1978
1.8615
2.3816

(c) Cross section: z = 10

Figure 4.6: Three contour plots of electrostatic potential through the bipolar channel at z = −10,
z = 0 and z = 10 on the xy-plane when the mesh size is h = 0.2. As indicated in the color bar, blue
colors represent negative values, while red colors represent some small positive values. At the cross
section of z = −10 the inside of the channel is blue (i.e., negatively charged solution), whereas at
the cross section of z = 10 the inside of the channel is red (i.e., slightly positively charged solution).

shown in Fig. 4.4(b) and Fig. 4.6 by solving the system of equations (4.3). From Fig. 4.4(b),
we are able to see that the atomic surface charges are changed from positive to negative
along the z-direction. Such properties of the bipolar channel are also clearly manifested in
the cross-sectional results in Fig. 4.6. When we compare the contour plots in Fig. 4.6(a)
and Fig. 4.6(c), it is a little bit difficult to distinguish the colors of the outside of the

67

Table 4.4: Numerical errors and orders for the bipolar channel.

Φ

C1

C2

Mesh size
h = 0.4
h = 0.32
h = 0.2
h = 0.16
h = 0.4
h = 0.32
h = 0.2
h = 0.16
h = 0.4
h = 0.32
h = 0.2
h = 0.16

L∞
Error
1.3752E-2
8.7568E-3
3.5649E-3
2.2421E-3
3.8743E-3
2.5663E-3
8.5125E-4
5.4306E-4
2.0384E-3
1.3132E-3
4.8477E-4
2.7790E-4

L2
Order
−
2.0227
1.9121
2.0782
−
1.8458
2.3479
2.0143
−
1.9706
2.1203
2.4935

Error
3.1649E-3
2.0217E-3
8.7942E-4
5.3585E-4
1.0406E-3
6.4975E-4
2.5016E-4
2.5016E-4
5.3493E-4
3.2773E-4
1.3665E-4
8.0280E-5

Order
−
2.0086
1.7711
2.2201
−
2.1105
2.0308
2.0191
−
2.1967
1.8611
2.3837

bipolar channel. However, the channel inside obviously shows different colors, which indicates
that the variation of atomic charges influences the ion transport through the ionic diffusive
nanochannel.
The validation of our numerical schemes for the bipolar nanofluidic channel is given in
Table 4.4. Again, we compare numerical solutions of Eq. (4.3) and exact solutions (4.1) and
(4.2) in this bipolar channel at four different mesh sizes: h = 0.4, h = 0.32, h = 0.2 and
h = 0.16. As we can see in Table 4.4, the computational result apparently demonstrates a
good second-order convergence for this test problem.

4.3

Nanofluidic simulations

Having validated the numerical convergence of our proposed PNP algorithm, we explore
the nature of charged ionic diffusive nanofluidic channels under various physical conditions
including applied voltage, atomic charge distribution and bulk ion concentration. Specifically,
we investigate the electrostatic distribution and transport features of potassium chloride
68

through the nanochannels. Our advanced numerical solver of the standard PNP theory
provides concentration for two ion species K+ and Cl− and electrostatic potential at each
node in the computational domain. Then the obtained values are averaged over the xy-cross
section at each z. Using such calculation process, we investigate ion concentration distributions
and electrostatic potential profile along the channel direction (z-direction). We also illustrate
current-voltage (I-V) curve in which the current at the center of the channel length is given
by Eq. (12). Particularly, ion concentration distribution describes the movements of two
different ion species through a channel in detail and current-voltage characteristic clearly
shows electrical features of a nanochannel. We examine three kinds of nano-scaled channels,
namely a negatively charged channel, a bipolar channel and a double-well channel. We have
used the dielectric constants m = 2 and s = 80 for all simulations in this section. For real
physical simulations of biological ion channels and synthetic nanofluidic channels with atomic
detail, it is also important to make use of one unit system [220]. Pertinent units that are
consistent with fundamental physical laws and central physical constants should be employed
so that the computational results can be directly compared with experimental measurements.
Thus, we discuss the units of the PNP model in this simulation work in Chapter B.

4.3.1

A negatively charged ionic diffusive nanofluidic channel

In order to clarify the role of atomic charge in ionic diffusive nanofluidic systems, we first
consider a negatively charged channel described in Subsection 4.2.4 and the atomic charges
of the channel are specified in Table 3. Here, all of the ions placed to generate the atomic
surface charge are anions with the same magnitude of charge. We vary the external voltage
only at the electrode of the right reservoir and keep bulk ion concentrations of both reservoirs
same. Therefore, the applied voltage is a driving force to relocate potassium ions and chloride
69

ions via the ionic diffusive nanochannel. The atomic charges determine the ion selectivity of
the nanochannel and the bulk ion concentration affects the current migrated through the
channel.

Effect of applied voltage

(a) Electrostatic potential

(b) Ion concentration

Figure 4.7: Effect of applied voltage on electrostatics and dynamics through the negatively charged
nanofluidic channel. (a) Electrostatic potential profiles and (b) ion concentration distributions along
the channel length (z−axis). Here ∆Φ is varied from 0V (squares), 0.2V (triangles), 0.4V (asterisks),
0.6V (diamonds), 0.8V (circles) to 1V (plus signs), where ∆Φ is the difference of the applied voltage
between two electrodes. Each surface atomic charge is Qk = −0.08ec and bulk concentration is
C 0 = 0.05M. Two dashed vertical lines indicate the ends of the nanochannel. As the applied voltage
difference gets larger, the potential at the right part of the inner channel is increased. Consequently,
more K+ ions (dashed line) are accumulated at the left part of the inner channel. However, there is
little change in the concentration of Cl− ion (solid line).

First, we study the impact of applied external voltage on behavior of ions traveling within
a negatively charged nanochannel. We set the ion concentrations of K+ and Cl− in both bulk
regions as C 0 = 0.05M and each atomic surface charge as Qk = −0.08ec . The voltage applied
at left electrode of the system, Φ0L , is fixed at 0V and the voltage applied at right electrode
of the system, Φ0R , is increased gradually from 0V to 1V. The ∆Φ represents the difference
70

between Φ0R and Φ0L , where ∆Φ = Φ0R − Φ0L . Figure 4.7 illustrates the electrostatic potential
and ion concentrations along the z-axis in response to the external voltage difference. In Fig.
4.7(a), as Φ0R gets increased, the electrostatic potential at the right part of the inner channel
becomes dramatically higher. As a result, more cations are accumulated at the left part of
the inner channel as shown in Fig. 4.7(b). In contrast, the concentration of Cl− shows little
difference from the bulk concentration. This transport phenomenon is originated from the
fact that the negative atomic charges of the channel electrostatically repel anions and attract
cations. This electrostatic effect also makes the solution within the channel almost unipolar,
which can be more obviously observed with current-voltage characteristic.

(b) Qk = −0.08ec

(a) Qk = 0ec

Figure 4.8: Ionic current for each ion species versus the applied potential difference ∆Φ through (a)
the uncharged nanochannel (Qk = 0ec ) and (b) the negatively charged nanochannel (Qk = −0.08ec ).
Here, the bulk ion concentrations at both reservoirs are fixed at C 0 = 0.05M. When Qk = 0ec , both
current-voltage graphs are linear and the positive current (squares) is roughly double of the negative
current (triangles). When Qk = −0.08ec , the positive current-voltage graph (squares) is nonlinear
and the negative current-voltage graph (triangles) is almost always zero.

By comparing the current-voltage (I-V) curve of the negatively charged channel with
that of an uncharged channel, we are able to clarify the effect of atomic surface charge
71

on a nanometer channel. In these two graphs, the current values are obtained using Eq.
(12) at the circular xy-cross section inside the channel when z = 0˚
A. The current at every
z-location through the channel is almost constant at steady state and hence the current at
the center of the channel axis is usually chosen to elucidate the current-voltage relation. Fig.
4.8(a) gives the relation between current of each ion species and external voltage within
the same dimensional nanofluidic channel without atomic charges (Qk = 0ec ). Both of the
ionic currents are proportional to the applied voltage and, especially, the K+ current is
roughly double of the Cl− current. In this case, the nanochannel obeys the Ohm’s law and is
non-selective. In contrast, the set of negative atomic charges on the channel surface destroys
the linearity of the positive current-voltage characteristic, and also generates a large difference
between two ionic currents as shown in Fig. 4.8(b). This nonlinearity or deviation from
the Ohm’s law is closely related to the non-proportionality between the potential change
within the channel and the applied voltage [57]. Since the negatively charged channel wall
hinders the access of Cl− ions at the channel gate, the negative current is almost zero for
every applied voltage. However, since the negative surface charge strongly attracts K+ ions
and the external potential difference encourages more ions to pass through the channel, the
cation current gets increased. From this numerical observation, we can conclude that an ionic
diffusive nanochannel with unipolar atomic charge distribution generates a charged current
mostly composing of counterions and the current through this channel can be increased by
providing more external voltage.

Effect of atomic charge
Next, we examine the effect of atomic surface charge on ionic flow through the negatively
charged channel. Except for the magnitude of Qk , we fix the bulk ion concentrations for
72

Figure 4.9: Effect of atomic surface charge on ionic current through the negatively charged
nanochannel. The ionic currents in response to the change of the atomic charge Qk are depicted.
Six different charges Qk = 0ec , Qk = −0.02ec , Qk = −0.04ec , Qk = −0.06ec , Qk = −0.08ec and
Qk = −0.1ec are simulated under the fixed bulk concentration C 0 = 0.05M and applied voltage
difference ∆Φ = 1V. Herein, |Qk | is the magnitude of the charge of the each atom. Stronger atomic
charges amplify the K+ current (squares) sharply, but reduce the Cl− current (triangles) to near
zero.

both ion species at C 0 = 0.05M and the applied voltage difference at ∆Φ = 1V. A stronger
negative atomic charge, i.e., a larger value of |Qk |, promotes the accumulation of the cations
inside the channel and prevents the anions from entering the channel gate. As a result, the
K+ current is dramatically increased; however, the Cl− current is decreased to almost zero.
The positive current is not directly proportional to the atomic charge magnitude because of
a stronger ionic diffusion induced by the larger concentration gradient [57]. From this study,
it is interesting to note that the channel current can be controlled by the atomic charge
amplitude. When the amplitude (∆I) of the current reaches a suitable threshold, almost all
ions, which have the same sign of charge with the channel atomic charge, cannot penetrate
through the inlet of the nanochannel. Therefore, our proposed nanofluidic channel possesses
a near perfect ion selectivity.

73

Effect of bulk ion concentration

Figure 4.10: Effect of bulk ion concentration on total current through the negatively charged
nanochannel. The graphs of total channel current against the external voltage difference (I-V)
are shown at five different bulk concentrations C 0 = 0.01M (squares), C 0 = 0.05M (triangles),
C 0 = 0.1M (diamonds), C 0 = 0.2M (circles) and C 0 = 0.4M (plus signs) when Qk = −0.08ec .
As the bulk ion concentration is increased, the total current gets higher and, moreover, the I-V
characteristic becomes near linear.

Another important aspect, which is also necessary to understand the transport within
an ionic diffusive nanofluidic channel, is the bulk ion concentration. The electrical double
layer produces a unique difference between a nanofluidic channel and a microfluidic channel.
The bulk ion concentration is a crucial factor to determine the thickness of the double layer.
When the double layers overlap inside a nanochannel, the aqueous solution confined in the
channel expresses charge opposite to the surface charge of the nanochannel [55]. Fig. 4.10
shows the total current as a function of the applied voltage difference for five different bulk
ion concentrations C 0 = 0.01M, C 0 = 0.05M, C 0 = 0.1M, C 0 = 0.2M and C 0 = 0.4M. Here,
the atomic charge Qk is assumed to be −0.08ec . The increased bulk ion concentration induces
more cations to penetrate through the channel, which results in the dramatic increase of the

74

total channel current. However, a higher bulk ion concentration reduces the double layer and,
additionally, the channel surface becomes neutralized by the pulled counterions within the
layer [57]. Consequently, the I-V curve is almost linear, in other words, obeys Ohm’s law. It
is pretty evident that a charged channel with high bulk ion concentration behaves like an
uncharged channel due to the decrease in the Debye length.

(a) The normalized ionic current

(b) The normalized channel current

Figure 4.11: The effect of bulk ion concentration on the normalized current through the negatively
charged nanochannel. (a) The normalized ionic current and (b) the normalized channel current with
respect to the increase in bulk ion concentration when the atomic charge Qk = −0.08ec and the
applied voltage difference ∆Φ = 1V are fixed. The negatively charged channel with a sufficiently
larger bulk ion concentration behaves like a uncharged channel because the normalized values
becomes one.

The normalized current is considered to assure that a higher bulk ion concentration
weakens the role of atomic surface charge. The normalized current is computed by dividing
the current through the negatively charged channel by that through the uncharged channel
so that we can observe the difference between two channels. Excluding the atomic charges,
both channels have the same voltage difference and the bulk ion concentration gets gradually
increased at the same time. We calculate the normalized value of each ion species current and
75

then add them together to obtain the total current. Fig. 4.11(a) represents the normalized
ionic current for each ion species and Fig. 4.11(b) shows the normalized total channel current.
In both figures, the normalized values get close to one as the bulk ion concentration is
increased. This feature implies that the charged channel under a high bulk ion concentration
does not demonstrate much atomic charge effect on the nanofluidic transport phenomena.

4.3.2

A bipolar ionic diffusive nanofluidic channel

As the second subject of the nanofluidic simulations, we examine a bipolar ionic diffusive
channel whose dimensions and atomic charge constitution are described in Subsection 4.2.5.
In this channel, the first half of the channel is positively charged and the second half is
negatively charged. Therefore, it is interesting to explore the transport properties of the
bipolar nanochannel because a bipolar channel can behave like a p-n junction.

Effect of applied voltage
We first consider three types of voltage bias across the channel length. One case is called
no bias (∆Φ = 0V) when both electrodes of the system have zero voltage. Another case is
named as a forward bias (∆Φ > 0) when the voltage applied at the right electrode of the
system is greater. The other case, on the contrary, is referred to as a reverse bias (∆Φ < 0)
when the voltage at the left electrode of the system is greater. To see the impact of external
voltage, we fix the bulk ion concentration of KCl at 0.1M, and the atomic charges on left
and right halves, respectively, at 0.08ec and −0.08ec . Fig. 4.12 compares the electrostatic
potential profiles and ion concentration distributions along the z-direction at three applied
voltage differences, including ∆Φ = 0V, ∆Φ = 1V and ∆Φ = −1V.
As shown in Fig. 4.12(i-a), when ∆Φ = 0V, the electrostatic potential is high under the
76

(i-a) Electrostatic potential

(i-b) Concentration

(ii-a) Electrostatic potential

(ii-b) Concentration

(iii-a) Electrostatic potential

(iii-b) Concentration

Figure 4.12: (a) Electrostatic and (b) ion concentration profiles of the bipolar nanofluidic channel
along the channel axis at (i) no bias ∆Φ = 0V, (ii) forward bias ∆Φ = 1V and (iii) reverse bias
∆Φ = −1V. Fix the bulk ion concentration C 0 = 0.1M and the atomic surface charges |Qk | = 0.08ec .

77

positive atomic charges, but it is low under the negative atomic charges. Subsequently, the
ion concentration is plotted in the opposite way as illustrated in Fig. 4.12(i-b). Generating
the potential gap ∆Φ between two ends of the system brings about two peculiar phenomena
within the bipolar ionic diffusive nanochannel. At the forward bias with ∆Φ = 1V, the
electrostatic potential is gradually increased as in Fig. 4.12(ii-a), but at the reverse bias with
∆Φ = −1V, it is sharply decreased as in Fig. 4.12(iii-a). These two results demonstrate
excellent consistency with the ion concentration curves in the way that the flux is invariable
along the channel axis at steady state and thus the main factor altering the potential is the
ion distribution [56]. As plotted in Fig. 4.12(ii-b), under the forward bias, both ion species
are attracted to the junction, so the peak value of the ion concentration is greater than that
under no bias. However, under the reverse bias, both ion species are moved away from the
junction and each one produces a small pile at the atomic surface charge with opposite sign as
presented in Fig. 4.12(iii-b). Accordingly, the forward bias brings about an ion accumulation
zone at the channel junction, whereas the reverse bias creates an ion depletion zone there.
These concentration patterns influence on current-voltage relation of the bipolar nanofluidic
channel.
Figure 4.13 shows the current-voltage relation of each ion species through the bipolar
channel. While both ionic current graphs are almost zero at every reverse bias (∆Φ < 0), they
are monotonically increased under forward bias as the potential difference becomes bigger.
This simulation phenomenon is closely related to the two novel transport phenomena within
a bipolar nanochannel: the ion-depletion zone under reverse bias and the ion-accumulation
zone under forward bias. We already manifest these features in Fig. 4.12(ii-b) and (iii-b),
respectively. The ion-depletion zone terminates the flow inside the bipolar nanochannel,
but the ion-accumulation zone encourages more ions to penetrate through the nanochannel.
78

Figure 4.13: Ionic current versus applied voltage difference in the bipolar nanofluidic channel.
The atomic charge |Qk | = 0.08ec and the bulk ion concentration C 0 = 0.1M are assumed. Under
reverse bias, both ion species cannot pass through the channel. However, under forward bias, the
currents of both ion species get enlarged. Especially, the cation current (squares) increases faster.

In particular, the positive ionic current enhances more abruptly as in Fig. 4.13. Another
remarkable discovery from the current-voltage characteristic is that a higher applied voltage
reduces the gradient of the curve

dI
. This feature is relevant to the polarization
d(∆Φ)

caused by the concentration gradient, which is analytically discussed in the literature [56].
To sum up, managing the external voltage facilitates turning on and off the current through
the bipolar nanochannel.

Effect of atomic charge
Next, we consider how atomic surface charge affects the dynamics inside the nano-scaled
channel. Fig. 4.14 depicts the channel currents in response to the applied voltage difference
at three different surface charges |Qk | = 0.04ec , |Qk | = 0.08ec and |Qk | = 0.12ec . In this
numerical experiment, we set the bulk ion concentration at two reservoirs to be 0.1M. In every
case of atomic charge, the total current gets enlarged as ∆Φ gets larger under forward bias,
79

Figure 4.14: Effect of atomic surface charge on the current through the bipolar nanofluidic
channel. Three sets of atomic charges, that is, |Qk | = 0.04ec (squares), |Qk | = 0.08ec (triangles)
and |Qk | = 0.12ec (diamonds), are studied. Here, the bulk ion concentration C 0 is fixed at 0.1M.
All I-V curves increase when ∆Φ varies from −1V to 1V. Greater magnitude of the atomic charge
results in a higher channel current with forward bias, but insufficient atomic charge may generate a
current leakage because it weakens the depletion zone with reverse bias.

but the rate of change of the current with respect to the voltage difference gets moderated.
The highest atomic charge, that is, |Qk | = 0.12ec , generates the greatest amplitude of the
current curve. In contrast, the lowest atomic charge amplitude, that is, |Qk | = 0.04ec does
not fully draw both ion species to either sides under reverse bias, and subsequently diminishes
the ion-depletion zone at the middle of the channel. Thus, the channel current becomes
nonzero at some negative voltage differences. From this numerical experiment, it is beneficial
to discover that the channel current within a bipolar ionic diffusive nanofluidic channel can
be perfectly regulated if the atomic charges on the channel surface satisfy an appropriate
threshold.

80

Figure 4.15: Effect of bulk ion concentration on the channel current through the bipolar nanofluidic
channel. Three sets of bulk ion concentrations, such as C 0 = 0.05M (squares), C 0 = 0.1M (triangles)
and C 0 = 0.2M (diamonds), are explored in terms of the current-voltage relation. The atomic
charges are given by |Qk | = 0.08ec . As the bulk ion concentration gets higher, the amplitude and
gradient of the current-voltage relationship gets maximized.

Effect of bulk ion concentration
We also test our bipolar ionic diffusive channel at two more different bulk ion concentrations,
namely, C 0 = 0.05M and C 0 = 0.2M, in order to observe the impact of bulk ion concentration
on the transport feature. Then we compare the total current-voltage graphs with that at
C 0 = 0.1M as described in Fig. 4.15. Every current-voltage curve nearly vanishes when ∆Φ
is negative, but significantly increases when ∆Φ is positive. As the bulk ion concentration is
doubled, more ions are accumulated at the junction of the bipolar nanofluidic channel under
forward bias, so the total current gets almost doubled. Moreover, the amplitude (∆I) and
gradient of the current alteration with respect to the external voltage difference

dI
d (∆Φ)

are increased as well. To this end, it is expected to surmise that bulk ion concentration
promotes the quantity of the total current through the bipolar nanochannel. Amazingly,
our computational outcomes are in a good agreement with other numerical studies in the
81

literature [56].

4.3.3

A double-well ionic diffusive nanofluidic channel

Finally, we investigate the transport phenomena within a double-well nanofluidic channel,
which is named after the shape of the electrostatic potential curve across the channel length.
The electrostatic potential within a cylindrical nanofluidic channel may have several potential
wells by modifying atomic charge distribution on the channel surface. Actually, the irregular
partial charges of biological ion channel proteins induce electrostatic potential wells. For
instance, Gramicidin A channel, one of the most well-known biological channels, is a doublewell transmembrane ion channel [220]. In this numerical observation, we design a new
cylindrical channel with the same dimensions by varying atomic surface charges so the the
electrostatic potential profile has a double-well structure. As illustrated in Fig. 4.16, only the
middle section (from z = −7.67˚
A to z = 7.67˚
A) of the channel axis is positively charged, but
the other parts of the channel are negatively charged. The detailed atomic charge distribution
of this channel is given in Table 5.

Figure 4.16: The schematic diagram of a three-dimensional double-well nanofluidic channel. The
channel length is divided into three parts. The first and last parts are negatively charged, but the
middle part is positively charged. The red dots indicate atoms with negative charge Qk = −0.12ec
and the green dots indicate atoms with positive charge Qk = 0.04ec .

82

Effect of applied voltage

(a) Electrostatic potential

(b) Ion concentration

Figure 4.17: Effect of applied voltage on electrostatics and dynamics through the double-well
nanofluidic channel. (a) Electrostatic potential profiles and (b) ion concentration distributions
along the channel length (z−axis). Here ∆Φ is varied from 0V (squares), −0.2V (triangles), −0.4V
(asterisks), −0.6V (diamonds), −0.8V (circles) to −1V (plus signs), where ∆Φ is the difference of
the applied voltage between two ends of the system. The bulk concentration is C 0 = 0.05M for
both ions K+ and Cl− . Two dashed vertical lines indicate the ends of the cylindrical channel. Each
electrostatic potential graph shows two potential wells, which brings about two piles of K+ ions
along the channel axis. Moreover, higher applied voltage at the left end of the system, Φ0L , breaks
the symmetry of the potential wells. The left well becomes weaker and the right well becomes
stronger.

At first, we alter only the applied voltage, but fix bulk ion concentration at C 0 = 0.05M. In
this study, Φ0R is set to be 0V and Φ0L is increased gradually from 0V to 1V. Fig. 4.17 presents
the electrostatic potential profiles and ionic concentration distributions along the channel
length at five different external voltages, that is, ∆Φ = 0V, ∆Φ = −0.2V, ∆Φ = −0.4V,
∆Φ = −0.6V, ∆Φ = −0.8V and ∆Φ = −1V. With no external voltage bias (∆Φ = 0V), the
electrostatic potential has two symmetric wells, which brings about two symmetric piles of
the K+ ions. The increase of the voltage applied at the left end of the system leads for the
electrostatic potential to get dramatically higher on the left hand side of the inner channel.
83

As a result, the left potential well is moderated as in Fig. 4.17(a). Parallel to this variation,
the positive ion concentration shows a dramatic change on the left hand side (See Figure
4.17(b)). On the contrary, the light change in the potential at the right hand side of the
channel inside corresponds to the light change in the concentration profile of cation on the
right. The negative atomic charges located around the channel gate electrostatically inhibits
the entering of anions, so the concentration distribution of Cl− is almost always zero as
shown in Fig. 4.17(b).

Effect of bulk ion concentration

Figure 4.18: Effect of bulk ion concentration on the channel current through the double-well
nanofluidic channel. The graphs of total channel current against the external voltage difference (I-V)
are studied at four different bulk concentrations: C 0 = 0.01M (squares), C 0 = 0.05M (triangles),
C 0 = 0.1M (diamonds), and C 0 = 0.2M (circles). Higher bulk ion concentration encourages the
total current to increase and the I-V characteristic to become linear.

We preform another numerical test with our designed double-well channel. In this test,
we only change bulk ion concentration, which is already proved to be an important factor to
determine the ion conductance with previous two ion diffusive nanofluidic channels. Herein,
84

we compare the current-voltage curves at four different bulk ion concentrations, including
C 0 = 0.01M, C 0 = 0.05M, C 0 = 0.1M and C 0 = 0.2M. As in Fig. 4.18, the increase in the
total current through the double-well nanochannel is derived from the increase in the bulk
ion concentration because more ions are passes from the reservoir to the channel inside. Like
the negatively charged channel, the I-V relation becomes linear as the bulk ion concentration
gets multiplied because higher bulk ion concentration weakens the double layer within the
channel.
The numerical observations of the electrostatic potential, concentration and current
through our double-well ionic diffusive nanochannel are qualitatively consistent with those
observed from both numerical simulations [220] and experimental measurements [25] of the
Gramicidin A channel. Therefore, the present study manifests a great possibility of the
atomic design of three dimensional nanofluidic channels in that the synthetic nanofluidic
channels with atomic surface charge distribution can be used to study biological channels,
which is particularly valuable when the structure is not available.

4.4

Concluding remarks

Recently the dynamics and transport of synthetic nanofluidic channels have received great
attention. Parallel to this trend, related experimental techniques and theoretical methods
have been substantially promoted in the past two decades [14, 15, 202]. Nanofluidic channels
are utilized for a vast variety of scientific and engineering applications, including separation,
detection, analysis and synthesis of not only chemicals but also biomolecules. Additionally,
inorganic nanochannels are manufactured to imitate biological channels, which is of great
significance in elucidating ion selectivity and ion current controllability in response to an

85

applied field in membrane channels [153, 185]. Molecular and atomic mechanisms are the
key ingredients in the design and fabrication of nanofluidic channels. However, atomic
details are scarcely considered in nanofluidic modeling and simulation, to our knowledge.
Previous simulations of transport in nanofluidic channels have been rarely carried out with
three-dimensional (3D) realistic physical geometry.
Present work introduces atomistic design of 3D realistic ionic diffusive nanofluidic channels.
The proposed mathematical model and numerical methods are employed for 3D realistic
simulations of such nanofluidic systems. Three distinct nanofluidic channels, including a
negatively charged nanochannel, a bipolar nanochannel and a double-well nanochannel, are
constructed to explore the capability and impact of atomic charges near the channel interface
on the channel fluid flow. We design a cylindrical nanofluidic channel of 49˚
A in length and
10˚
A in diameter. Several charged atoms of about 1.8˚
A gap are equally located at 1.5˚
A
from the channel surface in order to regulate nanofluidic patterns. For a negatively charged
channel, all of the atoms on the channel surface have the negative sign with the same charge,
whereas for a bipolar channel, half of them has the negative sign and the other half has the
positive sign, except that all the charges have the same magnitude. Finally, a double-well
channel has positively charged atoms at the middle and negatively charged atoms on the
remaining part of the channel surface. Each end of the channel is connected to a reservoir
of KCl solution and both reservoirs have the same bulk ion concentration. Asymmetry
in the applied electrostatic potentials at the electrodes gives rise to current through these
three nanochannels. We perform a number of numerical experiments to explore electrostatic
potential, ion concentration and current through each type of ionic diffusive channel under
the influence of applied voltage, atomic charge and bulk ion concentration.
The negatively charged channel generates a unipolar charged current because the negative
86

atomic surface charge attracts counterions, but repels coions. The current measure within
this nano-sized channel increases when external voltage, magnitude of atomic charges on the
channel surface and/or bulk ion concentration increase. However, the bulk ion concentration
has a limitation in its growth. A larger bulk ion concentration shortens the Debye length;
thus, the charged channel may behave like an uncharged one showing the Ohm’s law in the
current-voltage relation. The bipolar channel can create accumulation or depletion of both
ions in response to the current direction. When the right electrode has a higher voltage,
namely, with forward bias, both ions are stored at the junction of the channel length. On
the contrary, when the left electrode has a higher voltage, that is, with reverse bias, both
ions are moved away from the junction. External voltage applied at the ends of the system,
atomic surface charge and bulk ion concentration affect the amplitude and gradient of the
current-voltage characteristic. At last, the special atomic charge distribution of the doublewell channel produces the electrostatic potential profile with two potential wells. Increasing
applied voltage at the left hand side of the system results in an obvious change both in the
left potential well and in the left part of the K+ concentration.
The present study concludes that properties and quantity of the current though an ionic
diffusive nanochannel can be effectively manipulated by carefully altering applied voltage,
atomic surface charge and/or bulk ion concentration. Our computational observations
are compared well with those of experimental measurements and theoretical analysis in
the literature. Since the physical size of our channel model is close to that of realistic
transmembrane channels, the proposed nanofluidic model can be utilized not only for ionic
diffusive nanofluidic design and simulations, but also for the prediction of membrane channel
properties when the structure of the channel protein is not available or changed due to the
mutation.
87

Chapter 5
Thesis achievements and future work
The primary contributions of this dissertation are the mathematical modeling and simulations
of a mechanoelectrical transducer (MET) channel prototype and three types of ionic diffusive
nanofluidic channels. Both channel models are designed at the microscopic scale and,
especially, the atomic charges are employed to explain the gating or surface charge effect.
The present thesis study mainly focuses on the electrostatic and transport properties of these
channels under several physical conditions, such as applied potential, bulk ion concentration
and properties of an atomic charge. In order to explore the dynamics within both channel
models, we conduct various three dimensional numerical experiments using an electrodiffusion
continuum model, Poisson-Nernst-Planck (PNP) theory. Since these two channel models
have similar physical features, they have common mathematical interests and computational
challenges. At first, we propose the generalized Poisson-Nernst-Planck model to include nonelectrostatic interactions, which are significant factors to determine the transport phenomena
within a small nanometer-sized channel. Furthermore, several excellent computational skills
are employed to overcome singularities and interface problems.
Most of the materials of this dissertation work are based on the following publications:
• Jinkyoung Park, Kelin Xia and Guo-Wei Wei, “Atomic scale design and threedimensional simulation of ionic diffusive nanofluidic channels”, submitted to Microfluidics and Nanofluidics

88

• Jinkyoung Park and Guo-Wei Wei, “A molecular level prototype for mechanoelectrical
transducer in mammalian hair cells, Journal of Computational Neuroscience, 35, 231-241,
(2013)

5.1
5.1.1

Contributions
On the modeling and simulating a MET channel in mammalian hair cells

In ion channel modeling, one of the most difficult, but vigorous, challenges is to obtain
structural information about a target channel in molecular detail. The molecular composition
of a biological ion channel is essential to comprehend its major functions including ion selectivity and gating mechanism. In particular, the computational study of the mechanosensitive
channel in hair cell has several constraints to elucidate the intrinsic principles of mechanical
gating mechanism due to the lack of its molecular building blocks. In this respect, a molecular
level channel prototype, which we propose in the present work, contributes to the studies on
the MET channel in mammalian hair cells with the following aspects.
First of all, it is the first time to design a MET channel prototype in atomic detail, to our
knowledge. In our MET prototype, we incorporate an appropriately charged atom, called
a blocker, into the molecular structure of Gramicidin A (GA) channel. The GA channel
is very suitable to employ as the channel part of our MET prototype because of its wellestablished molecular building blocks. Additionally, its structure is comparatively simple and
safe to add an ion without worrying about any unexpected complication. Herein, the blocker
phenomenologically plays a role of the channel gate because its positive charge prevents the

89

inward flow of positive ions through this cation channel. The blocker displacement expresses
hair bundle deflection, which is believed to activate the MET channel. This molecular-level
MET channel model gives a great opportunity to account for the gating mechanism closely.
Secondly, we perform realistic three dimensional simulations of the MET channel transport
using the PNP theory. Our group has developed novel mathematical algorithms and useful
simulation skills, including Dirichlet to Neumann mapping and the matched interface and
boundary method, in order to investigate several proteins and ion channels in the past. In this
study, we employ the already proven PNP solver to demonstrate the electrostatic potential
profile and current-displacement characteristic though our MET channel prototype. These
simulation results in three dimensions can explain the gating principles of the MET channel.
Third, our numerical findings confirm that tip links, which connect hair bundles, are
involved in opening and closing the mechanosensitive transduction channels in hair cells. We
conduct a number of numerical experiments to analyze electrostatic potential and current
under the influence of an external voltage, blocker charge and displacement, and bulk ion
concentration. These computational results validate the usefulness and robustness of our
proposed MET channel model. Critically, when we compare our results with experimental
measurements in the literature, there is a remarkable consistency in terms of the relation
between channel open probability and relative displacement. This supports the idea that the
tip link is directly connected to the activation of the channel.

90

5.1.2

On the modeling and simulating charged nanofluidic channels

Nanofluidic channels are studied for diverse scientific and engineering purposes, including
separation, detection, analysis and synthesis of chemicals and biomolecules. Although
molecular and atomic mechanisms are crucial in designing nanofluidic channels, these details
are rarely explained in modeling and simulations to date. In addition, most transport
simulations in nanofluidic systems are carried out in one or two dimensional settings. Thus,
the second subject of the present study contributes to theoretical nanofluidic studies with
the following aspects.
First of all, we introduce the atomistic design of a three dimensional realistic ionic diffusive
nanofluidic channel. In our ionic diffusive nanofluidic system, atomic charges are utilized to
generate the surface charge of a channel. A number of charged atoms with the same size are
equally distributed around the channel structure boundary. One of the greatest strengths of
our nanochannel design is the simplicity to create a variety of surface charge distribution.
Since individual ions are used to express the fixed charges in the channel structure, we only
need to vary the sign and magnitude of each atomic charge. This design method enables
treating the surface charge effect more delicately.
Secondly, we propose a variational multiscale framework to facilitate the microscopic
description of nanochannel structures, including atomic surface charge, and the macroscopic
continuum treatment of the solvent and mobile ions. In this work, we generalize the PoissonNernst-Plank model for solid-state nanofluidic channel simulations. We employ a characteristic
function S to identify the channel-solution interface using the predetermined nanofluidic
channel structure. Then we incorporate the non-electrostatic interaction term U into the
91

total energy functional so that the generalized PNP theory is able to simulate very narrow
nanochannel. Further, we demonstrate the second order accuracy of our numerical algorithms
to solve the standard PNP equations with synthetic ionic diffusive nanochannels, including
the negatively charged channel and the bipolar channel.
Finally, we validate the usefulness and great feasibility of our channel model by simulating
three types of nanofluidic channels, namely, a negatively charged channel, a bipolar channel
and a double-well channel. Our simulation results, including electrostatic potential profile, ion
concentration distribution and current-voltage relation through these ionic diffusive channels,
show a great agreement with other observations in the literature. These features are very
helpful not only to understand the physical properties, but also to predict the functions for
each channel type. Especially, the transport behavior through the double-well nanofluidic
channel resembles that through a realistic biological ion channel, Gramicidin A channel. This
finding opens a great possibility to simulate ion channels when their structures are unknown.

5.2

Future work

Although the Poisson-Nernst-Planck theory is the one of the most widely used microscopic
approaches to describe electrokinetic processes in both biological ion channels and artificial
nanofluidic channels, it has several limitations. Especially, it ignores non-electrostatic
interactions between ions, ion-water and ion-protein, and steric effects which are critical to
determine the transport patterns due to a strong confinement or a dilute solution [12, 26,
52, 114]. As much as we understand these inter-particle interactions and atomic physical
features, we are able to design a more reasonable mathematical theory and to simulate ion
transport through a small pore more quantitatively.

92

For our future work, modeling and simulation can be substantially developed by incorporating the principles of the following well-established frameworks in the literature. For
example, Chen and his colleagues set up the free energy functional of short-range interactions
between proton-proton, proton-other mobile ions, proton-water and proton-protein in order to
describe the generalized correlation in the proton transport [33]. Each interaction is expressed
as a product of the number densities of related particles and the corresponding interaction
kernel. Similarly, we can incarnate non-electrostatic interactions between ions, ion-water and
ion-protein. As another fruitful computational framework, the density functional theory and
its variations have been employed to describe not only equilibrium, but also time-dependent
microscopic dynamics [211]. Especially, the time-dependent density functional successfully
elucidates the steric effects and intermolecular correlations [111]. Further, its combination
with the PNP theory describes well the electrolyte dynamics through a dilute solution [128],
and the chemical potentials of ions and the size selectivity of biological ion channels [87, 167].
By including this powerful mathematical insight to our dynamics model, we can simulate
more realistic nanofluidic systems.

93

APPENDICES

94

A
A.1

Mathematical algorithms
Dirichlet to Neumann mapping

Numerical treatment of the Dirac delta function in the Poisson equation may reduce the
order of accuracy as well as restrict the mesh size [85]. The Green’s function formulation for
the singular charges has been adopted for a higher-order interface scheme in the literature
[42, 217, 221, 220, 227]. First, the electrostatic potential Φ(r) is divided into the regular part
Φ(r) and the singular part Φ(r) so that the delta function is eliminated from Eq. (2.8). Here,
the singular part Φ(r) is composed of an essential solution of the Poisson equation Φ∗ (r)
and a harmonic function Φ(r) [85]. The singular part is only defined in the ion-exclusion
region Ωm because the fixed charges occur only in this region [221]. Specifically, the Green’s
function is
Φ∗ (r)

Nm

=
k=1

Qk
m |r − rk |

(1)

and Φ(r) satisfies the Laplace equation with a Dirichlet boundary condition


 ∇2 Φ(r) = 0,

r ∈ Ωm

(2)


 Φ(r) = −Φ∗ (r), r ∈ Γ.

Finally, this decomposition gives the Poisson equation of the regular part Φ (r) without charge
singularities which is
Ns

−∇ ·

(r)∇Φ(r) = 4π(1 − χ(r))

qα Cα (r),
α=1

where a factor of 4π is due to the selection of unit [220].
95

r ∈ Ω,

(3)

Furthermore, the following two jump conditions at the solvent-channel interface Γ should
be enforced in order to ensure the well-posedness of the problem [85].

[Φ(r)] = 0

(4)

[ (r)∇Φ(r)] · n = m ∇ Φ∗ (r) + Φ(r) · n,

(5)

where n is the outward normal direction from the interface.

A.2

Matched interface and boundary (MIB) method

The resulting equations (2) and (3) from the DNM algorithm involve material interface Γ
and three dimensional irregular geometry. Additionally, their solutions are still non-trivial if
the jump conditions (4) and (5) are to be properly enforced for arbitrarily complex interface
geometries. If the jump conditions are ignored, discontinuous coefficients and complex
geometry may reduce the order of convergence of the numerical scheme to at most one
[220]. Our group has proposed and developed an excellent numerical scheme, the matched
and interface boundary (MIB) method, to treat not only biological ion channels but also
nanofluidic systems with a higher order of accuracy [34, 32, 220]. In this section, we will
briefly elucidate the fundamental principles about the MIB method.
At first, every node in a computation domain is classified as either an irregular point
or a regular point. A point near the interface or boundary is called an irregular point, but
otherwise, it is called a regular point. The standard second order central finite difference
method is enough to interpolate at regular points. However, the application of the seven-point
finite difference scheme at an irregular point may result in divergence or reduce the order of
convergence due to the discontinuity of the functions involved. Thus, we need to carefully
96

deal with these points to obtain the convergence at the desired order of accuracy, especially
the second order of accuracy, in this dissertation work. If a discretization of Φ in Eq. (2)
involves any value from the region Ωs , it is replaced by the function value at the intersecting
node between the interface Γ and a Cartesian grid line. The difference scheme of Eq. (2.19)
is precisely modified according to the local geometry and the type of boundary condition and,
additionally, the detailed equations are discussed in the references [217, 220]. To calculate
Φ(r) numerically at irregular points from the interface problem (3), the MIB algorithm
employs fictitious values and auxiliary points.

Figure 1: An illustration of the MIB method for an irregular geometry. The kth grid line
intersects the interface at (xi , yo , zk ) (star), which involves two irregular points (circle) (xi , yj , zk )
and (xi , yj+1 , zk ). Here (xi , yo , zk+2 ) and (xi , yo , zk+1 ) (square) are two auxiliary points needed to
discretize ( Φz )z and (xi−1 , yo , zk ) and (xi−2 , yo , zk ) (square) are two auxiliary points needed to
discretize ( Φx )x . The other points (triangle) indicate the points selected to interpolate Φ along the
y-grid lines at (xi , yo , zk ), (xi , yo , zk+1 ) and (xi , yo , zk+2 ).

For example, if the interface crosses the kth mesh line as illustrated in Fig. 1, then
(xi , yj , zk ) and (xi , yj+1 , zk ) are irregular points. The discretizations of ( Φy )y at these two
irregular points contain two fictitious values fi,j,k and fi,j+1,k to extend the domain as

97

follows:

( Φy )y =

( Φy )y =

Φi,j−1,k
i,j− 1
2 ,k

−(

i,j− 12 ,k

+

i,j+ 21 ,k
2

)Φi,j,k +

f
i,j+ 21 ,k i,j+1,k

at

xi , yj , zk

Φ
i,j+ 23 ,k i,j+2,k

at

xi , yj+1 , zk .

(∆y)
f
i,j+ 21 ,k i,j,k

−(

i,j+ 12 ,k

+

i,j+ 32 ,k
2

)Φi,j+1,k +

(∆y)

Herein, the fictitious values are determined by the jump conditions (4) and (5) at the point
of intersection (xi , yo , zk ) between the interface and the y-mesh line.
Since the normal direction from the interface Γ changes according to the position of the
intersecting point, we employ a local coordinate (ξ, η, ζ) at this point, where ξ is the normal
direction, η is the tangential direction and ζ is the binormal direction. This introduction of a
local coordinate system leads the interface problem to become a one-dimensional problem [217].
At the point (xi , yo , zk ), the standard coordinate system is converted to a new coordinate
system using the following matrix equation
 
ξ 
 
 
η  = T
 
 
ζ

  
 
sin φ sin θ cos φ x
x  sin φ cos θ
  
 
  
 
y  =  − sin θ
 ,
cos θ
0 
  
 y 
  
 
z
− cos φ cos θ − cos φ sin θ sin φ
z

where φ and θ are, respectively, the azimuth and zenith angles determined by the normal
direction n. This transformation also results in three jump conditions at each direction:

s
m
s
m
[ Φξ ] = T11 ( s Φsx − m Φm
x ) + T12 ( s Φy − m Φy ) + T13 ( s Φz − m Φz )

(6)

s
m
[Φη ] = T21 (Φsx − Φm
x ) + T22 (Φy − Φy )

(7)

s
m
s
m
[Φζ ] = T31 (Φsx − Φm
x ) + T32 (Φy − Φy ) + T33 (Φz − Φz ),

(8)

98

where Tij is the element in the ith row and jth column of the matrix T . [ Φξ ] is from
the equation (5), but [Φη ] and [Φζ ] are obtained by differentiating the equation (4) [217].
m
m
These jump conditions involve six partial derivatives Φsx , Φsy , Φsz , Φm
x , Φy and Φz . When the

geometry is non-smooth or complex, some of them may be quite challenging to calculate. We
obtain four jump conditions from (4) through (8) that we can use, but the number of the
fictitious values we need to find is only two. Therefore, we do not need to find all of the
above six partial derivatives. In this respect, we want to remove two partial derivatives which
are the most difficult to discretize by two jump conditions [217]. In Fig. 1, Φsx and Φsz are
chosen to be eliminated, which leads to

m
s
m
a1 [ Φξ ] + a2 [Φη ] + a3 [Φζ ] = c1 Φm
x + c2 Φy + c3 Φy + c4 Φz ,

(9)

where a1 = T21 T33 , a2 = s (T31 T13 − T33 T11 ) and a3 = − s T13 T21 . Then the discretizations
of the remaining four partial derivatives require several auxiliary points. Finally, solving a
system of equations (4) and (9) produces the values of fi,j,k and fi,j+1,k . More comprehensive
details are available in the literature [217, 224]. In summary, the MIB scheme manages
interface problems in a well-organized way by using fictitious values, simple Cartesian
coordinate, typical finite difference methods and lower order physical jump conditions [221].

A.3

Iterations of Poisson and Nernst-Planck equations

In our modeling, we employ the following iterative procedure which is also concisely described
in the flowchart (See Fig. 2). First, the boundary value problem of Eq. (2) is required to solve
for Φ only at the first step because this equation does not involve ionic concentrations Cα .
Next, the equations (2.19) and (3) are solved iteratively with suitable boundary conditions of
99

Figure 2: Flowchart of the numerical implementation to solve the standard PNP system

applied voltage and bulk ion concentration until the designated tolerances for Φ and Cα are
reached. To ensure that the numerical process is convergent, Φ and Cα are updated by a
successive over relaxation (SOR)-like iterations

Φnew = 1 − wp Φnew + wp Φold

(10)

Cαnew = (1 − wc ) Cαnew + wc Cαold ,

(11)

where wp and wc are appropriately selected between 0 and 1. This procedure has been already
proved to be efficient for ion channel problems with the second order of accuracy [220].
After the electrostatic potential Φ and the ionic concentration Cα are obtained by solving
Eqs. (2.19), (2) and (3), the electric current is computed at each circular cross section inside

100

the channel along the channel axis [220]. For each ionic species α, its current is calculated by

Iα = qα

Dα
S

∂Cα
qα
∂Φ
+
Cα
dxdy,
∂z
kB T
∂z

(12)

where S is the cross section in the xy-plane. Further, the total channel current is the sum of
two ionic currents.

101

B

Unit system

For real physical simulations with atomic detail, such as biological ion channels and synthetic
nanofluidic channels, the PNP equations cannot be directly treated as a dimensionless
problem [220]. Pertinent units, which are consistent with fundamental physical laws and
central physical constants, should be employed. Consequently, the quantities computed from
the PNP system can be directly compared with experimental observations and measurements.
To this end, we discuss the units of the PNP model in our simulation work.
Table 1: Centimeter-Gram-Second system of units

Abbr.
Unit
cm
centimeter
esu
electrostatic unit
mol
mole
K
Kelvin

Represents
distance
charge
quantity
temperature

In electrostatics, the essential units are based on the centimeter-gram-second system
shown in Table 1 [102]. In our work, we employ the Gaussian units introduced in Table 2
so that we construct the dimensionless electrostatic potential u(r) and the dimensionless
¯ k [102, 220]. Moreover, the required parameters are Avagadro’s number
charges q¯α and Q
NA = 6.0220450 × 1023 , Boltzmann’s constant kB = 1.3806620 × 10−16 erg/K and elementary
charge ec = 4.80320425 × 10−10 esu.
Table 2: Gaussian units

Abbr.
Unit
Represents
Equivalent Expressions
˚
A
angstrom
distance
10−8 cm
l
liter
volume
cm3
M
molar
concentration
mol/l
dyn
dyne
force
esu2 /cm2
erg
erg
energy
dyn·cm

102

In our simulation system, the radii of all atoms and the structural dimensions are given
in terms of angstrom (˚
A). Then the applied voltage is in the unit of volt (V) and the
concentration for each ion species is provided in molar (M).
For the unit conversion as mentioned in our earlier work [220], we define the dimensionless
electrostatic potential u(r) =

ec Φ(r)
¯ k = Qk at kth
and the dimensionless partial charge Q
kB T
ec

atomic component of the channel region. For each ion species α, the dimensionless charge and
the dimensionless concentration are, respectively, defined as q¯α =

qα
e2
and C¯α (r) = c Cα (r).
ec
kB T

Then the electrical charges for K+ and Cl− are, respectively, denoted by q¯1 = q¯K+ = 1 and
q¯2 = q¯Cl− = −1. The Dirichlet boundary conditions for the electrostatic potential and ionic
concentration are modified in the following way:

u0L

Φ0L ec
=
,
300kB T

u0R

Φ0R ec
=
300kB T

2

A e2c
C0N ˚
and C¯α0 = α A
.
1000kB T

Finally, the PNP equations (3) and (2.19) are transformed into

Ns




q¯α C¯α (r),
 − ∇ · ( (r)∇u(r)) = 4π

r∈Ω
(13)

α=1





 ∇ · Dα (r) ∇C¯α (r) + q¯α C¯α (r)∇u(r) = 0,

r ∈ Ωs

together with the Laplace equation for the microscopic channel composition


 ∇2 u(r) = 0,

r ∈ Ωm

(14)


 u(r) = −u∗ (r), r ∈ Γ,

Nm

where

u∗ (r)

=
k=1

¯k
Q
. In our simulations of two subjects, i.e. the MET prototype
m |r − rk |

and ionic diffusive nanofluidic channels, we solve the above three equations with the interface
103

jump conditions [u(r)] and [ (r)u(r)] from Eqs. (4) and (5), and proper boundary conditions.
In addition, we utilize the dimensionless values of the dielectric function (r) which are
empirically determined [102].

(r) =






s

= 80,

r ∈ Ωs
(15)



 m = 2, r ∈ Ωm .
It is also significant to deal with boundary conditions for the electrostatic potential and
concentrations carefully. Theoretically, the electrostatic potential u (r) satisfies the far field
boundary condition, that is,

lim u (r) = 0. In our computation, we need mixed boundary

x→±∞

conditions. Specifically, the Poisson equation uses Dirichlet boundary condition at two
electrodes of ∂Ω where the external voltage is applied [220], i.e., u(r) = u0L at anode and
u(r) = u0R at cathode. Neumann boundary conditions are considered on the other parts of
∂Ω. In practice, the boundary conditions except the applied voltages are somewhat irrelevant
as long as boundaries are sufficiently far from the channel pore [220, 221].
Mixed boundary conditions are also applied to Nernst-Planck equations for all ion species.
Since the Nernst-Planck equation is only defined in the ion inclusion region Ωs , it is enough
to consider the boundary of Ωs . On ∂Ω ∩ ∂Ωs , C¯α = C¯α0 , where C¯α0 represents the bulk ion
concentration for species α. However, at the interface Γ = Ωs ∩ Ωm , the concentration satisfies
the following zero-flux condition −Dα (r) ∇C¯α (r) + q¯α C¯α (r)∇u(r) = 0, where 0 = (0, 0, 0)
is a null vector.

104

C

Atomic charge distribution

C.1

A negatively charged channel

Table 3: Positions and charges of all atomic charges in the negatively charged channel.
k

x(˚
A)

y(˚
A)

z(˚
A)

Qk (ec )

k

x(˚
A)

y(˚
A)

z(˚
A)

Qk (ec )

k

x(˚
A)

y(˚
A)

z(˚
A)

Qk (ec )

k

x(˚
A)

y(˚
A)

z(˚
A)

Qk (ec )

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962

0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962

-23.0000
-23.0000
-23.0000
-23.0000
-23.0000
-23.0000
-23.0000
-23.0000
-21.2963
-21.2963
-21.2963
-21.2963
-21.2963
-21.2963
-21.2963
-21.2963
-19.5926
-19.5926
-19.5926
-19.5926
-19.5926
-19.5926
-19.5926
-19.5926
-17.8889
-17.8889
-17.8889
-17.8889
-17.8889
-17.8889
-17.8889
-17.8889
-16.1852
-16.1852
-16.1852
-16.1852
-16.1852
-16.1852
-16.1852
-16.1852
-14.4815
-14.4815
-14.4815
-14.4815
-14.4815
-14.4815
-14.4815
-14.4815
-12.7778
-12.7778
-12.7778
-12.7778
-12.7778
-12.7778
-12.7778
-12.7778

-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08

57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112

6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962

0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
-4.5962
-6.5000
-4.5962

-11.0741
-11.0741
-11.0741
-11.0741
-11.0741
-11.0741
-11.0741
-11.0741
-9.3704
-9.3704
-9.3704
-9.3704
-9.3704
-9.3704
-9.3704
-9.3704
-7.6667
-7.6667
-7.6667
-7.6667
-7.6667
-7.6667
-7.6667
-7.6667
-5.9630
-5.9630
-5.9630
-5.9630
-5.9630
-5.9630
-5.9630
-5.9630
-4.2593
-4.2593
-4.2593
-4.2593
-4.2593
-4.2593
-4.2593
-4.2593
-2.5556
-2.5556
-2.5556
-2.5556
-2.5556
-2.5556
-2.5556
-2.5556
-0.8519
-0.8519
-0.8519
-0.8519
-0.8519
-0.8519
-0.8519
-0.8519

-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08

113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168

6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962

0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962

0.8519
0.8519
0.8519
0.8519
0.8519
0.8519
0.8519
0.8519
2.5556
2.5556
2.5556
2.5556
2.5556
2.5556
2.5556
2.5556
4.2593
4.2593
4.2593
4.2593
4.2593
4.2593
4.2593
4.2593
5.9630
5.9630
5.9630
5.9630
5.9630
5.9630
5.9630
5.9630
7.6667
7.6667
7.6667
7.6667
7.6667
7.6667
7.6667
7.6667
9.3704
9.3704
9.3704
9.3704
9.3704
9.3704
9.3704
9.3704
11.0741
11.0741
11.0741
11.0741
11.0741
11.0741
11.0741
11.0741

-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08

169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224

6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962

0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962

12.7778
12.7778
12.7778
12.7778
12.7778
12.7778
12.7778
12.7778
14.4815
14.4815
14.4815
14.4815
14.4815
14.4815
14.4815
14.4815
16.1852
16.1852
16.1852
16.1852
16.1852
16.1852
16.1852
16.1852
17.8889
17.8889
17.8889
17.8889
17.8889
17.8889
17.8889
17.8889
19.5926
19.5926
19.5926
19.5926
19.5926
19.5926
19.5926
19.5926
21.2963
21.2963
21.2963
21.2963
21.2963
21.2963
21.2963
21.2963
23.0000
23.0000
23.0000
23.0000
23.0000
23.0000
23.0000
23.0000

-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08

105

C.2

A bipolar channel

Table 4: Positions and charges of all atomic charges in the bipolar channel.
k

x(˚
A)

y(˚
A)

z(˚
A)

Qk (ec )

k

x(˚
A)

y(˚
A)

z(˚
A)

Qk (ec )

k

x(˚
A)

y(˚
A)

z(˚
A)

Qk (ec )

k

x(˚
A)

y(˚
A)

z(˚
A)

Qk (ec )

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962

0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962

-23.0000
-23.0000
-23.0000
-23.0000
-23.0000
-23.0000
-23.0000
-23.0000
-21.2963
-21.2963
-21.2963
-21.2963
-21.2963
-21.2963
-21.2963
-21.2963
-19.5926
-19.5926
-19.5926
-19.5926
-19.5926
-19.5926
-19.5926
-19.5926
-17.8889
-17.8889
-17.8889
-17.8889
-17.8889
-17.8889
-17.8889
-17.8889
-16.1852
-16.1852
-16.1852
-16.1852
-16.1852
-16.1852
-16.1852
-16.1852
-14.4815
-14.4815
-14.4815
-14.4815
-14.4815
-14.4815
-14.4815
-14.4815
-12.7778
-12.7778
-12.7778
-12.7778
-12.7778
-12.7778
-12.7778
-12.7778

0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08

57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112

6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962

0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
-4.5962
-6.5000
-4.5962

-11.0741
-11.0741
-11.0741
-11.0741
-11.0741
-11.0741
-11.0741
-11.0741
-9.3704
-9.3704
-9.3704
-9.3704
-9.3704
-9.3704
-9.3704
-9.3704
-7.6667
-7.6667
-7.6667
-7.6667
-7.6667
-7.6667
-7.6667
-7.6667
-5.9630
-5.9630
-5.9630
-5.9630
-5.9630
-5.9630
-5.9630
-5.9630
-4.2593
-4.2593
-4.2593
-4.2593
-4.2593
-4.2593
-4.2593
-4.2593
-2.5556
-2.5556
-2.5556
-2.5556
-2.5556
-2.5556
-2.5556
-2.5556
-0.8519
-0.8519
-0.8519
-0.8519
-0.8519
-0.8519
-0.8519
-0.8519

0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08
0.08

113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168

6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962

0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962

0.8519
0.8519
0.8519
0.8519
0.8519
0.8519
0.8519
0.8519
2.5556
2.5556
2.5556
2.5556
2.5556
2.5556
2.5556
2.5556
4.2593
4.2593
4.2593
4.2593
4.2593
4.2593
4.2593
4.2593
5.9630
5.9630
5.9630
5.9630
5.9630
5.9630
5.9630
5.9630
7.6667
7.6667
7.6667
7.6667
7.6667
7.6667
7.6667
7.6667
9.3704
9.3704
9.3704
9.3704
9.3704
9.3704
9.3704
9.3704
11.0741
11.0741
11.0741
11.0741
11.0741
11.0741
11.0741
11.0741

-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08

169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224

6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962

0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962

12.7778
12.7778
12.7778
12.7778
12.7778
12.7778
12.7778
12.7778
14.4815
14.4815
14.4815
14.4815
14.4815
14.4815
14.4815
14.4815
16.1852
16.1852
16.1852
16.1852
16.1852
16.1852
16.1852
16.1852
17.8889
17.8889
17.8889
17.8889
17.8889
17.8889
17.8889
17.8889
19.5926
19.5926
19.5926
19.5926
19.5926
19.5926
19.5926
19.5926
21.2963
21.2963
21.2963
21.2963
21.2963
21.2963
21.2963
21.2963
23.0000
23.0000
23.0000
23.0000
23.0000
23.0000
23.0000
23.0000

-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08
-0.08

106

C.3

A double-well channel

Table 5: Positions and charges of all atomic charges in the double well channel.
k

x(˚
A)

y(˚
A)

z(˚
A)

Qk (ec )

k

x(˚
A)

y(˚
A)

z(˚
A)

Qk (ec )

k

x(˚
A)

y(˚
A)

z(˚
A)

Qk (ec )

k

x(˚
A)

y(˚
A)

z(˚
A)

Qk (ec )

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962

0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962

-23.0000
-23.0000
-23.0000
-23.0000
-23.0000
-23.0000
-23.0000
-23.0000
-21.2963
-21.2963
-21.2963
-21.2963
-21.2963
-21.2963
-21.2963
-21.2963
-19.5926
-19.5926
-19.5926
-19.5926
-19.5926
-19.5926
-19.5926
-19.5926
-17.8889
-17.8889
-17.8889
-17.8889
-17.8889
-17.8889
-17.8889
-17.8889
-16.1852
-16.1852
-16.1852
-16.1852
-16.1852
-16.1852
-16.1852
-16.1852
-14.4815
-14.4815
-14.4815
-14.4815
-14.4815
-14.4815
-14.4815
-14.4815
-12.7778
-12.7778
-12.7778
-12.7778
-12.7778
-12.7778
-12.7778
-12.7778

-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12

57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112

6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962

0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
-4.5962
-6.5000
-4.5962

-11.0741
-11.0741
-11.0741
-11.0741
-11.0741
-11.0741
-11.0741
-11.0741
-9.3704
-9.3704
-9.3704
-9.3704
-9.3704
-9.3704
-9.3704
-9.3704
-7.6667
-7.6667
-7.6667
-7.6667
-7.6667
-7.6667
-7.6667
-7.6667
-5.9630
-5.9630
-5.9630
-5.9630
-5.9630
-5.9630
-5.9630
-5.9630
-4.2593
-4.2593
-4.2593
-4.2593
-4.2593
-4.2593
-4.2593
-4.2593
-2.5556
-2.5556
-2.5556
-2.5556
-2.5556
-2.5556
-2.5556
-2.5556
-0.8519
-0.8519
-0.8519
-0.8519
-0.8519
-0.8519
-0.8519
-0.8519

-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04

113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168

6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962

0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962

0.8519
0.8519
0.8519
0.8519
0.8519
0.8519
0.8519
0.8519
2.5556
2.5556
2.5556
2.5556
2.5556
2.5556
2.5556
2.5556
4.2593
4.2593
4.2593
4.2593
4.2593
4.2593
4.2593
4.2593
5.9630
5.9630
5.9630
5.9630
5.9630
5.9630
5.9630
5.9630
7.6667
7.6667
7.6667
7.6667
7.6667
7.6667
7.6667
7.6667
9.3704
9.3704
9.3704
9.3704
9.3704
9.3704
9.3704
9.3704
11.0741
11.0741
11.0741
11.0741
11.0741
11.0741
11.0741
11.0741

0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12

169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224

6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962

0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962
0.0000
4.5962
6.5000
4.5962
0.0000
-4.5962
-6.5000
-4.5962

12.7778
12.7778
12.7778
12.7778
12.7778
12.7778
12.7778
12.7778
14.4815
14.4815
14.4815
14.4815
14.4815
14.4815
14.4815
14.4815
16.1852
16.1852
16.1852
16.1852
16.1852
16.1852
16.1852
16.1852
17.8889
17.8889
17.8889
17.8889
17.8889
17.8889
17.8889
17.8889
19.5926
19.5926
19.5926
19.5926
19.5926
19.5926
19.5926
19.5926
21.2963
21.2963
21.2963
21.2963
21.2963
21.2963
21.2963
21.2963
23.0000
23.0000
23.0000
23.0000
23.0000
23.0000
23.0000
23.0000

-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12
-0.12

107

BIBLIOGRAPHY

108

BIBLOGRAPHY
[1] P. Abgrall and N. T. Nguyen. Nanofluidic devices and their applications. Analytical
chemistry, 80(7):2326–2341, 2008.
[2] H. Adalsteinsson, B. J. Debusschere, K. R. Long, and H. N. Najm. Components for
atomistic-to-continuum multiscale modeling of flow in micro-and nanofluidic systems.
Scientific Programming, 16(4):297–313, 2008.
[3] M. Aguilella-Arzo, J. J. Garc´ıa-Celma, J. Cervera, A. Alcaraz, and V. M. Aguilella.
Electrostatic properties and macroscopic electrodiffusion in OmpF porin and mutants.
Bioelectrochemistry, 70(2):320–327, 2007.
[4] A. Aksimentiev and K. Schulten. Imaging α-hemolysin with molecular dynamics: Ionic
conductance, osmotic permeability, and the electrostatic potential map. Biophysical
journal, 88(6):3745–3761, 2005.
[5] W. L. Ash, M. R. Zlomislic, E. O. Oloo, and D. P. Tieleman. Computer simulations of
membrane proteins. Biochimica et Biophysica Acta (BBA)-Biomembranes, 1666(1):158–
189, 2004.
[6] A. J. Bard and L. R. Faulkner. Electrochemical methods: fundamentals and applications,
volume 2. Wiley New York, 1980.
[7] G. Barreiro, C. R. W. Guimar˜aes, and R. B. de Alencastro. A molecular dynamics
study of an L-type calcium channel model. Protein engineering, 15(2):109–122, 2002.
[8] D. Bashford and D. A. Case. Generalized Born models of macromolecular solvation
effects. Annual Review of Physical Chemistry, 51(1):129–152, 2000.
[9] R. B. Bass, P. Strop, M. Barclay, and D. C. Rees. Crystal structure of Escherichia coli
MscS, a voltage-modulated and mechanosensitive channel. Science, 298(5598):1582–
1587, 2002.
[10] P. Bates, Z. Chen, Y. Sun, G.-W. Wei, and S. Zhao. Geometric and potential driving
formation and evolution of biomolecular surfaces. Journal of mathematical biology,
59(2):193–231, 2009.
[11] P. Bates, G.-W. Wei, and S. Zhao. Minimal molecular surfaces and their applications.
Journal of Computational Chemistry, 29(3):380–391, 2008.
[12] M. Z. Bazant, M. S. Kilic, B. D. Storey, and A. Ajdari. Towards an understanding
of induced-charge electrokinetics at large applied voltages in concentrated solutions.
Advances in colloid and interface science, 152(1):48–88, 2009.

109

[13] M. Z. Bazant, B. D. Storey, and A. A. Kornyshev. Double layer in ionic liquids:
Overscreening versus crowding. Physical Review Letters, 106(4):046102, 2011.
[14] M. Z. Bazant, K. Thornton, and A. Ajdari. Diffuse-charge dynamics in electrochemical
systems. Physical review E, 70(2):021506, 2004.
[15] P. Belgrader, M. Okuzumi, F. Pourahmadi, D. A. Borkholder, and M. A. Northrup. A
microfluidic cartridge to prepare spores for PCR analysis. Biosensors and bioelectronics,
14(10):849–852, 2000.
[16] M. Beurg, R. Fettiplace, J.-H. Nam, and A. J. Ricci. Localization of inner hair cell
mechanotransducer channels using high-speed calcium imaging. Nature neuroscience,
12(5):553–558, 2009.
[17] M. Beurg, J.-H. Nam, Q. Chen, and R. Fettiplace. Calcium balance and mechanotransduction in rat cochlear hair cells. Journal of neurophysiology, 104(1):18–34,
2010.
[18] M. Beurg, J.-H. Nam, A. Crawford, and R. Fettiplace. The actions of calcium on hair
bundle mechanics in mammalian cochlear hair cells. Biophysical journal, 94(7):2639–
2653, 2008.
[19] F. Bezanilla. Ion channels: from conductance to structure. Neuron, 60(3):456–468,
2008.
[20] I. Bir´o, S. Pezeshki, H. Weingart, M. Winterhalter, and U. Kleinekath¨ofer. Comparing
the temperature-dependent conductance of the two structurally similar E. coli Porins
OmpC and OmpF. Biophysical journal, 98(9):1830–1839, 2010.
[21] A. Borisyuk. Physiology and mathematical modeling of the auditory system. In
Tutorials in Mathematical Biosciences I, pages 107–168. Springer, 2005.
[22] D. Branton, D. W. Deamer, A. Marziali, H. Bayley, S. A. Benner, T. Butler, M. Di Ventra, S. Garaj, A. Hibbs, X. Huang, et al. The potential and challenges of nanopore
sequencing. Nature biotechnology, 26(10):1146–1153, 2008.
[23] B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and
M. Karplus. Charmm: A program for macromolecular energy, minimization, and
dynamics calculations. Journal of computational chemistry, 4(2):187–217, 1983.
[24] D. Burch and M. Z. Bazant. Size-dependent spinodal and miscibility gaps for intercalation in nanoparticles. Nano letters, 9(11):3795–3800, 2009.
[25] D. D. Busath, C. D. Thulin, R. W. Hendershot, L. R. Phillips, P. Maughan, C. D.
Cole, N. C. Bingham, S. Morrison, L. C. Baird, R. J. Hendershot, et al. Noncontact
dipole effects on channel permeation. I. Experiments with (5F-Indole) Trp Gramicidin
A channels. Biophysical journal, 75(6):2830–2844, 1998.
[26] H.-J. Butt, K. Graf, and M. Kappl. Physics and chemistry of interfaces. John Wiley &
Sons, 2006.
110

[27] J. Cervera, B. Schiedt, and P. Ramirez. A Poisson/Nernst-Planck model for ionic
transport through synthetic conical nanopores. EPL (Europhysics Letters), 71(1):35,
2005.
[28] G. Chang, R. H. Spencer, A. T. Lee, M. T. Barclay, and D. C. Rees. Structure of
the MscL homolog from Mycobacterium tuberculosis: a gated mechanosensitive ion
channel. Science, 282(5397):2220–2226, 1998.
[29] D. L. Chapman. A contribution to the theory of electrocapillarity. The London,
Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 25(148):475–481,
1913.
[30] D. Chen, Z. Chen, and G.-W. Wei. Quantum dynamics in continuum for proton
transport II: Variational solvent–solute interface. International journal for numerical
methods in biomedical engineering, 28(1):25–51, 2012.
[31] D. Chen and R. Eisenberg. Charges, currents, and potentials in ionic channels of one
conformation. Biophysical journal, 64(5):1405–1421, 1993.
[32] D. Chen and G.-W. Wei. Modeling and simulation of electronic structure, material
interface and random doping in nano-electronic devices. Journal of computational
physics, 229(12):4431–4460, 2010.
[33] D. Chen and G.-W. Wei. Quantum dynamics in continuum for proton transport –
Generalized correlation. The Journal of chemical physics, 136(13):134109, 2012.
[34] D. Chen and G.-W. Wei. Quantum dynamics in continuum for proton transport I:
Basic formulation. Communications in computational physics, 13(1):285, 2013.
[35] D. P. Chen, R. S. Eisenberg, J. W. Jerome, and C.-W. Shu. Hydrodynamic model of
temperature change in open ionic channels. Biophysical journal, 69(6):2304–2322, 1995.
[36] L. Chen and A. Conlisk. Electroosmotic flow and particle transport in micro/nano
nozzles and diffusers. Biomedical microdevices, 10(2):289–298, 2008.
[37] S. Chen and G. D. Doolen. Lattice Boltzmann method for fluid flows. Annual review
of fluid mechanics, 30(1):329–364, 1998.
[38] Z. Chen, N. A. Baker, and G. Wei. Differential geometry based solvation model II:
Lagrangian formulation. Journal of mathematical biology, 63(6):1139–1200, 2011.
[39] Z. Chen, N. A. Baker, and G.-W. Wei. Differential geometry based solvation model I:
Eulerian formulation. Journal of computational physics, 229(22):8231–8258, 2010.
[40] Z. Chen, S. Zhao, J. Chun, D. G. Thomas, N. A. Baker, P. W. Bates, and G. Wei.
Variational approach for nonpolar solvation analysis. The Journal of chemical physics,
137(8):084101, 2012.
[41] L.-J. Cheng and L. J. Guo. Nanofluidic diodes. Chemical Society Reviews, 39(3):923–938,
2010.
111

[42] I.-L. Chern, J.-G. Liu, W.-C. Wang, et al. Accurate evaluation of electrostatics for
macromolecules in solution. Methods and Applications of Analysis, 10(2):309–328, 2003.
[43] Y. S. Choi and S. J. Kim. Electrokinetic flow-induced currents in silica nanofluidic
channels. Journal of colloid and interface science, 333(2):672–678, 2009.
[44] T. Chou. Enhancement of charged macromolecule capture by nanopores in a salt
gradient. The Journal of chemical physics, 131(3):034703, 2009.
[45] K. T. Chu and M. Z. Bazant. Nonlinear electrochemical relaxation around conductors.
Physical Review E, 74(1):011501, 2006.
[46] S.-H. Chung and B. Corry. Three computational methods for studying permeation,
selectivity and dynamics in biological ion channels. Soft Matter, 1(6):417–427, 2005.
[47] S.-H. Chung, M. Hoyles, T. Allen, and S. Kuyucak. Study of ionic currents across a
model membrane channel using Brownian dynamics. Biophysical Journal, 75(2):793–809,
1998.
[48] S.-H. Chung and S. Kuyucak. Recent advances in ion channel research. Biochimica et
Biophysica Acta (BBA)-Biomembranes, 1565(2):267–286, 2002.
[49] R. D. Coalson and M. G. Kurnikova. Poisson–Nernst–Planck theory approach to
the calculation of current through biological ion channels. IEEE transactions on
nanobioscience, 4(1):81–93, 2005.
[50] D. Constantin and Z. S. Siwy. Poisson-Nernst-Planck model of ion current rectification
through a nanofluidic diode. Physical Review E, 76(4):041202, 2007.
[51] W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C.
Spellmeyer, T. Fox, J. W. Caldwell, and P. A. Kollman. A second generation force
field for the simulation of proteins, nucleic acids, and organic molecules. Journal of the
American Chemical Society, 117(19):5179–5197, 1995.
[52] B. Corry, S. Kuyucak, and S.-H. Chung. Tests of continuum theories as models of ion
channels. II. Poisson–Nernst–Planck theory versus Brownian dynamics. Biophysical
Journal, 78(5):2364–2381, 2000.
[53] E. R. Cruz-Chu, A. Aksimentiev, and K. Schulten. Ionic current rectification through
silica nanopores. The Journal of Physical Chemistry C, 113(5):1850–1862, 2009.
´ Cs´anyi, D. Boda, D. Gillespie, and T. Krist´of. Current and selectivity in a model
[54] E.
sodium channel under physiological conditions: Dynamic Monte Carlo simulations.
Biochimica et Biophysica Acta (BBA)-Biomembranes, 1818(3):592–600, 2012.
[55] H. Daiguji. Ion transport in nanofluidic channels. Chemical Society Reviews, 39(3):901–
911, 2010.
[56] H. Daiguji, Y. Oka, and K. Shirono. Nanofluidic diode and bipolar transistor. Nano
Letters, 5(11):2274–2280, 2005.
112

[57] H. Daiguji, P. Yang, and A. Majumdar. Ion transport in nanofluidic channels. Nano
Letters, 4(1):137–142, 2004.
[58] H. Daiguji, P. Yang, A. J. Szeri, and A. Majumdar. Electrochemomechanical energy
conversion in nanofluidic channels. Nano letters, 4(12):2315–2321, 2004.
[59] M. E. Davis and J. A. McCammon. Electrostatics in biomolecular structure and
dynamics. Chemical Reviews, 90(3):509–521, 1990.
[60] B. L. de Groot and H. Grubm¨
uller. Water permeation across biological membranes:
mechanism and dynamics of aquaporin–1 and GlpF. Science, 294(5550):2353–2357,
2001.
[61] P. Debye and E. H¨
uckel. The theory of electrolytes. I. Lowering of freezing point and
related phenomena. Phys. Z, 24:185–206, 1923.
[62] W. Denk, J. R. Holt, G. M. Shepherd, and D. P. Corey. Calcium imaging of single
stereocilia in hair cells: Localization of transduction channels at both ends of tip links.
Neuron, 15(6):1311–1321, 1995.
[63] B. Derjaguin and L. Landau. Theory of the stability of strongly charged lyophobic sols
and of the adhesion of strongly charged particles in solutions of electrolytes. Progress
in Surface Science, 43(1):30–59, 1993.
[64] T. J. Dolinsky, P. Czodrowski, H. Li, J. E. Nielsen, J. H. Jensen, G. Klebe, and N. A.
Baker. PDB2PQR: expanding and upgrading automated preparation of biomolecular
structures for molecular simulations. Nucleic acids research, 35(suppl 2):W522–W525,
2007.
[65] T. J. Dolinsky, J. E. Nielsen, J. A. McCammon, and N. A. Baker. PDB2PQR: an
automated pipeline for the setup of Poisson–Boltzmann electrostatics calculations.
Nucleic acids research, 32(suppl 2):W665–W667, 2004.
[66] B. N. Dominy and C. L. Brooks. Development of a generalized Born model parametrization for proteins and nucleic acids. The Journal of Physical Chemistry B, 103(18):3765–
3773, 1999.
[67] D. A. Doyle, J. M. Cabral, R. A. Pfuetzner, A. Kuo, J. M. Gulbis, S. L. Cohen, B. T.
Chait, and R. MacKinnon. The structure of the potassium channel: molecular basis of
K+ conduction and selectivity. science, 280(5360):69–77, 1998.
[68] C. Duan, W. Wang, and Q. Xie. Review article: Fabrication of nanofluidic devices.
Biomicrofluidics, 7(2):026501, 2013.
[69] R. Dutzler, E. B. Campbell, M. Cadene, B. T. Chait, and R. MacKinnon. X-ray
structure of a ClC chloride channel at 3.0 ˚
A reveals the molecular basis of anion
selectivity. Nature, 415(6869):287–294, 2002.
[70] J. C. Eijkel and A. Van Den Berg. Nanofluidics: what is it and what can we expect
from it? Microfluidics and Nanofluidics, 1(3):249–267, 2005.
113

[71] B. Eisenberg, Y. Hyon, and C. Liu. Energy variational analysis of ions in water and
channels: Field theory for primitive models of complex ionic fluids. The Journal of
chemical physics, 133(10):104104, 2010.
[72] R. Eisenberg and D. Chen. Poisson–Nernst–Planck theory of an open ionic channel.
Biophysical Journal, 64(2):A22–A22, 1993.
[73] D. Erickson, T. Rockwood, T. Emery, A. Scherer, and D. Psaltis. Nanofluidic tuning
of photonic crystal circuits. In Integrated Optoelectronic Devices 2007, pages 647513–
647513. International Society for Optics and Photonics, 2007.
[74] H. Farris, C. LeBlanc, J. Goswami, and A. Ricci. Probing the pore of the auditory hair
cell mechanotransducer channel in turtle. The Journal of physiology, 558(3):769–792,
2004.
[75] X. Feng, K. Xia, Z. Chen, Y. Tong, and G.-W. Wei. Multiscale geometric modeling
of macromolecules II: Lagrangian representation. Journal of computational chemistry,
34(24):2100–2120, 2013.
[76] X. Feng, K. Xia, Y. Tong, and G.-W. Wei. Geometric modeling of subcellular structures,
organelles, and multiprotein complexes. International journal for numerical methods in
biomedical engineering, 28(12):1198–1223, 2012.
[77] R. Fettiplace. Defining features of the hair cell mechanoelectrical transducer channel.
Pfl¨
ugers Archiv-European Journal of Physiology, 458(6):1115–1123, 2009.
[78] R. Fettiplace and C. M. Hackney. The sensory and motor roles of auditory hair cells.
Nature Reviews Neuroscience, 7(1):19–29, 2006.
[79] R. Fettiplace and K. X. Kim. The physiology of mechanoelectrical transduction channels
in hearing. Physiological reviews, 94(3):951–986, 2014.
[80] F. Fogolari and J. M. Briggs. On the variational approach to Poisson–Boltzmann free
energies. Chemical Physics Letters, 281(1):135–139, 1997.
[81] F. Fogolari, A. Brigo, and H. Molinari. The Poisson–Boltzmann equation for biomolecular electrostatics: a tool for structural biology. Journal of Molecular Recognition,
15(6):377–392, 2002.
[82] J. Fu, R. B. Schoch, A. L. Stevens, S. R. Tannenbaum, and J. Han. A patterned
anisotropic nanofluidic sieving structure for continuous-flow separation of DNA and
proteins. Nature nanotechnology, 2(2):121–128, 2007.
[83] M. Gad-el Hak. The fluid mechanics of microdevices - The Freeman scholar lecture.
Journal of Fluids Engineering, 121(1):5–33, 1999.
[84] W. Geng and G.-W. Wei. Multiscale molecular dynamics using the matched interface
and boundary method. Journal of computational physics, 230(2):435–457, 2011.

114

[85] W. Geng, S. Yu, and G. Wei. Treatment of charge singularities in implicit solvent
models. The Journal of chemical physics, 127(11):114106, 2007.
[86] D. Gillespie. Energetics of divalent selectivity in a calcium channel: the ryanodine
receptor case study. Biophysical journal, 94(4):1169–1184, 2008.
[87] D. Gillespie, W. Nonner, and R. S. Eisenberg. Coupling Poisson–Nernst–Planck and
density functional theory to calculate ion flux. Journal of Physics: Condensed Matter,
14(46):12129, 2002.
[88] P. G. Gillespie and U. M¨
uller. Mechanotransduction by hair cells: models, molecules,
and mechanisms. Cell, 139(1):33–44, 2009.
[89] M. K. Gilson, M. E. Davis, B. A. Luty, and J. A. McCammon. Computation of
electrostatic forces on solvated molecules using the Poisson–Boltzmann equation. The
Journal of Physical Chemistry, 97(14):3591–3600, 1993.
[90] G. Gouy. Constitution of the electric charge at the surface of an electrolyte. J. phys,
9(4):457–467, 1910.
[91] D. C. Grahame. The electrical double layer and the theory of electrocapillarity. Chemical
Reviews, 41(3):441–501, 1947.
[92] P. Grochowski and J. Trylska. Continuum molecular electrostatics, salt effects, and
counterion binding – a review of the Poisson–Boltzmann theory and its modifications.
Biopolymers, 89(2):93–113, 2008.
[93] W. Guan, S. X. Li, and M. A. Reed. Voltage gated ion and molecule transport
in engineered nanochannels: theory, fabrication and applications. Nanotechnology,
25(12):122001, 2014.
[94] M. Hacker, W. S. Messer II, and K. A. Bachmann. Pharmacology: principles and
practice. Academic Press, 2009.
[95] J. Han, S. Turner, and H. Craighead. Entropic trapping and escape of long DNA
molecules at submicron size constriction. Physical review letters, 83(8):1688, 1999.
[96] Y. He, D. Gillespie, D. Boda, I. Vlassiouk, R. S. Eisenberg, and Z. S. Siwy. Tuning
transport properties of nanofluidic devices with local charge inversion. Journal of the
American Chemical Society, 131(14):5194–5202, 2009.
[97] J. Hermans, H. J. Berendsen, W. F. Van Gunsteren, and J. P. Postma. A consistent
empirical potential for water–protein interactions. Biopolymers, 23(8):1513–1518, 1984.
[98] B. Hille. Ion channels of excitable membranes, volume 507. Sinauer Sunderland, MA,
2001.
[99] A. L. Hodgkin and A. F. Huxley. The dual effect of membrane potential on sodium
conductance in the giant axon of Loligo. The Journal of physiology, 116(4):497–506,
1952.
115

[100] A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane current
and its application to conduction and excitation in nerve. The Journal of physiology,
117(4):500, 1952.
[101] U. Hollerbach, D.-P. Chen, and R. S. Eisenberg. Two-and three-dimensional Poisson–
Nernst–Planck simulations of current flow through Gramicidin A. Journal of Scientific
Computing, 16(4):373–409, 2001.
[102] M. J. Holst. The Poisson-Boltzmann Equation: Analysis and Multilevel Numerical
Solution. PhD thesis, Numerical Computing Group, University of Illinois at UrbanaChampaign, 1994.
[103] B. Honig and A. Nicholls. Classical electrostatics in biology and chemistry. SCIENCENEW YORK THEN WASHINGTON-, pages 1144–1144, 1995.
[104] J. Howard and A. Hudspeth. Compliance of the hair bundle associated with gating of
mechanoelectrical transduction channels in the bullfrog’s saccular hair cell. Neuron,
1(3):189–199, 1988.
[105] G. Hu and D. Li. Multiscale phenomena in microfluidics and nanofluidics. Chemical
Engineering Science, 62(13):3443–3454, 2007.
[106] Y. Hyon, B. Eisenberg, and C. Liu. A mathematical model for the hard sphere repulsion
in ionic solutions. Commun. Math. Sci, 9(2):459–475, 2011.
[107] W. Im and B. Roux. Ion permeation and selectivity of OmpF porin: a theoretical study
based on molecular dynamics, Brownian dynamics, and continuum electrodiffusion
theory. Journal of molecular biology, 322(4):851–869, 2002.
[108] W. Im, S. Seefeld, and B. Roux. A Grand Canonical Monte Carlo–Brownian dynamics
algorithm for simulating ion channels. Biophysical Journal, 79(2):788–801, 2000.
[109] F. Jaramillo and A. Hudspeth. Localization of the hair cell’s transduction channels
at the hair bundle’s top by iontophoretic application of a channel blocker. Neuron,
7(3):409–420, 1991.
[110] J. W. Jerome. Analysis of charge transport. Springer, 1996.
[111] J. Jiang, D. Cao, D.-e. Jiang, and J. Wu. Time-dependent density functional theory
for ion diffusion in electrochemical systems. Journal of Physics: Condensed Matter,
26(28):284102, 2014.
[112] P. C. Jordan. Fifty years of progress in ion channel research. NanoBioscience, IEEE
Transactions on, 4(1):3–9, 2005.
[113] W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives. Development and testing of
the OPLS all-atom force field on conformational energetics and properties of organic
liquids. Journal of the American Chemical Society, 118(45):11225–11236, 1996.

116

[114] Y.-W. Jung, B. Lu, and M. Mascagni. A computational study of ion conductance in
the KcsA K+ channel using a nernst–planck model with explicit resident ions. The
Journal of chemical physics, 131(21):215101, 2009.
[115] R. Karnik, K. Castelino, R. Fan, P. Yang, and A. Majumdar. Effects of biological
reactions and modifications on conductance of nanofluidic channels. Nano letters,
5(9):1638–1642, 2005.
[116] R. Karnik, K. Castelino, and A. Majumdar. Field-effect control of protein transport in
a nanofluidic transistor circuit. Applied Physics Letters, 88(12):123114, 2006.
[117] R. Karnik, C. Duan, K. Castelino, H. Daiguji, and A. Majumdar. Rectification of ionic
current in a nanofluidic diode. Nano letters, 7(3):547–551, 2007.
[118] R. Karnik, R. Fan, M. Yue, D. Li, P. Yang, and A. Majumdar. Electrostatic control of
ions and molecules in nanofluidic transistors. Nano letters, 5(5):943–948, 2005.
[119] M. Karplus and J. A. McCammon. Molecular dynamics simulations of biomolecules.
Nature Structural & Molecular Biology, 9(9):646–652, 2002.
[120] J. J. Kasianowicz, E. Brandin, D. Branton, and D. W. Deamer. Characterization of
individual polynucleotide molecules using a membrane channel. Proceedings of the
National Academy of Sciences, 93(24):13770–13773, 1996.
[121] P. Kazmierczak and U. M¨
uller. Sensing sound: molecules that orchestrate mechanotransduction by hair cells. Trends in neurosciences, 35(4):220–229, 2012.
[122] P. Kazmierczak, H. Sakaguchi, J. Tokita, E. M. Wilson-Kubalek, R. A. Milligan,
U. M¨
uller, and B. Kachar. Cadherin 23 and protocadherin 15 interact to form tip-link
filaments in sensory hair cells. Nature, 449(7158):87–91, 2007.
[123] H. Kennedy, A. Crawford, and R. Fettiplace. Force generation by mammalian hair
bundles supports a role in cochlear amplification. Nature, 433(7028):880–883, 2005.
[124] H. J. Kennedy, M. G. Evans, A. C. Crawford, and R. Fettiplace. Fast adaptation
of mechanoelectrical transducer channels in mammalian cochlear hair cells. Nature
neuroscience, 6(8):832–836, 2003.
[125] F. Khalili-Araghi, J. Gumbart, P.-C. Wen, M. Sotomayor, E. Tajkhorshid, and K. Schulten. Molecular dynamics simulations of membrane channels and transporters. Current
opinion in structural biology, 19(2):128–137, 2009.
[126] F. Khalili-Araghi, E. Tajkhorshid, and K. Schulten. Dynamics of K+ ion conduction
through Kv1. 2. Biophysical journal, 91(6):L72–L74, 2006.
[127] M. S. Kilic, M. Z. Bazant, and A. Ajdari. Steric effects in the dynamics of electrolytes
at large applied voltages. I. Double-layer charging. Physical Review E, 75(2):021502,
2007.

117

[128] M. S. Kilic, M. Z. Bazant, and A. Ajdari. Steric effects in the dynamics of electrolytes at
large applied voltages. II. Modified Poisson–Nernst–Planck equations. Physical Review
E, 75(2):021503, 2007.
[129] M. G. Kurnikova, R. D. Coalson, P. Graf, and A. Nitzan. A lattice relaxation algorithm
for three–dimensional Poisson–Nernst–Planck theory with application to ion transport
through the Gramicidin A channel. Biophysical Journal, 76(2):642–656, 1999.
[130] G. Lamm. The Poisson–Boltzmann equation. Reviews in Computational Chemistry,
19:147–333, 2003.
[131] K. I. Lee, H. Rui, R. W. Pastor, and W. Im. Brownian dynamics simulations of ion
transport through the VDAC. Biophysical journal, 100(3):611–619, 2011.
[132] Y. Levin. Electrostatic correlations: from plasma to biology. Reports on progress in
physics, 65(11):1577, 2002.
[133] D. G. Levitt. Modeling of ion channels. The Journal of general physiology, 113(6):789–
794, 1999.
[134] J. Li, M. Gershow, D. Stein, E. Brandin, and J. Golovchenko. DNA molecules and
configurations in a solid-state nanopore microscope. Nature materials, 2(9):611–615,
2003.
[135] S. B. Long, E. B. Campbell, and R. MacKinnon. Crystal structure of a mammalian
voltage-dependent shaker family K+ channel. Science, 309(5736):897–903, 2005.
[136] C. D. Lorenz, P. S. Crozier, J. A. Anderson, and A. Travesset. Molecular dynamics
of ionic transport and electrokinetic effects in realistic silica channels. The Journal of
Physical Chemistry C, 112(27):10222–10232, 2008.
[137] C. D. Lorenz and A. Travesset. Charge inversion of divalent ionic solutions in silica
channels. Physical Review E, 75(6):061202, 2007.
[138] B. Lu and Y. Zhou. Poisson–Nernst–Planck equations for simulating biomolecular
diffusion–reaction processes II: size effects on ionic distributions and diffusion–reaction
rates. Biophysical journal, 100(10):2475–2485, 2011.
[139] B. Lu, Y. Zhou, M. Holst, and J. McCammon. Recent progress in numerical methods
for the Poisson–Boltzmann equation in biophysical applications. Commun Comput
Phys, 3(5):973–1009, 2008.
[140] B. Lu, Y. Zhou, G. A. Huber, S. D. Bond, M. J. Holst, and J. A. McCammon.
Electrodiffusion: A continuum modeling framework for biomolecular systems with
realistic spatiotemporal resolution. The Journal of chemical physics, 127(13):135102,
2007.
[141] R. Lynden-Bell and J. C. Rasaiah. Mobility and solvation of ions in channels. The
Journal of chemical physics, 105(20):9266–9280, 1996.
118

[142] D. Mackay, P. H. Berens, K. R. Wilson, and A. Hagler. Structure and dynamics of ion
transport through Gramicidin A. Biophysical journal, 46(2):229–248, 1984.
[143] A. D. MacKerell, D. Bashford, M. Bellott, R. Dunbrack, J. Evanseck, M. J. Field,
S. Fischer, J. Gao, H. Guo, S. a. Ha, et al. All-atom empirical potential for molecular
modeling and dynamics studies of proteins. The Journal of Physical Chemistry B,
102(18):3586–3616, 1998.
[144] R. MacKinnon. Potassium channels and the atomic basis of selective ion conduction
(Nobel Lecture). Angewandte Chemie International Edition, 43(33):4265–4277, 2004.
[145] C. Maffeo and A. Aksimentiev. Structure, dynamics, and ion conductance of the
phospholamban pentamer. Biophysical journal, 96(12):4853–4865, 2009.
[146] C. Maffeo, S. Bhattacharya, J. Yoo, D. Wells, and A. Aksimentiev. Modeling and
simulation of ion channels. Chemical reviews, 112(12):6250–6284, 2012.
[147] M. Manciu and E. Ruckenstein. On the chemical free energy of the electrical double
layer. Langmuir, 19(4):1114–1120, 2003.
[148] D. Marx and J. Hutter. Ab initio molecular dynamics: Theory and implementation.
Modern methods and algorithms of quantum chemistry, 1:301–449, 2000.
[149] J. A. McCammon, B. R. Gelin, and M. Karplus. Dynamics of folded proteins. Nature,
267:585–590, 1977.
[150] Y. Mei, C. Ji, and J. Z. Zhang. A new quantum method for electrostatic solvation
energy of protein. The Journal of chemical physics, 125(9):094906, 2006.
[151] D. Mijatovic, J. Eijkel, and A. Van Den Berg. Technologies for nanofluidic systems:
top-down vs. bottom-up – a review. Lab on a Chip, 5(5):492–500, 2005.
[152] N. Modi, P. R. Singh, K. R. Mahendran, R. Schulz, M. Winterhalter, and
U. Kleinekathofer. Probing the transport of ionic liquids in aqueous solution through
nanopores. The Journal of Physical Chemistry Letters, 2(18):2331–2336, 2011.
[153] N. Modi, M. Winterhalter, and U. Kleinekathoefer. Computational modeling of ion
transport through nanopores. Nanoscale, 4(20):6166–6180, 2012.
[154] M. R. Mohamadi, L. Mahmoudian, N. Kaji, M. Tokeshi, H. Chuman, and Y. Baba.
Nanotechnology for genomics & proteomics. Nano Today, 1(1):38–45, 2006.
[155] M. T. Neves-Petersen and S. B. Petersen. Protein electrostatics: a review of the equations and methods used to model electrostatic equations in biomolecules–applications
in biotechnology. Biotechnology annual review, 9:315–395, 2003.
[156] H. Ohmori. Mechano-electrical transduction currents in isolated vestibular hair cells of
the chick. The Journal of Physiology, 359(1):189–217, 1985.

119

[157] H. Ohshima and K. Furusawa. Electrical Phenomena at Interfaces: Fundamentals,
Measurements, and Applications, volume 76. CRC Press, 1998.
[158] B. Pan, J. Waguespack, M. E. Schnee, C. LeBlanc, and A. J. Ricci. Permeation
properties of the hair cell mechanotransducer channel provide insight into its molecular
structure. Journal of neurophysiology, 107(9):2408–2420, 2012.
[159] A. W. Peng, F. T. Salles, B. Pan, and A. J. Ricci. Integrating the biophysical and molecular mechanisms of auditory hair cell mechanotransduction. Nature communications,
2:523, 2011.
[160] J. L. Perry and S. G. Kandlikar. Review of fabrication of nanochannels for single phase
liquid flow. Microfluidics and Nanofluidics, 2(3):185–193, 2006.
[161] J. M. Perry, K. Zhou, Z. D. Harms, and S. C. Jacobson. Ion transport in nanofluidic
funnels. ACS nano, 4(7):3897–3902, 2010.
[162] S. Prakash, A. Piruska, E. N. Gatimu, P. W. Bohn, J. V. Sweedler, and M. A. Shannon.
Nanofluidics: systems and applications. Sensors Journal, IEEE, 8(5):441–450, 2008.
[163] Q. Pu, J. Yun, H. Temkin, and S. Liu. Ion-enrichment and ion-depletion effect of
nanochannel structures. Nano Letters, 4(6):1099–1103, 2004.
[164] A. Ricci, B. Kachar, J. Gale, and S. Van Netten. Mechano-electrical transduction: new
insights into old ideas. The Journal of membrane biology, 209(2-3):71–88, 2006.
[165] A. J. Ricci, A. C. Crawford, and R. Fettiplace. Tonotopic variation in the conductance
of the hair cell mechanotransducer channel. Neuron, 40(5):983–990, 2003.
[166] F. M. Richards. Areas, volumes, packing and protein structure. Annual review of
biophysics and bioengineering, 6:151–76, 1977.
[167] R. Roth and D. Gillespie.
95(24):247801, 2005.

Physics of size selectivity.

Physical Review Letters,

[168] B. Roux, T. Allen, S. Berneche, and W. Im. Theoretical and computational models of
biological ion channels. Quarterly reviews of biophysics, 37(01):15–103, 2004.
[169] G. Rutkai and T. Krist´of. Dynamic Monte Carlo simulation in mixtures. The Journal
of chemical physics, 132(10):104107, 2010.
[170] M. F. Sanner, A. J. Olson, and J.-C. Spehner. Reduced surface: an efficient way to
compute molecular surfaces. Biopolymers, 38(3):305–320, 1996.
[171] R. B. Schoch, A. Bertsch, and P. Renaud. pH-controlled diffusion of proteins with
different pI values across a nanochannel on a chip. Nano letters, 6(3):543–547, 2006.
[172] R. B. Schoch, J. Han, and P. Renaud. Transport phenomena in nanofluidics. Reviews
of Modern Physics, 80(3):839, 2008.

120

[173] R. B. Schoch and P. Renaud. Ion transport through nanoslits dominated by the effective
surface charge. Applied Physics Letters, 86(25):253111, 2005.
[174] R. Schulz and U. Kleinekath¨ofer. Transitions between closed and open conformations
of TolC: the effects of ions in simulations. Biophysical journal, 96(8):3116–3125, 2009.
[175] M. Schwander, B. Kachar, and U. M¨
uller. The cell biology of hearing. The Journal of
cell biology, 190(1):9–20, 2010.
[176] K. A. Sharp and B. Honig. Calculating total electrostatic energies with the nonlinear
Poisson–Boltzmann equation. Journal of Physical Chemistry, 94(19):7684–7692, 1990.
[177] K. A. Sharp and B. Honig. Electrostatic interactions in macromolecules: theory and
applications. Annual review of biophysics and biophysical chemistry, 19(1):301–332,
1990.
[178] K. Shirono, N. Tatsumi, and H. Daiguji. Molecular simulation of ion transport in silica
nanopores. The Journal of Physical Chemistry B, 113(4):1041–1047, 2009.
[179] S. W. Siu and R. A. B¨ockmann. Electric field effects on membranes: Gramicidin A as
a test ground. Journal of structural biology, 157(3):545–556, 2007.
[180] Z. Siwy, P. Apel, D. Baur, D. D. Dobrev, Y. E. Korchev, R. Neumann, R. Spohr,
C. Trautmann, and K.-O. Voss. Preparation of synthetic nanopores with transport
properties analogous to biological channels. Surface Science, 532:1061–1066, 2003.
[181] Z. Siwy, D. Dobrev, R. Neumann, C. Trautmann, and K. Voss. Electro-responsive
asymmetric nanopores in polyimide with stable ion-current signal. Applied Physics A,
76(5):781–785, 2003.
[182] Z. Siwy, Y. Gu, H. Spohr, D. Baur, A. Wolf-Reber, R. Spohr, P. Apel, and Y. Korchev. Rectification and voltage gating of ion currents in a nanofabricated pore. EPL
(Europhysics Letters), 60(3):349, 2002.
[183] Z. Siwy, E. Heins, C. C. Harrell, P. Kohli, and C. R. Martin. Conical-nanotube ioncurrent rectifiers: the role of surface charge. Journal of the American Chemical Society,
126(35):10850–10851, 2004.
[184] Z. S. Siwy. Ion-current rectification in nanopores and nanotubes with broken symmetry.
Advanced Functional Materials, 16(6):735–746, 2006.
[185] Z. S. Siwy and S. Howorka. Engineered voltage-responsive nanopores. Chemical Society
Reviews, 39(3):1115–1132, 2010.
[186] J. Song, R. Evans, Y.-Y. Lin, B.-N. Hsu, and R. Fair. A scaling model for electrowettingon-dielectric microfluidic actuators. Microfluidics and Nanofluidics, 7(1):75–89, 2009.
[187] M. Sotomayor and K. Schulten. Molecular dynamics study of gating in the mechanosensitive channel of small conductance MscS. Biophysical journal, 87(5):3050–3065, 2004.
121

[188] W. Sparreboom, A. Van Den Berg, and J. Eijkel. Principles and applications of
nanofluidic transport. Nature Nanotechnology, 4(11):713–720, 2009.
[189] W. Sparreboom, A. Van Den Berg, and J. Eijkel. Transport in nanofluidic systems: a
review of theory and applications. New Journal of Physics, 12(1):015004, 2010.
[190] D. Stein, M. Kruithof, and C. Dekker. Surface-charge-governed ion transport in
nanofluidic channels. Physical Review Letters, 93(3):035901, 2004.
[191] S. Sukharev, S. R. Durell, and H. R. Guy. Structural models of the MscL gating
mechanism. Biophysical journal, 81(2):917–936, 2001.
[192] J. Timko, D. Bucher, and S. Kuyucak. Dissociation of NaCl in water from ab initio
molecular dynamics simulations. The Journal of chemical physics, 132(11):114510,
2010.
[193] J. Tomasi, B. Mennucci, and R. Cammi. Quantum mechanical continuum solvation
models. Chemical reviews, 105(8):2999–3094, 2005.
[194] B. Tu, M. Chen, Y. Xie, L. Zhang, B. Eisenberg, and B. Lu. A parallel finite element
simulator for ion transport through three-dimensional ion channel systems. Journal of
computational chemistry, 34(24):2065–2078, 2013.
[195] S. Turner, M. Cabodi, and H. Craighead. Confinement-induced entropic recoil of single
DNA molecules in a nanofluidic structure. Physical Review Letters, 88(12):128103,
2002.
[196] B. M. Venkatesan and R. Bashir. Nanopore sensors for nucleic acid analysis. Nature
nanotechnology, 6(10):615–624, 2011.
[197] E. J. W. Verwey, J. T. G. Overbeek, and J. T. G. Overbeek. Theory of the stability of
lyophobic colloids. Courier Dover Publications, 1999.
[198] V. Vlachy. Ionic effects beyond Poisson–Boltzmann theory. Annual review of physical
chemistry, 50(1):145–165, 1999.
[199] I. Vlassiouk, S. Smirnov, and Z. Siwy. Ionic selectivity of single nanochannels. Nano
letters, 8(7):1978–1985, 2008.
[200] I. Vlassiouk, S. Smirnov, and Z. Siwy. Nanofluidic ionic diodes. Comparison of analytical
and numerical solutions. Acs Nano, 2(8):1589–1602, 2008.
[201] B. Wallace. Recent advances in the high resolution structures of bacterial channels:
Gramicidin A. Journal of structural biology, 121(2):123–141, 1998.
[202] J. Wang, M. Lin, A. Crenshaw, A. Hutchinson, B. Hicks, M. Yeager, S. Berndt, W.Y. Huang, R. B. Hayes, S. J. Chanock, et al. High-throughput single nucleotide
polymorphism genotyping using nanofluidic dynamic arrays. BMC genomics, 10(1):561,
2009.
122

[203] Y. Wang, K. Pant, Z. Chen, G. Wang, W. F. Diffey, P. Ashley, and S. Sundaram.
Numerical analysis of electrokinetic transport in micro-nanofluidic interconnect preconcentrator in hydrodynamic flow. Microfluidics and nanofluidics, 7(5):683–696,
2009.
[204] Y.-C. Wang, A. L. Stevens, and J. Han. Million-fold preconcentration of proteins and
peptides by nanofluidic filter. Analytical chemistry, 77(14):4293–4299, 2005.
[205] J. D. Weeks, D. Chandler, and H. C. Andersen. Role of repulsive forces in determining
the equilibrium structure of simple liquids. The Journal of chemical physics, 54(12):5237–
5247, 1971.
[206] G.-W. Wei. Differential geometry based multiscale models. Bulletin of mathematical
biology, 72(6):1562–1622, 2010.
[207] G.-W. Wei. Multiscale, multiphysics and multidomain models I: Basic theory. Journal
of Theoretical and Computational Chemistry, 12(08), 2013.
[208] G.-W. Wei, Q. Zheng, Z. Chen, and K. Xia. Variational multiscale models for charge
transport. SIAM Review, 54(4):699–754, 2012.
[209] P. K. Weiner and P. A. Kollman. AMBER: Assisted model building with energy
refinement. A general program for modeling molecules and their interactions. Journal
of Computational Chemistry, 2(3):287–303, 1981.
[210] D. Wu and A. J. Steckl. High speed nanofluidic protein accumulator. Lab on a Chip,
9(13):1890–1896, 2009.
[211] J. Wu and Z. Li. Density-functional theory for complex fluids. Annu. Rev. Phys. Chem.,
58:85–112, 2007.
[212] K. Xia, X. Feng, Z. Chen, Y. Tong, and G.-W. Wei. Multiscale geometric modeling of
macromolecules I: Cartesian representation. Journal of computational physics, 257:912–
936, 2014.
[213] K. Xia, K. Opron, and G.-W. Wei. Multiscale multiphysics and multidomain modelsflexibility and rigidity. The Journal of chemical physics, 139(19):194109, 2013.
[214] K. Xia, M. Zhan, and G.-W. Wei. The matched interface and boundary (MIB) method
for multi-domain elliptic interface problems. J. Comput. Phys, 230:8231–8258, 2011.
[215] R. Yan, W. Liang, R. Fan, and P. Yang. Nanofluidic diodes based on nanotube
heterojunctions. Nano letters, 9(11):3820–3825, 2009.
[216] S. Yu, W. Geng, and G. Wei. Treatment of geometric singularities in implicit solvent
models. The Journal of chemical physics, 126(24):244108, 2007.
[217] S. Yu and G. Wei. Three-dimensional matched interface and boundary (MIB) method
for treating geometric singularities. Journal of Computational Physics, 227(1):602–632,
2007.
123

[218] Z. Yuan, A. L. Garcia, G. P. Lopez, and D. N. Petsev. Electrokinetic transport and
separations in fluidic nanochannels. Electrophoresis, 28(4):595–610, 2007.
[219] S. Zhao and G. Wei. High-order FDTD methods via derivative matching for Maxwell’s
equations with material interfaces. Journal of Computational Physics, 200(1):60–103,
2004.
[220] Q. Zheng, D. Chen, and G.-W. Wei. Second–order Poisson–Nernst–Planck solver for
ion transport. Journal of computational physics, 230(13):5239–5262, 2011.
[221] Q. Zheng and G.-W. Wei. Poisson–Boltzmann–Nernst–Planck model. The Journal of
chemical physics, 134(19):194101, 2011.
[222] Z. Zheng, D. J. Hansford, and A. T. Conlisk. Effect of multivalent ions on electroosmotic
flow in micro-and nanochannels. Electrophoresis, 24(17):3006–3017, 2003.
[223] K. Zhou, J. M. Perry, and S. C. Jacobson. Transport and sensing in nanofluidic devices.
Annual Review of Analytical Chemistry, 4:321–341, 2011.
[224] Y. Zhou, M. Feig, and G.-W. Wei. Highly accurate biomolecular electrostatics in
continuum dielectric environments. Journal of Computational Chemistry, 29(1):87–97,
2008.
[225] Y. Zhou, B. Lu, G. A. Huber, M. J. Holst, and J. A. McCammon. Continuum
simulations of acetylcholine consumption by acetylcholinesterase: A Poisson–Nernst–
Planck approach. The Journal of Physical Chemistry B, 112(2):270–275, 2008.
[226] Y. Zhou, S. Zhao, M. Feig, and G.-W. Wei. High order matched interface and boundary
method for elliptic equations with discontinuous coefficients and singular sources.
Journal of Computational Physics, 213(1):1–30, 2006.
[227] Z. Zhou, P. Payne, M. Vasquez, N. Kuhn, and M. Levitt. Finite–difference solution
of the Poisson–Boltzmann equation: Complete elimination of self-energy. Journal of
computational chemistry, 17(11):1344–1351, 1996.

124