# Electroneutrality breakdown in nanopore arrays

J. Pedro de Souza1 Amir Levy2 Martin Z. Bazant1,3 1Department of Chemical
Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts
02139 USA  2Department of Physics, Massachusetts Institute of Technology,
Cambridge, Massachusetts 02139 USA  3Department of Mathematics, Massachusetts
Institute of Technology, Cambridge, Massachusetts 02139 USA

###### Abstract

The electrostatic screening of charge in one-dimensional confinement leads to
long-range breakdown in electroneutrality within a nanopore. Through a series
of continuum simulations, we demonstrate the principles of electroneutrality
breakdown for electrolytes in one dimensional confinement. We show how
interacting pores in a membrane can counteract the phenomenon of
electroneutrality breakdown, eventually returning to electroneutrality.
Emphasis is placed on applying simplifying formulas to reduce the
multidimensional partial differential equations into a single ordinary
differential equation for the electrostatic potential. Dielectric mismatch
between the electrolyte and membrane, pore aspect ratio, and confinement
dimensionality are studied independently, outlining the relevance of
electroneutrality breakdown in nanoporous membranes for selective ion
transport and separations.

##  I Introduction

The conduction of ions in nanochannels forms the basis of electrical
signalling in biology Hille  _et al._ (2001); Eisenberg (1996); Doyle  _et
al._ (1998); Catterall (2000); Sakmann (2013); Maffeo  _et al._ (2012) and of
promising technologies in desalination Shannon  _et al._ (2010); Goh  _et al._
(2016), ionic separations Lu  _et al._ (2020); Li  _et al._ (2020), and
nanofluidic transistors Karnik  _et al._ (2007, 2005); Daiguji (2010). As ions
become confined to nanopores, they experience strong chemical and
electrostatic interactions with the pore walls, leading to membrane
selectivity based on charge or chemical interactions. Nanoporous membranes can
even enter the regime where the double layers emanating from each charged
surface begin to overlap, leading to strong electrokinetic coupling of fluid
flow, electric field, and ionic fluxes Gross and Osterle (1968); Peters  _et
al._ (2016); Biesheuvel and Bazant (2016); Alizadeh and Mani (2017); Qiao and
Aluru (2005a, b). Understanding the electrostatics of charges in confinement
is crucial to determining the flux, selectivity, and driving force
relationship for engineering applications Elimelech and Phillip (2011);
Faucher  _et al._ (2019) and for understanding biological pore systems Branton
_et al._ (2008); Haque  _et al._ (2013); Gillespie  _et al._ (2002); Chen  _et
al._ (1997).

The unique physics of ionic screening in one dimensional confinement leads to
the phenomenon of electroneutrality breakdown, where the number of counter-
charges within a pore does not perfectly counterbalance the number of fixed
charges on the pore walls Levy  _et al._ (2020). One dimensional (1D)
confinement refers to confinement onto a line, for example, in a cylindrical
nanopore connecting two reservoirs of fixed concentration. In essence,
electroneutrality breakdown signifies that a fraction of the electric field
must escape through the pore walls into the dielectric matrix constituting the
membrane. The screening charge does not exist locally within the pore, but
rather is distributed over the membrane surface in the reservoirs, outside of
the membrane domain. Uniquely, in 1D confinement, the loss of
electroneutrality can extend to macroscopic scales (beyond length $L=10$
$\mu$m), since for strong confinement $\kappa_{D}R\rightarrow 0$ the system
only approaches electroneutrality as $\log(L/R)\rightarrow\infty$, where $L$
is the length of a pore, $R$ is its radius, and $\kappa_{D}$ is the inverse
Debye length. The long-range breakdown of local electroneutrality in 1D
confinement is surprising, and many models of electrokinetic phenomena have
assumed pore-wide electroneutrality, even in the limit of strong double layer
overlap Peters  _et al._ (2016); Secchi  _et al._ (2016); Biesheuvel and
Bazant (2016).

Here, we perform continuum simulations of 1D nanopores using COMSOL
Multiphysics, to confirm the occurrence of electroneutrality breakdown. We
show that the results for the screening charge within an isolated 1D pore can
be captured quantitatively with analytical formulas. Furthermore, we perform
simulations of a periodic array of channels, a multipore membrane, with
varying spacing between the channels. There, we find how the interactions of
closely-spaced channels can lead to the return of electroneutrality in the
system. In the limit of strongly interacting channels, the variations in
potential in the axial direction dominate the distribution of ionic charges
within the membrane. Effectively, when channels are too close to each other,
the electric field lines cannot emanate through the membrane domain.

The continuum simulation results show the ensemble interactions of the
channels with each other play a role in ionic conduction through nanochannels.
In membrane applications, the interactions mean that the greatest ionic
selectivity and per-channel-conductivity can be achieved when channels are
close together for the regime of strong double layer overlap. Further, the
results point to ensemble interactions between 1D-confined channels, and the
importance of electric field spilling into the dielectric matrix when channels
are isolated. The competition between channel interactions and
electroneutrality breakdown ultimately affects the conductance and slectivity
behavior of arrays of nanochannels in the low concentration limit.

Figure 1: The simulation configuration in COMSOL (a) with the membrane, pore,
and electrolyte (b) for an isolated pore that does not feel its periodic
neighbors and (c) a periodic arrangement of pores on a square lattice. (d) A
cross-section of the system through the center of the cylindrical pore is
shown to describe the domains and equations applied in the simulations. (e) A
sample of continuum simulation results showing the progress from non-
interacting to interacting channels for small charge densities. For all
channels, the channel radius is 1 nm and the length is 100 nm. The membrane
dielectric constant is $\epsilon_{m}=10\epsilon_{0}$ and the electrolyte
dielectric constant is $\epsilon_{w}=80\epsilon_{0}$. The salt concentration
is 1 mM. In order from left to right, the spacing between the channel centers
$\ell$ is 129 nm, 77 nm, 46 nm, 28 nm, and 17 nm. Also in order, the amount of
charge within the pore versus on the channel walls, $\mid
Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid$ is: 12%,13%, 17%, 26%, and 43%. As the
channels become closer together and more interacting, the system returns
closer to electroneutrality, which is evidenced by the higher $\phi$ values
within the pore.

##  II Theoretical framework

###  II.1 Outline of equations

As explained in more detail in Appendix A, the ionic flux and selectivity out
of equilibrium can be assumed to be related to the solution of the Poisson-
Boltzmann (PB) equation in equilibrium, assuming that both reservoirs on each
side of the membrane have the same concentration of electrolyte and fluid flow
is neglected Gross and Osterle (1968); Peters  _et al._ (2016). As emphasized
in Fig. 1, the simulation is composed of a square membrane domain,
$\Omega_{m}$, with side length $\ell$, through which a cylindrical pore of
radius $R$ and length $L$ connects two reservoirs of fixed concentration. The
reservoirs and pore constitute the electrolyte domain, signified as
$\Omega_{w}$, where $w$ is chosen to signify water.

As exhibited in Fig. 1(a-d), the system of equations being solved in the
electrolyte domain $\Omega_{w}$ is:

$\epsilon_{w}\nabla^{2}\phi=-\rho_{e}$ (1)

with dielectric constant $\epsilon_{w}$, electrostatic potential $\phi$, and
charge density $\rho_{e}$ where the ionic concentrations are Boltzmann-
distributed:

$\rho_{e}=\sum_{i}z_{i}ec_{i,b}\exp\left(-\frac{z_{i}e\phi}{k_{B}T}\right),$
(2)

where $z_{i}$ are the ion valencies, $e$ is the elementary charge, $c_{i,b}$
is the bulk reservoir concentration for ion $i$, $k_{B}$ is the Boltzmann
constant, and $T$ is the absolute temperature. Here we neglect any packing
Bazant  _et al._ (2009), correlation Bazant  _et al._ (2011), or charge
regulation Trefalt  _et al._ (2016) effects in our model. We assume a 1:1
solution of salt with concentration $c_{0}$ in units of number density, such
that the PB equation is reduced to:

$\epsilon_{w}\nabla^{2}\phi=2ec_{0}\sinh\left(\frac{e\phi}{k_{B}T}\right)$ (3)

In the membrane domain, we solve the Laplace equation:

$\epsilon_{m}\nabla^{2}\phi=0$ (4)

assuming that the membrane is a perfect dielectric material with dielectric
constant $\epsilon_{m}$.

At the membrane/electrolyte domain interface, Maxwell’s equation is enforced:

$\mathbf{\hat{n}}\cdot\left[-\epsilon_{m}\nabla\phi+\epsilon_{w}\nabla\phi\right]\Bigg{\rvert}_{s}=q_{s}$
(5)

where $\mathbf{\hat{n}}$ is the unit normal pointing from the electrolyte to
the membrane domain and $q_{s}$ is the surface charge density. To isolate the
electrostatic potential variations due to fixed charge on the channel walls,
$q_{s}$ is assumed to be zero on the membrane/electrolyte reservoir
interfaces, but is nonzero at the pore walls.

At the lateral boundaries of each cell, symmetry conditions are applied;
namely, the electric field at the boundary is zero
$\mathrm{\hat{n}}\cdot\nabla\phi=0$. The boundary conditions are identical to
assuming a periodic array of channels with regular spacing, consisting of unit
cells identical to the simulation box, as illustrated in Fig. 1(c). At the top
and bottom boundaries of the simulation box, Dirichlet conditions are applied
$\phi=0$. A useful check to ensure that the simulation box is large enough is
to ensure that all the integrated charges in the electrolyte domain are equal
and opposite to the integrated amount of fixed charges on the pore walls. For
this study, a reservoir height of 20 $\lambda_{D}$ is sufficient to meet this
criterion, where the Debye length, $\lambda_{D}$ is given by:

$\lambda_{D}={\kappa_{D}}^{-1}=\sqrt{\frac{\epsilon_{w}k_{B}T}{2e^{2}c_{0}}}$
(6)

Figure 2: The results for an isolated pore with overlapping and non-
overlapping double layers within the pore. (a) plot of the potential profile
for overlapping double layers at $c_{0}=1$ mM ($\kappa_{D}R\approx 0.1$) (a)
versus less overlapping double layers at $c_{0}=100$ mM ($\kappa_{D}R\approx
1$)(b). The same parameters are used as in Fig. 1(e), except the center to
center distance between pores is 500 nm such that the periodic channels are
not interacting. For (a) 12% of the charge is contained within the pore,
whereas for (b) 91% of the charge is contained within the pore.(c) The charge
within the pore versus on the pore walls as a function of $\kappa_{D}R$ for
the same channel in parts (a) and (b). Electroneutrality breakdown occurs in
the region of strong double layer overlap $\kappa_{D}R\rightarrow 0$. For (c)
the markers are the COMSOL simulations, whereas the line is the application of
the approximate formula in equation 10. (d-e) The integrated ionic charge as a
function of the lateral position for $c_{0}=$ 1 mM and 100 mM. The charge is
distributed over a wide area $O(L)$ extending beyond the pore mouth when
electroneutrality is broken ($c_{0}$=1 mM), but is more localized when
electroneutrality is maintained ($c_{0}$=100 mM) . (e) Quantification of end
effects for two different concentrations. The electrostatic potential is
plotted as a function of the $z$ coordinate, evaluated at the center axis of
the channel. End effects are not significant for isolated channels.

###  II.2 Approximate formulas

While we compute the full results of the PB equation in 3D, we compare the
results to mathematical simplified formulas. As outlined in more detail in
Appendices B and C, we can reduce our partial differential equation system of
the 3D PB and Laplace equations into ordinary differential equations with
appropriate boundary conditions. For small potentials, we then linearize the
equations and get simple analytical formulas for the number of ionic charges
within the membrane. In order to quantify the extent of electroneutrality
breakdown, we take the ratio for the integrated amount of charge within the
pore and the integrated amount of charge on the pore walls: $\mid
Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid$. In the limit of electroneutrality, we
get $\mid Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid\rightarrow 1$, whereas in the
limit of complete electroneutrality breakdown, we get $\mid
Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid\rightarrow 0$.

Figure 3: Extent of electroneutrality breakdown for channels of different
center-to-center separation distances on a lattice. (a) $\mid
Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid$ versus $\kappa_{D}R$ for $\ell/L=5$,
$\ell/L=1$, $\ell/L=0.1$. (c) $\mid Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid$
versus $\ell/L$ for $c_{0}=0.01,1,100$ mM for $\ell/L=0.02$ to $\ell/L=0.5$.
For both (a) and (c), the markers are the COMSOL simulations, whereas the
solid lines are the application of the approximate formula in equation 10 and
the dotted lines are the application of the approximate formula in equation
14. (b) The integrated ionic charge as a function of the lateral position for
$\ell/L=$ 5 and 0.05, with $c_{0}=$ 1 mM. The charge is distributed over a
wide area when channels are isolated, but is localized when the channels are
closely spaced and strongly interacting. (d) Quantification of end effects for
two different lattice spacings with $c_{0}=$ 1 mM. The electrostatic potential
is plotted as a function of the $z$ coordinate, evaluated at the center axis
of the channel. End effects are significant when the channels are interacting.

No end effects: Ignoring end effects, the inner potential can be solved for in
terms of only the radial coordinate:

$\frac{\epsilon_{w}}{r}\frac{d}{dr}\left(r\frac{d\phi}{dr}\right)=2ec_{0}\sinh\left(\frac{e\phi}{k_{B}T}\right)$
(7)

with boundary conditions given by:

$\frac{d\phi}{dr}(r=R)=\frac{q_{s}}{\epsilon_{w}}-\frac{\epsilon_{m}}{\epsilon_{w}}\frac{\phi(r=R)}{RM_{L/R}},\quad\frac{d\phi}{dr}\left(r=\frac{\ell}{2}\right)=0$
(8)

where the constant $M_{L/R}$ is:

$M_{L/R}=\log\left(\frac{2L}{\pi
R}\right)-\gamma+\frac{\mathrm{K}_{1}(\frac{\pi\ell}{2L})}{\mathrm{I}_{1}(\frac{\pi\ell}{2L})}.$
(9)

with $\gamma=0.577$ signifying the Euler-Mascheroni constant. Here, we assume
circular shape of a unit cell for analytical simplicity, and the derivation is
presented in Appendix B. Linearizing the equations and solving, we find the
folllowing relationship for the amount of charge inside versus the amount of
charge fixed on the pore walls:

$\left\rvert\frac{Q_{\mathrm{in}}}{Q_{\mathrm{out}}}\right\rvert=\frac{1}{1+\frac{\epsilon_{m}}{\epsilon_{w}}\frac{2}{{\kappa_{D}}^{2}R^{2}M_{L/R}}}$
(10)

Figure 4: Role of dielectric mismatch on extent of electroneutrality
breakdown. (a)Results for isolated pore with same properties as in Fig. 1 (e),
but with $\epsilon_{m}=\epsilon_{0}$, $\epsilon_{m}=10\epsilon_{0}$, and
$\epsilon_{m}=100\epsilon_{0}$. (b) $\mid
Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid$ versus $\kappa_{D}R$ for varying
$\epsilon_{m}$. (b) $\mid Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid$ versus
$\ell/L$ with $c_{0}=1$ mM for varying $\epsilon_{m}$.For both (b) and (c),
the markers are the COMSOL simulations, whereas the solid lines are the
application of the approximate formula in equation 10 and the dotted lines are
the application of the approximate formula in equation 14.

Looking at the form of the above equation, one recognizes that
electroneutrality breakdown is promoted as $\kappa_{D}R\rightarrow 0$.
Furthermore, the effect has a weak dependence on the length of non-interacting
nanopores, since the length appears logarithmically in $M_{L/R}$. However, the
electrostatic interactions of the cylinders with each other can cause a return
to electroneutrality as $\ell/L\rightarrow 0$. Furthermore, the amount of
charge within the pore is strongly controlled by the dielectric constant of
the membrane (and not necessarily the inner pore dielectric constant).

No radial variations: On the other hand, when channels become strongly
interacting ($\ell/L\rightarrow 0$), the potential variations happen
principally in the axial direction along the pore, governed by the
dimensionless ratio $\lambda_{D}/L$. In order to capture the end effects of
the channel, which become important on the scale of the Debye length, we write
down a homogenized equation for the potential in linear response, which
neglects radial variations in the potential by integrating over the lateral
dimensions of a unit cell:

$\bar{\epsilon}\frac{d^{2}\phi}{dz^{2}}=-\frac{2q_{s}}{R}+2ec_{0}\sinh\left(\frac{\phi
e}{k_{B}T}\right)$ (11)
$\bar{\epsilon}=\epsilon_{w}-\epsilon_{m}+\epsilon_{m}\left(\frac{\ell}{\pi
R}\right)^{2}$ (12)

The linear response boundary conditions are given by:

$\frac{d\phi}{dz}\left(z=0,L\right)=\pm\frac{\epsilon_{w}}{\bar{\epsilon}}\frac{\ell^{2}}{\pi
R^{2}}\frac{\phi\left(z=0,L\right)}{\lambda_{D}}.$ (13)

Solving the set of equations, we get the following fraction of charge inside
the pore.

$\left\rvert\frac{Q_{\mathrm{in}}}{Q_{\mathrm{out}}}\right\rvert=\frac{1-(\gamma-p)\tanh(1/\gamma)}{1+p\tanh(1/\gamma)}$
(14)

with constants

$\gamma=2\sqrt{\frac{\bar{\epsilon}}{\epsilon_{w}}}\frac{\lambda_{D}}{L},\quad
p=\left(\frac{\bar{\epsilon}}{\epsilon_{w}}\frac{\pi R^{2}}{\ell^{2}}\frac{\pi
R^{2}}{(\ell^{2}-\pi R^{2})}\right)^{1/2}$ (15)

For strongly interacting channels, the electroneutrality condition is
controlled by the ratio of $\lambda_{D}/L$. As $\lambda_{D}/L\rightarrow 0$,
the system returns to overall electroneutrality within the pore.

Following the theoretical argument, a comprehensive set of numerical
simulations is presented to validate the above formulas in their regime of
validity, focusing on the linear regime with small but finite surface charge,
$q_{s}\rightarrow 0$. The standard conditions chosen for the simulations,
unless otherwise stated, are a pore length of 100 nm, a pore radius of 1 nm, a
membrane permittivity of $\epsilon_{m}=10\epsilon_{0}$, an electrolyte
permittivity of $\epsilon_{w}=80\epsilon_{0}$. Finally, a comparison is also
made to the nonlinear solution of the full 3D equations.

The results are presented with a number of dimensionless numbers, so a reader
can easily interpret the plots. First, the dimensionless number $\kappa_{D}R$
indicates the extent of double layer overlap within the channel.
$\kappa_{D}R\rightarrow 0$ indicates strong double layer overlap while
$\kappa_{D}R\rightarrow\infty$ indicates thin double layers relative to the
pore radius. $\ell/L$ is the ratio of the center to center spacing between
channels to the length of the channels. In plotting, the potential is
normalized by the charge per unit length of a cylinder:

$\tilde{\phi}=\frac{\phi\epsilon_{w}R}{2q_{s}\lambda_{D}^{2}}$ (16)

One can roughly interpret these graphs as $\tilde{\phi}\approx 1$ means local
electroneutrality in a give cross section of the pore, and
$\tilde{\phi}\approx 0$ as local electroneutrality within the pore. Further,
the depth-integrated charge density is plotted as a function of lateral
position, to illustrate the extent of screening charges at the
membrane/reservoir interfaces:

$\bar{\rho}_{e}(x,y)=\frac{\int\rho_{e}(x,y,z)dz}{L}.$ (17)

##  III Results and Discussion

First, Fig. 1(e) summarizes the main trends seen in the simulations, where the
potential is plotted as a function of position around the pore, as the spacing
between pores is modified. When channels are isolated ($\ell=129$ nm), the
potential within the channel is fairly constant. The variations in the
potential in the radial direction are dominant. For the parameters chosen
($\kappa_{D}R=0.1$, $c_{0}=$ 1 mM), the amount of charge within the channel is
only 12% of the charge on the pore walls. However, as the size of a unit cell
is reduced ($\ell=17$ nm) the variations in the potential become dominant in
the axial direction, and edge effects become more pronounced. Furthermore, the
potential does not vary significantly in the radial direction from the channel
center axis. The close channel spacing also limits the amount of field that
can escape out of the pore, meaning that more charges are present within the
channel, 43%.

Figure 5: Role of aspect ratio on extent of electroneutrality breakdown.
(a)Results for isolated pore with same properties as in Fig. 1 (e), but with
$L=10$ nm, $L=100$ nm, and $L=1000$ nm, also with $c_{0}=$ 1 mM. (b) $\mid
Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid$ versus $\kappa_{D}R$ for varying $L$,
with $\ell/L=5$. (b) $\mid Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid$ versus
$\ell/L$ with $c_{0}=1$ mM for varying $L$. For both (b) and (c), the markers
are the COMSOL simulations, whereas the solid lines are the application of the
approximate formula in equation 10 and the dotted lines are the application of
the approximate formula in equation 14.

The trends in the potential profiles as a function of the salt concentration
are shown in Fig. 2. Fig. 2(a) corresponds to $\kappa_{D}R=0.1$ or $c_{0}$=1
mM, and 2(b) corresponds to $\kappa_{D}R=1$ of $c_{0}=$ 100 mM, illustrating
that electroneutrality is restored as the concentration is increased. The
solid line in Fig. 2(c) given by equation 10 coincides quite closely with the
results of the continuum simulations for the isolated channels. For these
cases, neglecting end effects is a reasonable approximation for determining
the number of charges within the pore, which is confirmed via the plot of the
potential as a function of position along the center axis of the pore in Fig.
2(f). In Fig. 2(d-e), the depth-integrated charge density, $\bar{\rho}_{e}$,
is plotted as a function of the $x$ and $y$ coordinate. When electroneutrality
is broken within the pore, the screening charge is distributed over the
membrane surface over a distance that is on the order of the channel length,
$L$. However, if electroneutrality is maintained, the screening charge is
localized within the channel, and does not extend very far beyond the pore
mouth.

Figure 6: Role of dimensionality of confinement by inspecting a slit pore
geometry. (a)Results for isolated pore with same properties as in Fig. 1 (e),
but for a slit pore. (b) $\mid Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid$ versus
$\kappa_{D}R$ for varying $\ell/L$. (c) $\mid
Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid$ versus $\ell/L$ for varying $c_{0}$.
(d-e) The same as (b-c) but with $L=1000$ nm. For (b), (c), (d), and (e), the
markers are the COMSOL simulations, whereas the solid lines are the
application of the approximate formula in equation 18 and the dotted lines are
equation 14 with constants given by equations 20 and 21.

Next, we isolate the influence of channels interacting through the membrane
domain, via modifying the ratio $\ell/L$ in Fig. 3. As channels are closer
together, $\ell/L\rightarrow 0$, the amount of charge within the channel is
increased closer towards electroneutrality. We also see a clear difference
between the predictions neglecting end effects (equation 10) in solid lines
and the predictions neglecting radial variations (equation 14) in dotted
lines. As $\ell/L\ll 1$, the predictions of the model neglecting radial
variations become superior. However, when $\ell/L$ is 1 or greater, the
predictions of the model that neglects end effects are superior. Such a result
is expected, given the illustrative example in Fig. 1(e), where the importance
of end effects are visible on the plot when channels are close together. For
further confirmation, the magnitude of the electric field is plotted as a
function of position for varying channel spacings in Fig. S2, exhibiting the
importance of end effects for closely spaced channels. In Fig. 3(d), the
potential plotted as a function of position on the center axis of the nanopore
shows significant end effects for closely spaced channels ($\ell/L$=0.05).
Also, the closely spaced channels localize the screening charge near the pore
mouth, as evidenced by the depth-averaged charge density in Fig. 3(b).

Figure 7: Role of nonlinearity on extent of electroneutrality breakdown.
(a)Results for isolated pore with same properties as in Fig. 1 (e), but with
$q_{s}\rightarrow 0$, $q_{s}=0.001$ C/m2, and $q_{s}=0.01$ C/m2 with $c_{0}=$
1 mM. (b) $\mid Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid$ versus $\kappa_{D}R$ for
varying $q_{s}$. (b) $\mid Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid$ versus
$\ell/L$ with $c_{0}=1$ mM for varying $q_{s}$. For both (b) and (c), the
markers are the COMSOL simulations, whereas the solid line is the application
of the approximate formula in equation 10 and the dotted line is the
application of the approximate formula in equation 14.

The interactions of channels has a profound impact on the selectivity of a
multipore dielectric membrane. The ensemble interactions can lead to a return
to electroneutrality in strongly confined systems, leading to more charge
selectivity and higher channel conductance. The extent of electrokinetic
couplings are also maximized when the amount of ionic charge within the pore
is higher Peters  _et al._ (2016); van der Heyden  _et al._ (2006). Such a
design principle could be used to promote higher electrokinetic conversion
efficiencies in “blue energy” harvesting of salinity gradients Post  _et al._
(2007); He  _et al._ (2017).

Another important design parameter is the dielectric constant of the membrane
matrix. In Figure 4, we examine the influence of the membrane dielectric
constant on the extent of electroneutrality breakdown. As the membrane
dielectric constant decreases, the system moves closer towards
electroneutrality in the pore. Again, the predictions from equations 10 and 14
seem to describe the data quite well within their respective realms of
validity.

A critical question remains to understand the influence of the aspect ratio on
the extent of electroneutrality breakdown, the subject of Figure 5. Again
observing equation 10, for an isolated channel, the electroneutrality is
enforced as $\log(L/R)\rightarrow\infty$. Therefore, the electroneutrality
breakdown is only weakly affected by the length of the channel. As shown in
Figure 5(b), we see only small shifts in the extent of electroneutrality
breakdown with increasing channel length, $L$. Furthermore, the channels
exhibit similar behavior with varying channel spacing, $\ell$. Such results
arise from the exponentially long screening length in one dimensional
confinement Levy  _et al._ (2020).

As a point of comparison, it is instructive to perform the same analysis for
two dimensional confinement, or slit pore geometry. Here, we perform analogous
continuum simulations to the 1D confinement case. In Figure 6, we have
simulated a slit pore with width $2R$ and length $L$, with channel center to
channel center distances of $\ell$. We compare the simulation results to two
analytical formulas, one where we have neglected end effects, and one where we
have neglected radial variations in the potential, similar to the analysis for
1D confinement, as outlined below.

The extent of electroneutrality breakdown in 2D confinement, ignoring edge
effects, is:

$\left\rvert\frac{Q_{\mathrm{in}}}{Q_{\mathrm{out}}}\right\rvert=\frac{1}{1+\frac{\epsilon_{m}}{\epsilon_{w}}\frac{1}{{\kappa_{D}}^{2}R^{2}M_{L/R}}}$
(18)

with modified constant $M_{L/R}$ given by:

$M_{L/R}=\frac{L}{\pi R}\coth\left(\frac{\pi}{L}(\ell/2-R)\right)$ (19)

Observing the form of the 2D $M_{L/R}$ in this scenario, it is apparent that
electroneutrality is enforced as $L/R\rightarrow\infty$. This means that
electroneutrality will be much more strongly upheld in 2D confinement compared
to 1D confinement.

When $\ell/L$ becomes smaller, the potential variations in the axial direction
of the 2D slit pore become dominant, similar to 1D confinement. The extent of
electroneutrality breakdown neglecting radial variations is given identically
by equation 14, but with constants

$\gamma=2\sqrt{\frac{\bar{\epsilon}}{\epsilon_{w}}}\frac{\lambda_{D}}{L},\quad
p=\left(\frac{\bar{\epsilon}}{\epsilon_{w}}\frac{4R^{2}}{\ell^{2}}\right)^{1/2}$
(20)

where

$\bar{\epsilon}=\epsilon_{w}-\epsilon_{m}+\epsilon_{m}\left(\frac{\ell}{2R}\right).$
(21)

The progression towards electroneutrality breakdown for 2D confinement is
similar to 1D confinement when the channel spacing is close together.

In Fig. 6(b-e), we see that the propensity towards electroneutrality is much
stronger in 2D confinement than in 1D confinement. We find that, similar to 1D
confinement, electroneutrality is promoted when channels are strongly
interacting. We also find that the extent of electroneutrality breakdown is
extremely sensitive the to the length of the channel domain, especially as
compared to 1D confinement, as emphasized by comparing Fig. 6(d-e) to Fig.
6(b-c). Note that the normalization of the potential in Fig. 6(a) is adjusted
from the definition in equation 16 due to the difference in geometry.

So far, we have examined the linear regime with small but finite values of
$q_{s}$. Evidently, the nonlinearity in the equations will affect the validity
of our approximations in equations 10 and 14. In Fig. 7, we show that as the
charge density is increased into the nonlinear regime, the system moves closer
towards electroneutrality. The approximations we derived earlier are
insufficient to describe the extent of electroneutrality breakdown in the
nonlinear limit. It is particularly difficult to derive analytical
approximations in the nonlinear limit, so we do not explore such
approximations here. Even so, for isolated channels, the electroneutrality
breakdown emerges as a function of the ratio of the Gouy-Chapman length to the
pore radius in this limit Levy  _et al._ (2020).

Finally, we preview the possible manifestations of electroneutrality breakdown
that could be observed in experiments: single channel conductance and
transference number, using the simplified formulas in equation 10 and equation
14. We assume a KCl solution with fixed and equal mobilites, equal to the bulk
value of $D=D_{+}=D_{-}=2\times 10^{-9}$ m2/s. For an uncharged pore, the
fraction of current carried by each ion would be 50%. However, in this case,
we assume that the pore walls are negatively charged with a value of
$q_{s}=-0.001$ C/m2. Therefore, we approximate the anion and cation
concentrations within the channel as:

$c_{+}=c_{0}+\frac{2\mid q_{s}\mid}{eR}\frac{\mid Q_{\mathrm{in}}\mid}{\mid
Q_{\mathrm{out}}\mid},\quad c_{-}=c_{0}$ (22)

We approximate the overall channel conductance as:

$G=\frac{2\pi DR^{2}e^{2}c_{0}}{k_{B}TL}+\frac{2\pi DRe\mid
q_{s}\mid}{k_{B}TL}\frac{\mid Q_{\mathrm{in}}\mid}{\mid Q_{\mathrm{out}}\mid}$
(23)

which can be rendered dimensionless:

$\tilde{G}=G/\left(\frac{2\pi
DR^{2}e^{2}c_{\mathrm{ref}}}{k_{B}TL}\right)=\frac{c_{0}}{c_{\mathrm{ref}}}+\frac{\mid
q_{s}\mid}{Rec_{\mathrm{ref}}}\frac{\mid Q_{\mathrm{in}}\mid}{\mid
Q_{\mathrm{out}}\mid}$ (24)

where $c_{\mathrm{ref}}$ is arbitrarily chosen to be 1 mM. The corresponding
cation transference number, or fraction of current carried by the cation, is:

$t_{+}=\left(c_{0}+\frac{2\mid q_{s}\mid}{eR}\frac{\mid
Q_{\mathrm{in}}\mid}{\mid
Q_{\mathrm{out}}\mid}\right)/\left(2c_{0}+\frac{2\mid q_{s}\mid}{eR}\frac{\mid
Q_{\mathrm{in}}\mid}{\mid Q_{\mathrm{out}}\mid}\right).$ (25)

The dimensionless conductance and cation transference number are shown in Fig.
8, for a channel with the sample parameters: $L=$ 100 nm, $R=$ 1 nm,
$\epsilon_{m}=10\epsilon_{0}$, $\epsilon_{w}=80\epsilon_{0}$, and $T=$ 300 K.
The plateau in conductance does not occur when electroneutrality breakdown is
present. The decrease in conductance at low concentration has been
experimentally observed in Refs. Secchi  _et al._ (2016) and Yao  _et al._
(2019), but the effects were ascribed to surface reactions and electrokinetic
coupling, respectively. However, as channels become closely spaced or strongly
interacting, their behavior returns to the plateau behavior. Furthermore,
electroneutrality breakdown at low concentrations leads to a cation
transference number that does not saturate at $t_{+}=1$. Another practically
significant quantity is the capacitance of a conducting nanotube embedded in a
membrane dielectric medium, which is explored in the Supplemental Information.
Electroneutrality breakdown can be used to increase the effective capacitance
per unit pore area at low ionic concentrations. However, the capacitance per
total membrane area and per total membrane volume do not benefit from
electroneutrality breakdown, since the dense channel spacing reduces the
effectiveness of electroneutrality breakdown.

Figure 8:  The dimensionless conductance (a) and cation transference number
(b) through a negatively charged nanochannel for varying channel separation
distances. The solid lines are the predictions using equation 10 and the
dotted lines are the predictions using equation 14. The dotted lines are good
predictors of the extent of electroneutrality breakdown at small $\ell/L$, but
fail at large $\ell/L$, e.g. the blue dotted lines. The plateau in conductance
at low concentration is only apparent when $\mid
Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid=1$. As the ratio $\ell/L$ decreases, the
system moves closer to electroneutrality. When electroneutrality is broken,
the cation transference number does not saturate to 1 at low concentration.

Nanoconfined domains are most ubiquitous in biological membranes, where
protein channels selectively conduct specific ions. The selectivity filter in
ion channels are at the molecular scale, with a radii on the order of single
angstroms Malasics  _et al._ (2009). In our model thus far, the radius of the
channel signifies the accessible area for the ion centers, which becomes
negligible in true molecular confinement, $R\rightarrow 0$. Therefore, we
expect electroneutrality breakdown to be critical in describing the conduction
of ions through protein channels. In the absence of an electroneutrality
constraint, molecular separation in the selectivity filter would instead rely
on specific chemical interactions or energy barriers to enter the pore
Epsztein  _et al._ (2019). Closely spaced channels on the order of the
membrane thickness $\sim 10$ nm could interact electrostatically with each
other, leading to ensemble gating and ion conduction events.

The electrokinetic coupling and electrokinetic conversion from salinity
gradients or pressure gradients are most effective and efficient in the regime
of strong double layer overlap. Electroneutrality breakdown can adversely
impact the expected performance of such a process. A large density of channels
is not only desired for higher flux membranes, but to eliminate the
possibility of electroneutrality breakdown at low concentration.

In unstructured charged nanoporous media, such as porous rock or polymer
membranes, a large network of charged pores are connected and interact
strongly with each other. Applications include desalination, ionic
separations, and oil recovery Dydek and Bazant (2013); Deng  _et al._ (2015);
Schlumpberger  _et al._ (2015); Alkhadra  _et al._ (2019); Conforti and Bazant
(2020); Alkhadra  _et al._ (2020). In an interconnected porous medium
consisting of closely spaced pores, a similar homogenized model to equation 14
should be pursued Schmuck and Bazant (2015). If the medium length scale is
large relative to the Debye length, then electroneutrality will be maintained.

With our emphasis on the physics of electroneutrality breakdown, we have
neglected other chemical mechanisms, such as charge regulation Trefalt  _et
al._ (2016); Behrens and Grier (2001); Lin  _et al._ (2019a, b). Charge
regulation, or reactions to form or neutralize surface charge, is certainly
occurring at the pore walls– after all, the origin of surface charge is a
result of charge adsorption to the interface. Even if charge regulation is
present, the effect of electroneutrality breakdown should still be important.
In order to delineate from a chemical mechanism, one convincing evidence of
electric field escape from a channel would be to observe differences in
measureable quantities as a function of the spacing between channels. Only
measurements that determine single channel conductance, transference number,
or capacitance, as a function of the density of channels in a membrane could
distinguish the electroneutrality breakdown mechanism from other competing
chemical mechanisms. Furthermore, in Single Digit Nanopores channels, we can
expect energy barriers associated with ions’ dehydration to enter the channel.
Energy barriers in the pore domain can be easily added to the model explored
here, by changing the effective chemical potential for ions within the pore
Levy  _et al._ (2020). Non-ideal effects Bazant  _et al._ (2009) such as
packing effects Kilic  _et al._ (2007); Kornyshev  _et al._ (2007) or
electrostatic correlations Bazant  _et al._ (2011); Storey and Bazant (2012);
de Souza and Bazant (2020),are also not considered in this work.

##  IV Conclusions

The electroneutrality breakdown phenomenon is studied with an extensive set of
continuum simulations. The results are shown to agree with simplifying
analytical formulas within their regime of validity. Furthermore, the
practical influence of electroneutrality breakdown on channel conductance and
selectivity is discussed.

The experimental validation of screening in lower dimensions presents multiple
competing mechanisms which can obscure the presence of electroneutrality
breakdown. Even so, the set of simulation results presented here can guide
researchers to isolate electroneutrality breakdown for multipore systems. One
can expect wide variations of properties as a function of channel number
density per area when electric fields enter into the membrane domain. As
channels are placed closer together, they interact more strongly, changing the
transport properties.

##  V Acknowledgements

This research was supported as part of the Center for Enhanced Nanofluidic
Transport (CENT), an Energy Frontier Research Center funded by the U.S.
Department of Energy, Office of Science, Basic Energy Sciences under Award #
DE-SC0019112 . JPD acknowledges support from the National Science Foundation
Graduate Research Fellowship under Grant No. 1122374.

## References

  * Hille  _et al._ (2001) B. Hille _et al._ ,  _Ion channels of excitable membranes_ , Vol. 507 (Sinauer Sunderland, MA, 2001).
  * Eisenberg (1996) R. Eisenberg, Journal of Membrane Biology 150, 1 (1996).
  * Doyle  _et al._ (1998) D. A. Doyle, J. M. Cabral, R. A. Pfuetzner, A. Kuo, J. M. Gulbis, S. L. Cohen, B. T. Chait, and R. MacKinnon, science 280, 69 (1998).
  * Catterall (2000) W. A. Catterall, Neuron 26, 13 (2000).
  * Sakmann (2013) B. Sakmann,  _Single-channel recording_ (Springer Science & Business Media, 2013).
  * Maffeo  _et al._ (2012) C. Maffeo, S. Bhattacharya, J. Yoo, D. Wells, and A. Aksimentiev, Chemical Reviews 112, 6250 (2012).
  * Shannon  _et al._ (2010) M. A. Shannon, P. W. Bohn, M. Elimelech, J. G. Georgiadis, B. J. Marinas, and A. M. Mayes, in  _Nanoscience and technology: a collection of reviews from nature Journals_ (World Scientific, 2010) pp. 337–346.
  * Goh  _et al._ (2016) P. Goh, A. Ismail, and N. Hilal, Desalination 380, 100 (2016).
  * Lu  _et al._ (2020) J. Lu, H. Zhang, J. Hou, X. Li, X. Hu, Y. Hu, C. D. Easton, Q. Li, C. Sun, A. W. Thornton, _et al._ , Nature Materials , 1 (2020).
  * Li  _et al._ (2020) Z. Li, Y. Li, Y.-C. Yao, F. Aydin, C. Zhan, Y. Chen, M. Elimelech, T. A. Pham, and A. Noy, ACS nano 14, 6269 (2020).
  * Karnik  _et al._ (2007) R. Karnik, C. Duan, K. Castelino, H. Daiguji, and A. Majumdar, Nano letters 7, 547 (2007).
  * Karnik  _et al._ (2005) R. Karnik, R. Fan, M. Yue, D. Li, P. Yang, and A. Majumdar, Nano letters 5, 943 (2005).
  * Daiguji (2010) H. Daiguji, Chemical Society Reviews 39, 901 (2010).
  * Gross and Osterle (1968) R. J. Gross and J. Osterle, The Journal of chemical physics 49, 228 (1968).
  * Peters  _et al._ (2016) P. Peters, R. Van Roij, M. Z. Bazant, and P. Biesheuvel, Physical review E 93, 053108 (2016).
  * Biesheuvel and Bazant (2016) P. Biesheuvel and M. Bazant, Physical Review E 94, 050601 (2016).
  * Alizadeh and Mani (2017) S. Alizadeh and A. Mani, Langmuir 33, 6205 (2017).
  * Qiao and Aluru (2005a) R. Qiao and N. R. Aluru, Langmuir 21, 8972 (2005a).
  * Qiao and Aluru (2005b) R. Qiao and N. R. Aluru, Applied Physics Letters 86, 143105 (2005b).
  * Elimelech and Phillip (2011) M. Elimelech and W. A. Phillip, science 333, 712 (2011).
  * Faucher  _et al._ (2019) S. Faucher, N. Aluru, M. Z. Bazant, D. Blankschtein, A. H. Brozena, J. Cumings, J. Pedro de Souza, M. Elimelech, R. Epsztein, J. T. Fourkas, _et al._ , The Journal of Physical Chemistry C 123, 21309 (2019).
  * Branton  _et al._ (2008) D. Branton, D. W. Deamer, A. Marziali, H. Bayley, S. A. Benner, T. Butler, M. Di Ventra, S. Garaj, A. Hibbs, X. Huang, S. B. Jovanovich, P. S. Krstic, S. Lindsay, X. S. Ling, C. H. Mastrangelo, A. Meller, J. S. Oliver, Y. V. Pershin, J. M. Ramsey, R. Riehn, G. V. Soni, V. Tabard-Cossa, M. Wanunu, M. Wiggin, and J. A. Schloss, Nature biotechnology 26, 1146 (2008), 18846088[pmid].
  * Haque  _et al._ (2013) F. Haque, J. Li, H.-C. Wu, X.-J. Liang, and P. Guo, Nano Today 8, 56 (2013).
  * Gillespie  _et al._ (2002) D. Gillespie, W. Nonner, and R. S. Eisenberg, Journal of Physics: Condensed Matter 14, 12129 (2002).
  * Chen  _et al._ (1997) D. Chen, J. Lear, and B. Eisenberg, Biophysical Journal 72, 97 (1997).
  * Levy  _et al._ (2020) A. Levy, J. P. de Souza, and M. Z. Bazant, Journal of Colloid and Interface Science (2020).
  * Secchi  _et al._ (2016) E. Secchi, A. Niguès, L. Jubin, A. Siria, and L. Bocquet, Physical review letters 116, 154501 (2016).
  * Bazant  _et al._ (2009) M. Z. Bazant, M. S. Kilic, B. D. Storey, and A. Ajdari, Advances in colloid and interface science 152, 48 (2009).
  * Bazant  _et al._ (2011) M. Z. Bazant, B. D. Storey, and A. A. Kornyshev, Physical Review Letters 106, 046102 (2011).
  * Trefalt  _et al._ (2016) G. Trefalt, S. H. Behrens, and M. Borkovec, Langmuir 32, 380 (2016).
  * van der Heyden  _et al._ (2006) F. H. J. van der Heyden, D. Stein, K. Besteman, S. G. Lemay, and C. Dekker, Physical Review Letters 96, 224502 (2006).
  * Post  _et al._ (2007) J. W. Post, J. Veerman, H. V. Hamelers, G. J. Euverink, S. J. Metz, K. Nymeijer, and C. J. Buisman, Journal of Membrane Science 288, 218 (2007).
  * He  _et al._ (2017) Y. He, Z. Huang, B. Chen, M. Tsutsui, X. Shui Miao, and M. Taniguchi, Scientific Reports 7, 13156 (2017).
  * Yao  _et al._ (2019) Y.-C. Yao, A. Taqieddin, M. A. Alibakhshi, M. Wanunu, N. R. Aluru, and A. Noy, ACS nano 13, 12851 (2019).
  * Malasics  _et al._ (2009) A. Malasics, D. Gillespie, W. Nonner, D. Henderson, B. Eisenberg, and D. Boda, Biochimica et Biophysica Acta (BBA)-Biomembranes 1788, 2471 (2009).
  * Epsztein  _et al._ (2019) R. Epsztein, E. Shaulsky, M. Qin, and M. Elimelech, Journal of Membrane Science 580, 316 (2019).
  * Dydek and Bazant (2013) E. V. Dydek and M. Z. Bazant,  (2013), arXiv:1304.6598 .
  * Deng  _et al._ (2015) D. Deng, W. Aouad, W. A. Braff, S. Schlumpberger, M. E. Suss, and M. Z. Bazant, Desalination 357, 77 (2015).
  * Schlumpberger  _et al._ (2015) S. Schlumpberger, N. B. Lu, M. E. Suss, and M. Z. Bazant, Environmental Science & Technology Letters 2, 367 (2015).
  * Alkhadra  _et al._ (2019) M. A. Alkhadra, K. M. Conforti, T. Gao, H. Tian, and M. Z. Bazant, Environmental Science & Technology 54, 527 (2019).
  * Conforti and Bazant (2020) K. M. Conforti and M. Z. Bazant, AIChE Journal 66, e16751 (2020).
  * Alkhadra  _et al._ (2020) M. A. Alkhadra, T. Gao, K. M. Conforti, H. Tian, and M. Z. Bazant, Desalination 476, 114219 (2020).
  * Schmuck and Bazant (2015) M. Schmuck and M. Z. Bazant, SIAM Journal on Applied Mathematics 75, 1369 (2015).
  * Behrens and Grier (2001) S. H. Behrens and D. G. Grier, The Journal of Chemical Physics 115, 6716 (2001).
  * Lin  _et al._ (2019a) C.-Y. Lin, E. Turker Acar, J. W. Polster, K. Lin, J.-P. Hsu, and Z. S. Siwy, ACS nano 13, 9868 (2019a).
  * Lin  _et al._ (2019b) C.-Y. Lin, P.-H. Wong, P.-H. Wang, Z. S. Siwy, and L.-H. Yeh, ACS Applied Materials & Interfaces (2019b).
  * Kilic  _et al._ (2007) M. S. Kilic, M. Z. Bazant, and A. Ajdari, Physical review E 75, 021502 (2007).
  * Kornyshev  _et al._ (2007) A. A. Kornyshev, D. J. Lee, S. Leikin, and A. Wynveen, Reviews of Modern Physics 79, 943 (2007).
  * Storey and Bazant (2012) B. D. Storey and M. Z. Bazant, Physical Review E 86, 056303 (2012).
  * de Souza and Bazant (2020) J. P. de Souza and M. Z. Bazant, The Journal of Physical Chemistry C 124, 11414 (2020).

Figure S1: The configuration of multiple channels in a membrane (a) is reduced
to a single pore problem (b). The partial differential equation and boundary
conditions are listed for the given single pore problem.

##  Appendix A Reduction of the PNP system into the PB system

The PNP system of equations in the electrolyte domain that applies for
transport measurements is given by coupling the Poisson equation to the
Nernst-Planck equation for dilute species:

$\begin{split}&\epsilon_{w}\nabla^{2}\phi=-\rho_{e}\\\
&\mathbf{j_{i}}=-\frac{D_{i}c_{i}}{k_{B}T}\nabla\mu_{i}^{e}\\\
&\mu_{i}^{e}=k_{B}T\ln(c_{i})+z_{i}e\phi\end{split}$ (26)

with steady state species conservation equations given by:

$\nabla\cdot\mathbf{j_{i}}=0.$ (27)

Since the channels we are considering are long and thin, we can neglect the
change in flux in the normal direction to the pore walls Peters  _et al._
(2016), which is zero due to no penetration at the pore walls:

$j_{i,n}=-D_{i}\nabla_{n}c_{i}-\frac{D_{i}z_{i}ec_{i}}{k_{B}T}\nabla_{n}\phi=0$
(28)

Integrating with respect to the normal coordinate, we can define a Boltzmann
distribution with:

$c_{i}(r,z)=c_{i,v}(z)\exp\left(\frac{-z_{i}e\phi}{k_{B}T}\right).$ (29)

Since both the open ends of the channel are fixed at approximately the
reservoir concentration, we can choose the pre-exponential factor as
$c_{i},0$. At this point, we have reduced the dynamic Nernst-Planck equation
into the equilibrium Poisson-Boltzmann equation. Driving forces applied to the
system will linearly act upon the equilibrium configuration, as calculated in
this paper. For driving forces that are not too large, we can expect the
equilibrium conformation to also remain unperturbed by driving forces. Here,
we have neglected the coupling to the Stokes equation, but it can be easily
added to the analysis Peters  _et al._ (2016).

Figure S2:  Electric field magnitude normalized by the surface charge as a
function of position within a cross section of the membrane/channel system,
for the same conditions as Fig. 1(e). As the unit cell becomes smaller and
channels are closer together, the electric field lines are strongest at the
entrance regions of the channel. End effects can be important in this regime.

##  Appendix B Neglecting end effects

First, we derive the approximate formulas for the number of ions within the
channel by neglecting end effects, based on the work Levy  _et al._ (2020) but
extended to a nanopore array configuration. In order to simplify the formulas,
we consider a cylindrical unit cell that encapsulates a single nanochannel, as
shown in Fig. S1. In the membrane domain, because the Laplace equation is long
range, we must account for the variations in the axial direction:

$\frac{1}{r}\frac{d}{dr}\left(r\frac{d\phi_{\mathrm{out}}}{dr}\right)+\frac{d^{2}\phi_{\mathrm{out}}}{dz^{2}}=0$
(30)

With homogeneous boundary conditions $\phi_{1}=\phi_{2}=0$, the solutions have
the form:

$\begin{split}&\phi_{\mathrm{out}}(z,r)=\sum_{n=1}^{\infty}A_{n}\left[\mathrm{K}_{0}(\lambda_{n}r)+\alpha_{n}I_{0}(\lambda_{n}r)\right]\sin(\lambda_{n}z)\\\
&\alpha_{n}=\frac{\mathrm{K}_{1}(\frac{\lambda_{n}\ell}{2})}{\mathrm{I}_{1}(\frac{\lambda_{n}\ell}{2})},\quad\quad\lambda_{n}=\frac{n\pi}{L}\end{split}$
(31)

The normal derivative to the channel wall is:

$\frac{\partial\phi_{\mathrm{out}}}{\partial
r}(z,r)=\sum_{n=1}^{\infty}A_{n}\left[-\lambda_{n}\mathrm{K}_{1}(\lambda_{n}r)+\alpha_{n}\lambda_{n}I_{1}(\lambda_{n}r)\right]\sin(\lambda_{n}z).$
(32)

Far from the electrolyte/membrane interface, the first eigenvalue dominates
the solution. We can approximate each summation asymptotically by truncating
to the first term:

$\begin{split}&\phi_{\mathrm{out}}(z,R)\sim
A_{1}\left[\mathrm{K}_{0}(\lambda_{1}R)+\alpha_{1}I_{0}(\lambda_{1}r)\right]\sin(\lambda_{1}z)\\\
&\frac{\partial\phi_{\mathrm{out}}}{\partial r}(z,R)\sim
A_{1}\left[-\lambda_{1}\mathrm{K}_{1}(\lambda_{1}R)+\alpha_{1}\lambda_{1}I_{1}(\lambda_{1}R)\right]\sin(\lambda_{1}z)\end{split}$
(33)

Taking the limit as $R/L\rightarrow 0$, we get:

$\begin{split}&\phi_{\mathrm{out}}(z,R)\sim A_{1}\left[-\gamma-\log(\frac{\pi
R}{2L})+\alpha_{1}\right]\sin(\lambda_{1}z)\\\
&\frac{\partial\phi_{\mathrm{out}}}{\partial
r}(z,R)\sim-\frac{A_{1}}{R}\sin(\lambda_{1}z)\end{split}$ (34)

where $\gamma$ is the Euler constant. Taking the ratio of the potential and
its derivative at the external radius, we get:

$\frac{\partial\phi_{\mathrm{out}}}{\partial
r}(z,R)=\frac{-\phi_{\mathrm{out}}(z,R)}{RM_{L/R}}$ (35)

with $M_{L/R}$ defined by equation 9. If we apply the above relationship into
the electric flux boundary conditions, we get the boundary condition in
equation 8. Next, we reduce the Poisson-Boltzmann equation from a partial
differential equation to an ordinary differential equation by neglecting the
variation in potential in the electrolyte domain in the $z$ direction:

$\frac{\epsilon_{w}}{r}\frac{d}{dr}\left(r\frac{d\phi}{dr}\right)=2ec_{0}\sinh\left(\frac{e\phi}{k_{B}T}\right).$
(36)

Solving the approximate boundary conditions in the linear response, we can
arrive at the equation for the amount of charge in the pore relative to the
charge on the pore walls, given in equation 10.

##  Appendix C Neglecting radial variations

As shown in Fig. S2, when channels become closely spaced, the electric field
is concentrated at the ends of the channel, meaning that end effects become
more important. In order to turn our partial differential equation system into
an ordinary differential equation system valid in this regime, we average over
the $x$ and $y$ dimensions of the unit cell. Integrating the Poisson and
Laplace equation system over $x$ and $y$, dividing by $\pi R^{2}$, and
assuming no variations in $\phi$ in the $x$ or $y$ directions, we get
equations 11 and 12. Next, at linear response, we can relate the derivative of
the potential at the interface to the value of the potential itself at the
interface, since the potential should exponentially decay into the reservoir
domain. Mathematically, the relationship is given in equation 13. Note that
the approximations we have made here are only valid if the potential is
constant in each $z$-slice of our simulation. Such an approximation is
expected to fail at large channel separation distances, as shown throughout
our results.

##  Appendix D Unit cell shape

Here, we explore whether the shape of the unit cell affects the extent of
electroneutrality breakdown. We find that the results for a hexagonal unit
cell with the same center to center distance as a square unit cell exhibits
the same extent of electroneutrality breakdown, as shown in Fig. S3. While
certain differences may arise by fluctuations from the lattice position or
irregular lattice spacing of channels, we have provided some additional
evidence that our square array is a good approximation for a variety of
channel arrangements.

Figure S3:  The extent of electroneutrality breakdown for two different unit
cell shapes: square and hexagonal. The center to center distance between
channels, $\ell$ is maintained the same between the channels. The conditions
are identical to Fig. 3(c). As $\ell$ is varied, the hexagonal and square unit
cells appear to have similar extents of electroneutrality breakdown.  Figure
S4: The capacitance (a) normalized by the pore surface area (b) normalized by
the membrane/electrolyte surface area per unit cell and (d) normalized by the
membrane/pore volume. In each plot, the separation distance between channels
is varied, changing the extent of electroneutrality breakdown. The solid lines
are the predictions using equation 10 and the dotted lines are the predictions
using equation 14. The dotted lines are good predictors of the extent of
electroneutrality breakdown at small $\ell/L$, but fail at large $\ell/L$,
e.g. the blue dotted lines.

##  Appendix E Capacitance

As a final practical application of electroneutrality breakdown, we examine
the capacitance of a nanopore that experiences electroneutrality breakdown. As
a device, we imagine an electrically-connected nanopore embedded in a
dielectric membrane. Scaled up, a set of nanopores is arranged in an array
with center to center spacing of $\ell$, separated by a membrane dielectric
domain. Using the same parameters as in Fig. 8, we calculate the differential
capacitance within the channel at linear response:

$C=\frac{dq_{s}}{d\phi}=\frac{\epsilon_{w}R}{2\lambda_{D}^{2}}/\mid
Q_{\mathrm{in}}/Q_{\mathrm{out}}\mid$ (37)

valid for strongly overlapping double layers $\kappa_{D}R\rightarrow 0$, and
normalized by the surface area of the charged nanopore, $2\pi RL$. In Fig.
S4(a), the capcacitance is plotted, showing that for channels that have
electroneutrality breakdown $\ell/L>1$ at low concentrations, that the
capacitance saturates at low concentration. The value of the capacitance at
low concentration is higher for channels which exhibit electroneutrality
breakdown. If we renormalize the capacitance by the total surface area of the
membrane/electrolyte domain per unit cell, we arrive at the results in Fig.
S4(b). If we renormalize the capacitance by the volume of the membrane/pore
system, we get the results in Fig. S4(c). These renormalized capacitance
values exhibit the opposite trend–namely, as channels are closer together
($\ell/L\rightarrow 0$), the cumulative surface area and volume decrease,
leading to greater renormalized capacitance. Even though the single channels
have more capacitance per pore surface area due to electroneutrality
breakdown, the additional surface area or volume for the channels to be spaced
far apart from each other leads to smaller renormalized capacitance. One could
potentially design a membrane with significant density which still exhibits
some electroneutrality breakdown, leading to better performance than a denser
array of channels, as shown in the comparison between the green and red lines
in Fig. S4.


# Learning with Chemical versus Electrical Synapses
Does it Make a Difference?

Mónika Farsang1∗, Mathias Lechner2, David Lung1, Ramin Hasani2, Daniela Rus2,
Radu Grosu1 1CPS, Technische Universität Wien (TU Wien), Austria2CSAIL,
Massachusetts Institute of Technology (MIT), MA,
<EMAIL_ADDRESS>project has received funding from the
European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie grant agreement No 101034277.

###### Abstract

Bio-inspired neural networks have the potential to advance our understanding
of neural computation and improve the state-of-the-art of AI systems. Bio-
electrical synapses directly transmit neural signals, by enabling fast current
flow between neurons. In contrast, bio-chemical synapses transmit neural
signals indirectly, through neurotransmitters. Prior work showed that
interpretable dynamics for complex robotic control, can be achieved by using
chemical synapses, within a sparse, bio-inspired architecture, called Neural
Circuit Policies (NCPs). However, a comparison of these two synaptic models,
within the same architecture, remains an unexplored area. In this work we aim
to determine the impact of using chemical synapses compared to electrical
synapses, in both sparse and all-to-all connected networks. We conduct
experiments with autonomous lane-keeping through a photorealistic autonomous
driving simulator to evaluate their performance under diverse conditions and
in the presence of noise. The experiments highlight the substantial influence
of the architectural and synaptic-model choices, respectively. Our results
show that employing chemical synapses yields noticeable improvements compared
to electrical synapses, and that NCPs lead to better results in both synaptic
models.

##  I INTRODUCTION

Continuous-time recurrent neural networks (CT-RNNs) are well-suited for
robotic and control applications, because they are able to learn dependencies
and patterns occurring in CT processes [1, 2, 3, 4]. These networks are
described by ordinary differential equations (ODEs) representing the time
evolution of the activation of neurons over time, based on the inputs and the
internal dynamics of the network.

Most robotic and control applications are safety-critical. Their faulty
operation may cause human harm or physical damage. Interpretable machine
learning methods are thus a natural fit for robotic applications, as they
allow the system behavior to be studied and audited by humans.

Liquid-time-constant neural networks (LTCs), a class of biology-inspired
neural models with chemical synapses [4], have recently demonstrated
interpretable-control capabilities in autonomous-driving and drone-navigation
tasks [5, 6, 7], when organized in a sparse 4-layer Neural Circuit Policy
architecture (NCP) [8]. However, the individual contributions of the NCP
architecture and the synaptic model to interpretable control, is still poorly
understood.

In this paper, we first show that electrical biologic synapses underlay CT-
RNNs [2], and that chemical synapses underlay CT-RNNs with synapse-bound
activations, and input-dependent forget gates [9]. We then study how much
learning with electrical versus chemical synapses really makes a difference in
the network’s performance when they share the same type of architecture. For a
systematic comparison, we experiment with LTCs [4] and classic CT-RNNs [2], in
both NCP [5] and fully-connected architectures.

We evaluate the performance of these networks in the complex task of
autonomous car-navigation across different seasons. We first train the
networks using imitation learning [10], where they learn to imitate the
actions of a human driver through expert demonstrations. Subsequently, we
deploy the trained networks in closed-loop settings to examine their
adaptivity to previously unseen scenarios and their ability to react to
environmental feedback. These evaluations are conducted using the data-driven
photorealistic autonomous driving simulator VISTA [11]. This has been proven
to yield consistent results, when comparing trained agents in closed-loop
evaluations in simulation, and real-world tests.

Our main contributions are the following:

  * •

We show that saturated electrical synapses are the foundation of classic
continuous-time RNNs, and that saturated chemical synapses underlay CT-RNNs
with synapse-bound activations and input-dependent forget gates.

  * •

We conduct a comparative analysis between electrical and chemical synaptic
models, respectively, in sparse NCP and fully-connected wiring architectures,
evaluating their performance, interpretability and robustness.

  * •

Our results show that: (1) The architectural wiring choice significantly
impacts the network performance, with NCPs exhibiting better results with both
synaptic models, and (2) Employing chemical synapses yields noticeable
improvements compared to electrical synapses.

This paper is structured as follows: In Section II we discuss the bio-inspired
models of electrical and chemical synapses, and the wiring model of NCPs.
Section III provides a summary of related work on imitation learning, lane-
keeping, and applications of the models under investigation. In Section IV, we
outline the metrics of interest and detail the experimental setup. The results
are presented in Section V, followed by a concluding discussion in Section VI.

##  II SYNAPTIC MODELS

The dynamic behavior of neurons is captured in neuroscience with electrical
equivalent circuits (EECs) [12], [13]. In EECs, a neuron-$i$’s membrane is
regarded as a capacitor whose potential $x_{i}$ is determined by the
difference between intra- and extra-cellular concentrations of ions,
respectively. The membrane potential thus evolves as follows:

$C\,\dot{x}_{i}=i_{li}+i_{si},\quad i_{li}=g_{li}(e_{li}-x_{i})$ (1)

where $C$ is the membrane’s capacitance, and $i_{li}$ and $i_{si}$ are the
passive-leaking and dynamic-synaptic currents, traversing the capacitor,
respectively. Without loss of generality, we will assume that $C\,{=}\,1$ in
the rest of the paper. The leaking current, is determined by the Ohm’s
equation, where $e_{li}$ is the resting (or reversal) potential of the neuron,
and $g_{li}$ is its leaking conductance. They are considered to be constant.

####  II-1 Electrical synapses

occur in all organisms, especially in embryonic stage. Their prevalent form is
that of an Ohmic current $i_{s,ji}\,{=}g_{ji}(y_{j}\,{-}\,x_{i})$, where
$y_{j}$ is either the potential $x_{j}$ of a presynaptic neuron, or the
potential $u_{j}$ of an input. Constant $g_{ji}$ is the conductance of the
synapse. Hence, in electrical synapses, $x_{i}$ faithfully follows $y_{j}$.
Let $y\,{=}[x,u]$, $|x|\,{=}\,m$ and $|u|\,{=}\,n$. As a network, we thus
obtain:

$\dot{x}_{i}=g_{li}(e_{li}-x_{i})+\sum_{j=1}^{m+n}g_{ji}(y_{j}-x_{i})$ (2)

Rearranging this equation in terms of $x_{i}$ and $e_{li}$, one obtains:

$\dot{x}_{i}=-(g_{li}+\sum_{j=1}^{m+n}g_{ji})\,x_{i}+(g_{li}+\sum_{j=1}^{m+n}\frac{g_{ji}}{e_{li}}y_{j})e_{li}$
(3)

Equations (2)-(3) are the EECs of neural networks with electrical synapses,
only. However, neurons saturate their overall conductance, through diverse
cellular mechanisms. A more faithful and numerically stable representation is
thus:

$\dot{x}_{i}=-S(g_{li}+\sum_{j=1}^{m}g_{ji})\,x_{i}+T(g_{li}+\sum_{j=1}^{m}h_{ji}y_{j})e_{li}$
(4)

where $S$ and $T$ are sigmoidal and hyperbolic-tangent functions,
respectively, and $h_{ji}\,{=}\,g_{ji}\,{/}\,e_{li}$. This equation is
normalized, with conductances $g_{li},g_{ji}\,{\in}\,[0,1]$. Hence, saturation
is with $S$. Since $h_{ij}\,{\in}\,[-1,1]$ for $e_{li}\,{\in}\,[-1,1]$, we use
$T$.

Taking $\tau\,{=}\,S(g_{li}+\sum_{j=1}^{m}g_{ij})$ one obtains the classic CT-
RNNs [2], with constant gate $\tau$. Thus, classic CT-RNNs are neural networks
with electrical synapses, only!

####  II-2 Chemical synapses

occur in all organisms, especially in later stages. A current
$i_{s,ji}\,{=}g_{ji}S(a_{ji}y_{j}+b_{ji})(e_{ji}\,{-}\,x_{i})$ passes through
a postsynaptic neuron’s channels, when these open, as neurotransmitters
released by the presynaptic neuron start binding to the receptors of the
channels. Here, $g_{ji}$ is the maximum conductance of the channels,
$S(a_{ji}y_{j}+b_{ji})$ is the probability of the channels to be open, with
parameters $a_{ji}$ and $b_{ji}$, and $e_{ji}$ is the reversal potential of
the channels.

Note that in a chemical synapse, $x_{i}$ does not necessarily follow $y_{j}$.
The synaptic current $i_{s,ji}$ is 0 even if all channels are open, when
$e_{ji}\,{=}\,x_{i}$. A network is an ODEs system:

$\dot{x}_{i}=g_{li}(e_{li}-x_{i})+\sum_{j=1}^{m+n}g_{ji}S(a_{ji}y_{j}+b_{ji})(e_{ji}-x_{i})$
(5)

These are the EECs (or LTCs [4]) of chemical synapses. They also ignore
saturation aspects. Taking those into account, and rearranging terms as
before, one obtains SEECs (or SLTCs):

$\begin{array}[]{c}\dot{x}_{i}=-S(f_{i})\,x_{i}+T(g_{i})\,e_{li}\\\\[5.69054pt]
f_{i}=g_{li}+\sum_{j=1}^{m+n}g_{ji}S(a_{ji}y_{j}+b_{ji})\\\\[5.69054pt]
g_{i}=g_{li}+\sum_{j=1}^{m+n}h_{ji}S(a_{ji}y_{j}+b_{ji})\end{array}$ (6)

As in gated RNNs [9], $S(f_{i})$ is the forget gate. However, the other gate
is just $e_{li}$, the leaking potential. Moreover, each synapse has a
different activation function (different $a_{ji}$, $b_{ji}$).

####  II-3 Neural Circuit Policies

for short NCPs, consist of 4 layers [8]. The incoming information is received
by the sensory neurons. This information is forwarded to interneurons, playing
the role of a self-attention layer. Command neurons, which are connected in a
recurrent fashion, store the state of the NCP, and use this together with the
attention information, to compute actions. These are forwarded to the motor
neurons. The NCP design aims to create a sparse architecture with
interpretable dynamics. In contrast, commonly used fully-connected RNN
architectures, do not distinguish layers of neurons, and neurons are connected
all-to-all.

##  III RELATED WORK

Imitation learning is about training the agents to learn the connections
between perceptions and actions, by observing expert demonstrations of a task
[10]. Its application extends across various domains, including robotic
manipulation [14], locomotion [15], and autonomous navigation, both on the
ground [5, 16] and in aerial environments [17, 6, 7].

Autonomous lane-keeping solutions include end-to-end convolutional neural
networks (ConvNet) [18, 19], and more advanced solutions of ConvNet heads
stacked with RNNs to consider previous information during decision making [5].
Here, we use the latter variant, as similarly, we need a ConvNet head to
preprocess the raw visual input which will be fed into recurrent chemical or
electrical models.

The CT-RNN model of electrical synapses is a commonly used neural model, that
has been applied in walking of a six-legged agent [20], robot-to-human
handover [21], balancing a biped robot [22] and swarm robotics [23]. Recently,
the LTC model of chemical synapses has demonstrated the ability to
autonomously navigate both cars [5] and drones [6, 7].

An NCP wiring architecture built up from 19 LTC neurons, demonstrated better
robustness compared to an end-to-end ConvNet, a fully-connected 64-neuron
long-short term memory (LSTM) [24] and a fully-connected 64-neuron CT-RNN [5].
As this paper investigates autonomous lane-keeping too, our setup includes the
two bio-inspired models from the mentioned paper: LTC neurons in a 19-neuron
NCP architecture and 64 CT-RNNs all-to-all connected.

To the best of our knowledge, this paper is the first one to systematically
explore the role of chemical and electrical synapses, in the same types of
architectures. This combination enables us to provide a response to how much
the architectural or synaptic model, respectively, has an impact on the
network’s accuracy, robustness, and interpretability.

##  IV METHODS

Figure 1: The 4-layer NCP stacked with a convolutional head. The ConvNet
extracts information from the input image which then serves as sensory input
for the NCP. In a fully connected architecture, we replace the NCP with all-
to-all wirings.

We thoroughly analyze the impact of the architecture and the synaptic models
by conducting two types of experiments:

####  IV-1 Fix the network architecture

In other words, change the synaptic model while keeping the number of neurons
and their connections. This may result in an unfair performance comparison,
because for the same network architecture, the simpler synaptic model used by
CT-RNNs, offers much less trainable parameters. At the same time, the activity
of the corresponding neurons is straightforward to compare.

####  IV-2 Fix the number of network parameters

In other words, change the network architecture, such that it results in the
same number of trainable parameters. This way, the networks have an equal
learning capacity. However, as CT-RNNs require a much bigger network
architecture, they may lose the nice properties of interpretability and
compactness. There is a plain trade-off between using the same network
architecture, or the same number of trainable parameters.

As network bases, we use the LTC-NCP with 19 neurons, and the fully-connected
version of CT-RNN with 64 neurons, as in [5]. To close the knowledge gap
between these two, we analyze the CT-RNNs in NCP architecture, with two sizes:
19 neurons, which corresponds to the type (1) experiment, and 64 neurons,
which corresponds to the type (2) experiment, with approximately 1800
trainable parameters. Additionally, we include the LTCs in a fully-connected
architecture, for a thorough analysis of the effect of the architecture. These
setups are summarized in Table I.

###  IV-A Evaluation Metrics

We use several evaluation metrics to assess driving performance. They include
prediction error, driving smoothness, and crash likelihood in closed-loop
simulation on unknown trajectories. For interpretability, the metrics of
interest are the network attention, and the changes in neuron activity during
driving. In order to evaluate network robustness, in noisy conditions not
present during training, we compare the attention maps and the driving
ability, using the VISTA simulator. Our evaluation seeks to assess the
importance of the synaptic model and the network architecture.

TABLE I: Parameters of examined networks. The 19-neuron NCP has 12 interneurons, 6 command neurons and 1 motor neuron. The 64-neuron NCP is proportionally larger, with 42 interneurons, 21 command neurons and 1 motor neuron. ID |  Neural model |  Network arch. | #Neuron | #Synapse | #Param
---|---|---|---|---|---
1 | LTC | NCP | 19 | 444 | 1833
2 | LTC | fully | 19 | 1577 | 6365
3 | CT-RNN | NCP | 19 | 444 | 482
4 | CT-RNN | NCP | 64 | 1700 | 1828
5 | CT-RNN | fully | 64 | 8192 | 8320

####  IV-A1 Training and Validation Loss

We calculate the mean squared error (MSE) between the predicted road curvature
$\hat{y_{i}}$ and the ground-truth $y_{i}$ of the validation set. Previous
work [25] pointed out that MSE alone, may not show the full driving quality.
We thus use a second metric, weighting MSE values with the road steepness, at
each time step.

####  IV-A2 Crash-Likelihood

We aim to assess the performance of the networks in a closed-loop setting,
where the next camera input is influenced by the previous steering decision.
To this end, we measure the crash likelihood during multiple test runs. In the
simulator, a crash is counted, if the car leaves the road. To test the models
in even more extreme conditions, we add extra Gaussian noise to the camera
input, which was not present as a data augmentation technique in the training
process. First, we use mean $\mu\,{=}\,0$ and variance $\sigma^{2}\,{=}\,0.1$,
to see how they tackle this condition, and then we run experiments with
doubled noise variance $\sigma^{2}\,{=}\,0.2$.

####  IV-A3 Trajectory Smoothness

To measure the smoothness of the driving trajectory, we compute the Lipschitz
constant of the prediction over a test episode. A function $f$ is called
Lipschitz continuous with Lipschitz constant $L\,{\geq}\,0$ if it satisfies:
$|f(x_{1})\,{-}\,f(x_{2})|\,{\leq}\,L\,|x_{1}\,{-}\,x_{2}|$. If
$|x_{1}\,{-}\,x_{2}|$ is the time passed between two timesteps, then
$|f(x_{1})\,{-}\,f(x_{2})|$ is the corresponding change in the prediction. The
value of $L$ measures the maximum rate of change in the network’s output. Thus
smaller values indicate smoother driving.

####  IV-A4 Saliency Maps and Structural-Similarity Index

Getting insights into how a neural network makes its decision is very
important in end-to-end autonomous systems. A possible approach is to analyze
where the neural network pays the most attention during decision-making, that
is, find out which are the pixels in the input image that have the greatest
impact. We observe the distinct evolution of the identical ConvNet
architecture across all the models after training. To visualize the network’s
attention, we utilize the VisualBackprop algorithm [26], on our ConvNet heads.
The main idea behind this method is, that the most relevant information comes
from the feature maps of the deepest layers, which we bring backward layer by
layer. This goes repeatedly, until the first layer, after which we take the
absolute value of the attention, because the image standardization layer
shifted the pixel values, and both positive and negative values should be
considered with equal importance. To make one extra step, we aim to see the
robustness of their learned ConvNet heads, in the presence of input noise, by
comparing their saliency maps under normal and additional noise pairwise. The
structural-similarity index (SSIM) [27], compares luminance, contrast, and
structure of two signals, respectively. Signals are in our context images. In
Figure 2 we show how to compute the SSIM metric.

Figure 2: Structural-similarity index (SSIM) between saliency maps, under
normal and Gaussian-noise conditions.

####  IV-A5 Neural Activity

Another way to interpret, that is, to get insights into the network behavior,
besides saliency maps, is to examine neural activity over time. Particularly
focusing on whether the neural network’s individual neurons exhibit correlated
activity with straight road segments, and turns, along the driving path. This
is achieved by extracting the driving trajectory from the simulator, as seen
from a bird’s-eye view, and then projecting the output of each neuron onto the
car’s trajectory, separately. To measure the similarity between the road
curvature and the neural activity, we compute the cross-correlation of these
signals. As we are interested in both negative and positive correlations,
compared to no correlation, we take the absolute value of the correlations,
and average them across all neurons in the network.

###  IV-B Experimental Setup

We evaluate the performance of the NCP and all-to-all architectures, with
either electrical or chemical synapses, respectively, within an end-to-end
autonomous driving scenario. Here the agent receives the front camera input of
the car, and predicts the road curvature, which corresponds to the steering
action of the lane-keeping task. Our dataset consists of human driving
recordings, under different weather conditions [28]. Figure 3 shows samples
from this dataset.

Figure 3: Examples from the dataset: summer and winter

Predicting the road curvature instead of directly predicting the steering
angle is beneficial because the curvature value is independent from the
vehicle, while the steering angle is not. The steering angle can be computed
from the curvature:

$\alpha_{t}=S_{v}\>\textrm{arctan}(L_{v}\>y_{t})$ (7)

where $\alpha_{t}$ is the steering angle at timestep $t$, $S_{v}$ and $L_{v}$
are the steering-ratio and wheelbase properties of the vehicle, respectively,
and $y_{t}$ is the road curvature at timestep $t$.

Figure 4: Training and testing procedure.

For training purposes, we create an open-loop imitation learning task, from
summer and winter data. The training data contains input-output sequences
$\\{(x_{t},y_{t})\,|\,t\,{\in}[0..T]\\}$, where $x_{t}$ is the camera frame,
and $y_{t}$ is the road curvature at time $t$. The network receives the camera
frame $x_{t}$ at time step $t$ and predicts the road curvature $\hat{y}_{t}$.

As data augmentation, we use the following procedure. We first adjust the
brightness, by adding to the pixels a random value chosen from a uniform
distribution in $[-0.4,0.4]$. After this, the contrast and saturation of the
image are shifted, by multiplying each RGB channel by a contrast factor, and
the saturation channel in its HSV representation by a saturation factor,
respectively. These random values are taken from a uniformly distribution in
$[0.6,1.4]$. In addition to this, the training data is augmented by Guided
Policy Learning [29] using the VISTA simulator [30, 11]. This places the car
off-orientation on the track, and the model has to learn to recover the car
from these quite challenging positions.

The model is trained using the mean squared error (MSE) of the prediction
$\hat{y}_{t}$ and the label $y_{t}$. This is offline learning: the prediction
of the network does not influence the next input image of the sequence. This
can be seen on the left-hand side in Figure 4. We run the training on a
sequence of input images of length 32, with a batch size of 32, with a cosine
learning-rate decay schedule with an initial value of $10^{-3}$ over 100
epochs. We apply the AdamW optimizer [31] which in our case decouples the
weight decay with the value of $10^{-6}$ from the learning rate. The best
epoch is determined by the highest validation loss, without early stopping, in
order to achieve a better generalization [32].

After training the models, we test them in the data-driven VISTA simulation
platform, which synthesizes new possible trajectories from the dataset. This
way we can evaluate the generalization performance of the models, in a closed-
loop environment, where the prediction of the current time step affects the
camera input of the next time step. This approach is illustrated on the right-
hand side in Figure 4.

##  V RESULTS

###  V-A MSE in Open-Loop

Table II compares the neural networks in terms of MSE and Weighted MSE in the
open-loop setting. We observe that the CT-RNN NCP-64 model exhibited the
lowest MSE error, indicating the best fit to the data compared to the other
models. By assigning higher weights to predictions in steeper curves, LTC-
based models showed better performance to it by achieving the lowest Weighted
MSE values.

TABLE II: Prediction error in the validation set, for 5 trained models, per model-type. Smaller error is better performance. Model | MSE | Weighted MSE
---|---|---
LTC NCP-19 | 0.173 $\pm$ 0.021 | 0.011 $\pm$ 0.001
LTC fully-19 | 0.200 $\pm$ 0.027 | 0.010 $\pm$ 0.001
CT-RNN NCP-19 | 0.174 $\pm$ 0.012 | 0.015 $\pm$ 0.002
CT-RNN NCP-64 | 0.163 $\pm$ 0.017 | 0.012 $\pm$ 0.002
CT-RNN fully-64 | 0.219 $\pm$ 0.015 | 0.018 $\pm$ 0.002

###  V-B Crash-Likelihood in Closed-Loop

We tested the networks for 125 episodes for each season and calculated how
likely they drive off the road over these episodes. We found that the
chemical-synapse-based LTCs independently of the architecture, and the
electrical-synapse-based CT-RNNs in NCP wiring architectures drive without
crash in the closed-loop setting using the VISTA simulator. Moreover, they are
robust enough, to tackle the challenge of receiving input images with 0-mean
and $0.1$-variance Gaussian noise. In contrast, the fully-connected CT-RNN
network crashes without any input disturbances during winter and under
additional noise in both summer and winter tests.

We also aimed to push these models to their limits, and observed that by
doubling the noise variance, the models are not able to avoid the crashes
anymore. The trend of the crash likelihood shows that the 4-layer structure of
the NCP brings an advantage compared to the fully-connected variant.

Note that some of the evaluation metrics of the CT-RNN in fully-connected
architecture could be not relevant, as this was not able to drive without
crashes under normal circumstances during winter, and under 0.1-variance
Gaussian noise in both seasons. However, we included them in our full
analysis.

TABLE III: Crash likelihood in closed-loop simulation, under noise-free conditions, and under Gaussian noise with mean of 0 and variance of 0.1 and 0.2. Smaller is better. S is summer, W is winter. 125 runs per model in each condition. | Noise-free | $\sigma^{2}=0.1$ | $\sigma^{2}=0.2$
---|---|---|---
Model | S | W | S | W | S | W
LTC NCP-19 | 0.0 | 0.0 | 0.0 | 0.0 | 0.11 | 0.01
LTC fully-19 | 0.0 | 0.0 | 0.0 | 0.0 | 0.18 | 0.02
CT-RNN NCP-19 | 0.0 | 0.0 | 0.0 | 0.0 | 0.20 | 0.08
CT-RNN NCP-64 | 0.0 | 0.0 | 0.0 | 0.0 | 0.05 | 0.03
CT-RNN fully-64 | 0.0 | 0.02 | 0.09 | 0.06 | 0.68 | 0.25

###  V-C Trajectory Smoothness

Table IV provides the Lipschitz constants for networks in summer and winter
contexts, both in noise-free conditions and under the influence of 0-mean and
0.1-variance noise. Observe that LTC networks generally exhibit smoother
driving performance compared to CT-RNNs, in most scenarios. Moreover, when
employing the NCP architecture, LTCs achieve an exceptionally smooth
trajectory.

TABLE IV: Lipschitz constant of the prediction in the first 4 columns. Smaller is better. Similarity of the trajectory between normal conditions and under Gaussian noise variance of 0.1 in the last 2 columns. Higher is better. LTC-based networks can maintain a smooth driving trajectory in the presence of noise better than CT-RNN networks. Using an NCP architecture makes the neural network more resilient. Note that even though the CT-RNN fully-64 results show good characteristics, this model drives off the road under this amount of noise, while the others can drive without crashes. | Lipschitz constant of the output | Similarity of the trajectory
---|---|---
| Noise-free | Under noise $N(\mu=0,\sigma^{2}=0.1)$ | Under noise $N(\mu=0,\sigma^{2}=0.1)$
Model | Summer | Winter | Summer | Winter | Summer | Winter
LTC NCP-19 | 0.0081 $\pm$ 0.0028 | 0.0087 $\pm$ 0.0031 | 0.0115 $\pm$ 0.0028 | 0.0113 $\pm$ 0.0036 | 0.971 $\pm$ 0.043 | 0.980 $\pm$0.027
LTC fully-19 | 0.0076 $\pm$ 0.0022 | 0.0101 $\pm$ 0.0045 | 0.0124 $\pm$ 0.0028 | 0.0126 $\pm$ 0.0041 | 0.967 $\pm$ 0.046 | 0.976 $\pm$ 0.030
CT-RNN NCP-19 | 0.0107 $\pm$ 0.0097 | 0.0150 $\pm$ 0.0098 | 0.0143 $\pm$ 0.0093 | 0.0172 $\pm$ 0.0092 | 0.964 $\pm$ 0.054 | 0.968 $\pm$ 0.047
CT-RNN NCP-64 | 0.0109 $\pm$ 0.0067 | 0.0124 $\pm$ 0.0070 | 0.0124 $\pm$ 0.0060 | 0.0148 $\pm$ 0.0070 | 0.962 $\pm$ 0.057 | 0.967 $\pm$ 0.047
CT-RNN fully-64 | 0.0076 $\pm$ 0.0029 | 0.0107 $\pm$ 0.0089 | 0.0126 $\pm$ 0.0033 | 0.0122 $\pm$ 0.0049 | 0.911 $\pm$ 0.171 | 0.970 $\pm$ 0.040

###  V-D Network Interpretability

####  V-D1 Saliency maps

Figure 5 shows the saliency maps of the networks considered, for both summer
and winter. We found that the networks extract information differently in the
two seasons: the attention goes to the road itself in summer, and the contour
(the side) of the road in winter. Intuitively, the areas of the input image
containing the most relevant steering information should be highlighted: the
road, and its contours, without the unnecessary environmental features which
occur outside of the road. Note that the NCPs focus on the relevant parts,
alone, while the fully connected networks also consider irrelevant areas of
the input frames.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5: Summer-winter input images in (a), and their saliency maps of LTC
NCP-19 in (b), LTC fully-19 in (c), CT-RNN NCP-19 in (d), CT-RNN NCP-64 in (e)
and CT-RNN fully-64 in (f). Observe that an NCP architecture makes the models
more focused compared to their fully connected version. The LTC NCP-19 keeps
its attention on the road and its contours in summer, and on the road contours
in winter. It ignores the irrelevant information from the border of the image.

####  V-D2 Neural Activity

When considering the neural activity during driving, we observed that in NCP
architectures, neural activities change their dynamics based on the road
curvature, regardless of the synaptic type. This behavior holds for LTCs in
the fully-connected architecture, too, which means that this interpretability
feature holds regardless of the network architecture, for chemical synapses.
However, the electrical-synapse-based CT-RNN does not have this nice property
without the NCP structure. Its all-to-all wiring shows no patterns of
different activation during turns and straight sections of the road. For
illustration, Figure 6 shows the activity of one of the neurons in the
examined networks.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: Neural activity. (a) Road trajectory from above, arrows indicate the
driving direction (b)-(f) Models in the same order as before. The NCP
architecture brings interpretability to CT-RNNs as the changes in neural
activity correspond to the road’s trajectory, as illustrated in (d)-(e). LTC
in all-to-all architecture is still interpretable in this sense, as seen in
(c).

In addition, we also computed the average correlation between the neural
activity of each neuron, and the driving trajectory. The results in Table V
clearly show that the LTC neurons have the best interpretability in terms of
network dynamics, irrespective of the architecture used. In contrast, for the
classic CT-RNN models, the NCP architecture gets around 1.5-2 times increased
activity correlation compared to the all-to-all one. Hence, the NCP
architecture brings CT-RNNs very close to the chemical-synapse-based LTC
model.

TABLE V: Correlation between the neural activity and the driving trajectory. Higher is better in this case. The NCP architecture brings the CT-RNN model much closer to the LTC-type networks, in terms of interpretable neural activity. Model | Summer | Winter
---|---|---
LTC NCP-19 | 0.662 $\pm$ 0.290 | 0.636 $\pm$ 0.279
LTC fully-19 | 0.661 $\pm$ 0.307 | 0.698 $\pm$ 0.270
CT-RNN NCP-19 | 0.649 $\pm$ 0.316 | 0.630 $\pm$ 0.293
CT-RNN NCP-64 | 0.568 $\pm$ 0.313 | 0.524 $\pm$ 0.282
CT-RNN fully-64 | 0.395 $\pm$ 0.282 | 0.321 $\pm$ 0.254

###  V-E Network Robustness

Finally, we measured how much the saliency maps of the various networks,
changed in the presence of a noisy input. An SSIM value close to 1 indicates a
more robust attention mechanism, meaning that the additional noise has less of
an influence on the convolutional head of these networks. In Figure 7 we
present our results. These show that the LTC model in the NCP wiring
architecture is the most robust: it extracts nearly identical information from
the input images even when noise is present.

Figure 7: The structural-similarity index (SSIM) for the neural networks
considered, under Gaussian noise with mean 0, and variance 0.1 and 0.2, in
summer on the left-hand side and winter on the right-hand side. Higher is
better.

In order to evaluate how well the five models considered can maintain their
driving trajectory, we also calculated the similarity between their
trajectories with and without noise. These results are reported and discussed
in Table IV.

##  VI CONCLUSIONS

In this paper, we showed how neuroscience models of electrical and chemical
synapses result in continuous-time recurrent neural networks. These bio-
inspired neural networks hold the potential for interpretable solutions, which
are highly sought after in robotics systems.

We investigated the impact of synaptic choices in the architectures of sparse
and fully-connected wirings. Through extensive experiments conducted in a
photorealistic autonomous driving simulator, we evaluated their performance,
interpretability, and robustness under diverse conditions and in the presence
of noise.

Our findings revealed that the choice of wiring architecture had a significant
impact on the networks’ performance. Specifically, the sparse NCP wirings were
more effective with both synaptic models. Concerning the synaptic models
themselves, our research demonstrated that learning with chemical synapses
makes a difference. They not only led to a more robust behavior but also
resulted in better interpretability and significantly smoother driving
dynamics when compared to electrical synapses.

## References

  * [1] R. D. Beer and J. C. Gallagher, “Evolving dynamical neural networks for adaptive behavior,” Adaptive behavior, vol. 1, no. 1, pp. 91–122, 1992.
  * [2] K. Funahashi and Y. Nakamura, “Approximation of dynamical systems by continuous time recurrent neural networks,” Neural Networks, vol. 6, pp. 801–806, 1993.
  * [3] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, “Neural ordinary differential equations,” in Advances in Neural Information Processing Systems, vol. 31, pp. 1–13, Curran Associates, Inc., 2018.
  * [4] R. M. Hasani, M. Lechner, A. Amini, D. Rus, and R. Grosu, “Liquid time-constant networks,” in AAAI Conference on Artificial Intelligence, 2020.
  * [5] M. Lechner, R. M. Hasani, A. Amini, T. A. Henzinger, D. Rus, and R. Grosu, “Neural circuit policies enabling auditable autonomy,” Nature Machine Intelligence, vol. 2, pp. 642–652, 2020.
  * [6] C. J. Vorbach, R. M. Hasani, A. Amini, M. Lechner, and D. Rus, “Causal navigation by continuous-time neural networks,” in Neural Information Processing Systems, 2021.
  * [7] M. Chahine, R. Hasani, P. Kao, A. Ray, R. Shubert, M. Lechner, A. Amini, and D. Rus, “Robust flight navigation out of distribution with liquid neural networks,” Science Robotics, vol. 8, no. 77, p. eadc8892, 2023.
  * [8] M. Lechner, R. Hasani, M. Zimmer, T. A. Henzinger, and R. Grosu, “Designing worm-inspired neural networks for interpretable robotic control,” in 2019 International Conference on Robotics and Automation (ICRA), pp. 87–94, 2019.
  * [9] G.-B. Zhou, J. Wu, C.-L. Zhang, and Z.-H. Zhou, “Minimal gated unit for recurrent neural networks,” 2016.
  * [10] S. Schaal, “Is imitation learning the route to humanoid robots?,” Trends in Cognitive Sciences, vol. 3, pp. 233–242, 1999.
  * [11] A. Amini, T.-H. Wang, I. Gilitschenski, W. Schwarting, Z. Liu, S. Han, S. Karaman, and D. Rus, “Vista 2.0: An open, data-driven simulator for multimodal sensing and policy learning for autonomous vehicles,” in 2022 International Conference on Robotics and Automation (ICRA), pp. 2419–2426, IEEE, 2022.
  * [12] S. R. Wicks, C. J. Roehrig, and C. H. Rankin, “A dynamic network simulation of the nematode tap withdrawal circuit: predictions concerning synaptic function using behavioral criteria,” Journal of Neuroscience, vol. 16, no. 12, pp. 4017–4031, 1996.
  * [13] E. R. Kandel, J. H. Schwartz, T. M. Jessell, S. Siegelbaum, A. J. Hudspeth, S. Mack, et al., Principles of neural science, vol. 4.  McGraw-hill New York, 2000.
  * [14] E. Johns, “Coarse-to-fine imitation learning: Robot manipulation from a single demonstration,” in 2021 IEEE international conference on robotics and automation (ICRA), pp. 4613–4619, IEEE, 2021.
  * [15] X. B. Peng, E. Coumans, T. Zhang, T.-W. Lee, J. Tan, and S. Levine, “Learning agile robotic locomotion skills by imitating animals,” arXiv preprint arXiv:2004.00784, 2020.
  * [16] H. Karnan, G. Warnell, X. Xiao, and P. Stone, “Voila: Visual-observation-only imitation learning for autonomous navigation,” in 2022 International Conference on Robotics and Automation (ICRA), pp. 2497–2503, IEEE, 2022.
  * [17] Y. Fan, S. Chu, W. Zhang, R. Song, and Y. Li, “Learn by observation: Imitation learning for drone patrolling from videos of a human navigator,” in 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 5209–5216, IEEE, 2020.
  * [18] M. Bojarski, D. W. del Testa, D. Dworakowski, B. Firner, B. Flepp, P. Goyal, L. D. Jackel, M. Monfort, U. Muller, J. Zhang, X. Zhang, J. Zhao, and K. Zieba, “End to end learning for self-driving cars,” ArXiv, vol. abs/1604.07316, 2016.
  * [19] H. Xu, Y. Gao, F. Yu, and T. Darrell, “End-to-end learning of driving models from large-scale video datasets,” 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 3530–3538, 2016.
  * [20] R. D. Beer, “A dynamical systems perspective on agent-environment interaction,” Artificial intelligence, vol. 72, no. 1-2, pp. 173–215, 1995.
  * [21] H. Tang, B. H. Tan, and R. Yan, “Robot-to-human handover with obstacle avoidance via continuous time recurrent neural network,” in 2016 IEEE Congress on Evolutionary Computation (CEC), pp. 1204–1211, 2016.
  * [22] P. Hénaff, V. Scesa, F. B. Ouezdou, and O. Bruneau, “Real time implementation of ctrnn and bptt algorithm to learn on-line biped robot balance: Experiments on the standing posture,” ArXiv, vol. abs/2011.06910, 2020.
  * [23] R. Sendra-Arranz and Á. Gutiérrez, “Emergence of communication through artificial evolution in an orientation consensus task in swarm robotics,” in IFIP International Conference on Artificial Intelligence Applications and Innovations, pp. 515–526, Springer, 2023.
  * [24] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural computation, vol. 9, no. 8, pp. 1735–1780, 1997.
  * [25] F. Codevilla, A. M. Lopez, V. Koltun, and A. Dosovitskiy, “On offline evaluation of vision-based driving models,” in Proceedings of the European Conference on Computer Vision (ECCV), pp. 236–251, 2018.
  * [26] M. Bojarski, A. Choromańska, K. Choromanski, B. Firner, L. D. Jackel, U. Muller, and K. Zieba, “Visualbackprop: visualizing cnns for autonomous driving,” ArXiv, vol. abs/1611.05418, 2016.
  * [27] Z. Wang, A. Bovik, H. Sheikh, and E. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600–612, 2004.
  * [28] M. Lechner, R. Hasani, A. Amini, T.-H. Wang, T. A. Henzinger, and D. Rus, “Are all vision models created equal? a study of the open-loop to closed-loop causality gap,” arXiv preprint arXiv:2210.04303, 2022.
  * [29] S. Levine and V. Koltun, “Guided policy search,” in Proceedings of the 30th International Conference on Machine Learning (S. Dasgupta and D. McAllester, eds.), vol. 28 of Proceedings of Machine Learning Research, (Atlanta, Georgia, USA), pp. 1–9, PMLR, 17–19 Jun 2013.
  * [30] A. Amini, I. Gilitschenski, J. Phillips, J. Moseyko, R. Banerjee, S. Karaman, and D. Rus, “Learning robust control policies for end-to-end autonomous driving from data-driven simulation,” IEEE Robotics and Automation Letters, vol. 5, pp. 1143–1150, 2020.
  * [31] I. Loshchilov and F. Hutter, “Fixing weight decay regularization in adam,” ArXiv, vol. abs/1711.05101, 2017.
  * [32] A. Power, Y. Burda, H. Edwards, I. Babuschkin, and V. Misra, “Grokking: Generalization beyond overfitting on small algorithmic datasets,” ArXiv, vol. abs/2201.02177, 2022.