∎

11institutetext: Pilhwa Lee 22institutetext: Department of Mathematics, Morgan
State University
Baltimore, MD
22email<EMAIL_ADDRESS>
# Electrodiffusion with calcium-activated potassium channels in dendritic
spine

Pilhwa Lee

(Received: date / Accepted: date)

###### Abstract

We investigate calcium signaling feedback through calcium-activated potassium
channels of a dendritic spine by applying the immersed boundary method with
electrodiffusion. We simulate the stochastic gating of such ion channels and
the resulting spatial distribution of concentration, current, and membrane
voltage within the dendritic spine. In this simulation, the permeability to
ionic flow across the membrane is regulated by the amplitude of chemical
potential barriers. With spatially localized ion channels, chemical potential
barriers are locally and stochastically regulated. This regulation represents
the ion channel gating with multiple subunits, the open and closed states
governed by a continuous-time Markov process. The model simulation
recapitulates an inhibitory action on voltage-sensitive calcium channels by
the calcium-activated potassium channels in a stochastic manner as a _non-
local_ feedback loop. The model predicts amplified calcium influx with more
closely placed channel complexes, proposing a potential mechanism of
differential calcium handling by channel distributions. This work provides a
foundation for future computer simulation studies of dendritic spine motility
and structural plasticity.

###### Keywords:

the immersed boundary method electrodiffusion dendritic spine voltage-
sensitive calcium channel calcium-activated potassium channel continuous-time
Markov process

##  1 Introduction

Dendritic spines are small protrusions in the postsynapse and dendritic trees
of neurons (Shepherd 1996). They are crucial in learning and memory (Harris
and Kater 1993, Seung and Kim 2002, Yasuda et al. 2003, Matus 2004, Bloodgood
and Sabatini 2007). In developmental stages, they are highly motile in the
course of neuronal wiring and pruning (Bonhoeffer and Yuste 2002, Tada and
Seung 2006, Bosch and Hayashi 2011). However, dendritic spine motility and
potentially involving structural plasticity are persistently observed in
matured stages (Lamprecht et al. 2006) and quite evident in aging and
neurodegeneration, neuronal injury, as well as psychiatric disorders (Koleske
2013, Kim et al. 2015, Rice et al. 2015). In the presence of several ion
channels in the dendritic spine (Sabatini and Svoboda 2000, Bloodgood and
Sabatini 2007b), relevant modeling with electrodiffusion in the micro- and
nano-domains is needful to reconstruct physiological ionic transports and
understand consequential motility and associated neuronal functions (Holcman
and Yuste 2015).

Previously we have studied the electrodiffusion with stochastic voltage-
sensitive calcium channels in a one-dimensional setting using an immersed
boundary method for electrodiffusion (Lee et al. 2013). The interface
conditions for membrane permeability to each ionic species are replaced by
chemical potential barriers in a unified Cartesian domain without explicit
dissection of the computational domain between intracellular and extracellular
domains. This article extends this approach to a two-dimensional domain
considering the spatial effects of stochastic channel gating and channel
distribution. We focus on voltage-sensitive calcium channels and calcium-
activated potassium channels. In the configuration of the spatial expression
of those two types of channels, they are mostly coupled together as “complex”
with about 15 nm distance in the nanodomain (Berkefeld et al. 2006). There are
two types of calcium-activated potassium channels; small-conductance (SK) and
large-conductance (BK) (Sah and Faber 2002). These potassium channels give
feedback to the voltage-sensitive calcium channels with
repolarization/hyperpolarization to inhibit further calcium influx (Griffith
et al. 2016) and regulate in the synaptic strength and short-term plasticity
among many other roles (He et al. 2014, Lujan et al. 2018, Gutzmann et al.
2019).

The chemical potential barriers are chosen so that the membrane is permeable
to Ca2+ and K+ locally in space when the voltage-sensitive calcium channels
and the calcium-activated potassium channels are open. To recapitulate the
back-propagating action potential from the soma, the membrane is also
regulated uniformly in space to be semi-permeable to Na+ for the membrane
depolarization. The membrane is uniformly semi-permeable to Cl- all the time
of the simulation.

We place five voltage-sensitive calcium channels and five calcium-activated
potassium channels (BK) in the dendritic spine head. According to continuous-
time Markov processes, the opening and closing of the channels are modeled by
lowering and raising the heights of the chemical potential barriers of calcium
and potassium. The voltage-sensitive calcium channels have four independent
subunits with inactivation from the intracellular local calcium concentration
(Cox 2014). The channels are open only when all four subunits are in the open
state. Voltage-sensitive rate constants govern the transitions between the
open and closed states of the subunits, and the transition to the inactivated
state is governed by the intracellular calcium concentration. The calcium-
activated potassium channels also have four independent subunits with closed
and open states dependent on the membrane voltage, and each subunit is
activated by the intracellular calcium (Cox 2014).

As a whole, we have incorporated _an immersed boundary method with
electrodiffusion_ for ionic transport in the dendritic spine of synapse on the
two-dimensional domain, coupled to continuous-time Markov processes for the
ion channel gating of voltage-sensitive calcium channels and calcium-activated
potassium channels, possibly distributed non-uniformly in the dendritic spine.
The model simulation recapitulates an inhibitory action on voltage-sensitive
calcium channels by calcium-activated potassium channels in a stochastic
manner as a _non-local_ feedback loop. The model also predicts enhanced
calcium influx with more closely placed channel complexes, proposing a
mechanism of differential calcium handling by the _locality_ of channel
distribution.

The paper is organized in the following way; in Section 2, we describe the
mathematical formulation of the electrodiffusion of ion species in the
immersed boundary formalism, the ion-channel gating as a continuous-time
Markov process, and the resulting regulation of the chemical potential
barriers that model ion-channel selectivity. In Section 3, we present a two-
dimensional study with spatially localized channels. In Section 4, we summary
the model results and mention future works.

##  2 MATHEMATICAL FORMULATION

In this section, we consider a fixed two-dimensional computational domain with
dissolved ions. Immersed within the domain is a closed membrane, which is
fixed in place. The membrane may be permeable or impermeable to each ionic
species, the permeability being controlled in a graded manner by its chemical
potential barrier. We have the following notations:

$D_{i}$: diffusion coefficient of the $i^{\rm th}$ ion species
$q$: the unit electrical charge (charge on a proton)
$qz_{i}$: charge of the $i^{\rm th}$ species
$K_{\rm B}$: Boltzmann constant
$T$: absolute temperature (degrees in Kelvin)
$\Omega_{\rm E}$: Eulerian domain
$\Omega_{\rm L}$: Lagrangian domain
$\varepsilon$: dielectric constant

The notations for the variables are the following:
$\mathbf{x}=(x_{1},x_{2})$: Eulerian coordinate
$\mathbf{x}=\mathbf{X}(s)$: Lagrangian description
$\psi_{i}(\mathbf{x},t)$: chemical potential of the $i^{\rm th}$ ion species
$\Psi(\mathbf{x})$: chemical potential kernel
$A_{i}(s,t)ds$: contribution of arc $(s,s+ds)$ of the membrane to the chemical
potential of $i^{\rm th}$ ion species
$c_{i}(\mathbf{x},t)$: concentration of the $i^{\rm th}$ ion species
$\mathbf{J}_{i}(\mathbf{x},t)$: flux per unit area of the $i^{\rm th}$ ion
species
$\varphi(\mathbf{x},t)$: electrical potential
$\rho_{0}(\mathbf{x})$: background electrical charge density

###  2.1 The chemical potential

The chemical potential is expressed in the following way:

$\psi_{i}(\mathbf{x},t)=\int_{\Omega_{\rm
L}}\Psi(\mathbf{x}-\mathbf{X}(s),\mathbf{t},\mathbf{n})A_{i}(s,t)ds.$ (1)

Here, $\mathbf{X}(s)$ is the configuration of the immersed boundary, where $s$
is a Lagrangian parameter. The function $A_{i}(s,t)$ describes the
contribution of the membrane at $\mathbf{X}(s)$ to the chemical potential
barrier for the $i^{\rm th}$ ionic species. The chemical potential kernel
$\Psi$ defines how the contribution $A_{i}(s,t)ds$ is to be spread out in
space in the neighborhood of $\mathbf{X}(s)$. In regulating membrane
permeability to the $i^{\rm{th}}$ species of ion, the chemical potential
amplitude $A_{i}(s,t)$ is modulated locally or globally. When it is locally
controlled, the localized area represents the domain of ion channels for the
$i^{\rm{th}}$ ionic species.

In general, any bell-shaped function with compact support can be available for
the chemical potential kernel. For the one-dimensional kernel, a smoothed
Dirac delta function $\phi$ of the second-order moment with compact support is
used following the function constructed by Peskin peskin . The chemical
potential kernel on the two-dimensional domain is as follows:

$\Psi_{w}(\mathbf{x},\mathbf{t},\mathbf{n})=\frac{1}{w^{2}}\phi(\frac{\mathbf{x}\cdot\mathbf{t}}{w})\phi(\frac{\mathbf{x}\cdot\mathbf{n}}{w}),$
(2)

where $w$ is a scaling factor such that $\Psi_{w}$ has a support of a square
of edge $4w$. The coordinates used in the $\phi$ are in the local frame from
the tangential and normal directions, $\mathbf{t}$ and $\mathbf{n}$ concerning
the membrane at $\mathbf{X}(s)$.

###  2.2 The electrostatic potential and the electrodiffusion: Poisson-Nernst-
Planck equations

The electrical potential is a solution of the Poisson equation:

$\displaystyle-\nabla\cdot(\epsilon\nabla\varphi)=\sum_{i}qz_{i}c_{i}+\rho_{0},$
(3)

where $\rho_{0}$ represents the background electrical charge density. The two-
dimensional domain is prescribed to be periodic in each direction, and the
necessary condition for the existence of the solution of the Poisson equation
requires the global electroneutrality. In the immersed boundary method with
electrodiffusion, local electroneutrality is also satisfied except for the
space charge layer around the membrane (Lee 2007, Lee et al. 2010).

The electrodiffusion equations are formulated in the following way:

$\displaystyle\frac{\partial c_{i}}{\partial t}+\nabla\cdot\mathbf{J}_{i}$
$\displaystyle=$ $\displaystyle 0,$ (4) $\displaystyle\mathbf{J}_{i}$
$\displaystyle=$ $\displaystyle-D_{i}(\nabla
c_{i}+c_{i}\frac{\nabla(\psi_{i}+qz_{i}\varphi)}{K_{\rm B}T}).$ (5)

Eq.(4) is the conservation law (the continuity equation) for the $i^{\rm th}$
ionic species. In this equation, $c_{i}$ is the concentration, and
$\mathbf{J}_{i}$ is the flux per unit area of this ion species. Eq.(5) gives
the flux per unit area as a sum of three terms: diffusion, drift caused by the
chemical potential, and drift caused by the electrical potential.

###  2.3 Continuous-time Markov process for the stochastic ion channel gating

We will reformulate the kinetic models of voltage-sensitive calcium channels
and calcium-activated potassium channels in Cox 2014 as continuous-time Markov
processes to represent their stochastic channel gating activities on the
spatial domain. For the ion channel with the center at $\mathbf{X}(s)$, the
involving membrane voltage, $V_{\rm m}(s)$, the intracellular concentration,
$c_{i}(s)$, and ionic current, $I_{i}(s)$ are collected from the following:

$\displaystyle V_{\rm m}(s)$ $\displaystyle=$ $\displaystyle\int_{\Omega_{\rm
E}}(\Psi_{w}(\mathbf{x}-\mathbf{X}(s)+2w\mathbf{n})-\Psi_{h}(\mathbf{x}-\mathbf{X}(s)-2w\mathbf{n}))\varphi(\mathbf{x})d\mathbf{x},$
(6) $\displaystyle c_{i}(s)$ $\displaystyle=$ $\displaystyle\int_{\Omega_{\rm
E}}\Psi_{w}(\mathbf{x}-\mathbf{X}(s)+2w\mathbf{n})c_{i}(\mathbf{x})d\mathbf{x},$
(7) $\displaystyle I_{i}(s)$ $\displaystyle=$ $\displaystyle
qz_{i}\int_{\Omega_{\rm
E}}\Psi_{w}(\mathbf{x}-\mathbf{X}(s)+2w\mathbf{n})\mathbf{J}_{i}(\mathbf{x})\cdot\mathbf{n}d\mathbf{x},$
(8)

where $4w$ is mostly the width of the membrane, and $\mathbf{n}$ is the normal
unit vector at $\mathbf{X}(s)$ outward from the membrane.

####  2.3.1 The voltage-sensitive calcium channel gating

The transition between closed and open states of each subunit is expressed as
follows:

$\lx@xy@svg{\hbox{\raise 0.0pt\hbox{\kern
8.16667pt\hbox{\ignorespaces\ignorespaces\ignorespaces\hbox{\vtop{\kern
0.0pt\offinterlineskip\halign{\entry@#!@&&\entry<EMAIL_ADDRESS>0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{{\rm
C_{S}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern
7.85826pt\raise 10.80554pt\hbox{{}\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern
0.0pt\raise-1.75pt\hbox{$\scriptstyle{\alpha(V_{\rm m})}$}}}\kern
3.0pt}}}}}}\ignorespaces{\hbox{\kern 32.16667pt\raise
4.30554pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern
32.16667pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern
3.0pt\raise 0.0pt\hbox{$\textstyle{{\rm
O_{S}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern
8.11754pt\raise-10.80554pt\hbox{{}\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern
0.0pt\raise-1.75pt\hbox{$\scriptstyle{\beta(V_{\rm m})}$}}}\kern
3.0pt}}}}}}\ignorespaces{\hbox{\kern
8.16667pt\raise-4.30554pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces}}}}\ignorespaces,$
(9)

where $\rm C_{S}$ and $\rm O_{S}$ represent closed and open states of a
subunit. The rate constants of ${\alpha}$ and ${\beta}$ are functions of
membrane voltage $V_{\rm m}$.

We describe the states of subunits and an ion channel in the on/off way:

$\displaystyle\chi_{i}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}1&{\rm subunit~{}open}\\\ 0&{\rm
subunit~{}closed}\end{array}\right.$ (12) $\displaystyle S$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}{\rm
C}_{i-1}&i-1{\rm~{}subunits~{}are~{}open}\\\ {\rm O}&{\rm
ion~{}channel~{}open}\\\ {\rm I}&{\rm inactivated}~{}{\rm
state}\end{array}\right.$ (16)

where $\chi_{i}$ indicates the open/closed state of $i^{\rm th}$ subunit, and
$S$ the state of the ion channel. When we express the transition probability
between the open and closed states of each subunit based on Eq.(9),

$\displaystyle P(\chi_{i}(t+dt)=1|\chi_{i}(t)=0)=\alpha(V_{m})dt,$ (17)
$\displaystyle P(\chi_{i}(t+dt)=0|\chi_{i}(t)=1)=\beta(V_{m})dt,$ (18)
$\displaystyle\alpha(V_{m})=\alpha_{0}e^{q_{\rm forward}V_{\rm m}/26.7},$ (19)
$\displaystyle\beta(V_{m})=\beta_{0}e^{-q_{\rm backward}V_{\rm m}/26.7},$ (20)

where $\alpha_{0}$ = 3.0 ms-1, $\beta_{0}$ = 0.241 ms-1, $q_{\rm
forward}=1.16$, and $q_{\rm backward}=1.94$, where $V_{\rm m}$ is in the unit
of mV. Eq.(17) represents the probability of $i^{\rm th}$ subunit to take the
transition from the closed state at $t$ to the open state at $t+dt$ in the
infinitesimal time interval $dt$. For the individual ion channel gating, a
continuous-time Markov process is applied peskin2 . The ion channel is assumed
to have four independent subunits; each has open and closed states.

The channel is open only when all four subunits are in the open state, and
when the channel is not inactivated. The diagram for the Markov process with
the discrete states is presented in Fig.(2). In the discrete states labeled
$C_{i}$, $i-1$ represents the number of subunits in the open state. The state
with all four subunits open, however, is given the special symbol O. The ion
channel is open when all those subunits are open, i.e., when it is in the
state O. In the state of I, the channel is inactivated and closed.

####  2.3.2 The calcium-activated potassium channel gating

The channel is also constituted with four independent subunits (Fig. 3). Each
subunit is activated by calcium, with binding and unbinding rate constants
$K_{\rm o}[{\rm Ca}^{2+}]$ and $K_{\rm-o}$ in the open state, and $K_{\rm
c}[{\rm Ca}^{2+}]$ and $K_{\rm-c}$ in the closed state. The constants $K_{\rm
o}$, $K_{\rm-o}$, $K_{\rm c}$, and $K_{\rm-c}$ are 1.0 nM-1s-1, 1.065 ms-1,
1.0 nM-1s-1, 11.917 ms-1, respectively. With $i$ subunits activated by
calcium, the open and closed states are determined by the forward and backward
rate constants, $K_{i}$ and $K_{-i}$:

$\displaystyle K_{i}(V_{m})=\alpha_{0,i}e^{q_{\rm forward}V_{\rm m}/26.7},$
(21) $\displaystyle K_{-i}(V_{m})=\beta_{0,i}e^{-q_{\rm backward}V_{\rm
m}/26.7},$ (22)

where $\alpha_{0,0}$ through $\alpha_{0,4}$ are 5.5, 8,0, 2.0, 884, 900 s-1,
and $\beta_{0,0}$ through $\beta_{0,4}$ are 8.669, 1.127, 0.0252, 1.013,
0.1257 ms-1, respectively. The rate constants $q_{\rm forward}$ and $q_{\rm
backward}$ are 1.16 and 1.94.

###  2.4 Regulation of the chemical potentials

As described in the previous section, the continuous-time Markov process for
the voltage-sensitive calcium channel is applied with the membrane voltage and
the intracellular calcium concentration adjacent to each ion channel. The
chemical potential for Ca2+ is modulated in the on/off way with the channel
state variable $S$ from the Markov process. In the two-dimensional scheme, the
chemical potential of calcium is _locally_ regulated. Let the width of the ion
channel be $w_{\rm ch}$. For the $j^{\rm th}$ ion channel in the state
$S_{j}$, with the Lagrangian parameter $s_{j}$ for the location of the center
of the channel, we make spatial regulation of the chemical potential in the
following way:

$A_{{\rm Ca}^{2+}}(s,t)=\left\\{\begin{array}[]{ll}A_{{\rm Ca}^{2+},{\rm
open}}&S_{j}={\rm O},s\in[s_{j}-w_{\rm ch}/2,s_{j}+w_{\rm ch}/2]\\\ A_{{\rm
Ca}^{2+},{\rm closed}}&{\rm otherwise}\end{array}\right.$ (23)

where $A_{{\rm Ca}^{2+},{\rm closed}}$ and $A_{{\rm Ca}^{2+},{\rm open}}$ are
specified chemical potential amplitudes to make the membrane mostly
impermeable or semi-permeable to calcium ions (Table 2). For the domain of the
membrane without ion channels, the chemical potential for calcium is fixed.
The regulation for calcium-activated potassium channels is mostly the same,
but the regulated chemical potential is specified to potassium ions.

###  2.5 Numerical implementation

The Poisson equation, Eq. (3) is solved by Fourier transformation with
electrical density given by the computed concentration for each ionic species.
The electrodiffusion equations (Eqs. 4 and 5) are solved for the concentration
of each ionic species by the backward Euler method with a second-order Godunov
upwind method (Lee et al. 2010). We employ the Krylov subspace iterations
(GMRES) using the PETSc library (Balay et al. 1997, Balay et al. 2019a,b)
preconditioned by the algebraic multigrid methods in the _hypre_ library
(HYPRE, Falgout and Yang 2002). The numerical algorithm for the ion channel
gating is based on the Monte Carlo method, determining whether to transit, and
secondly the state to transit if needed. For the details, see Lee et al. 2013.

##  3 RESULTS AND DISCUSSION

For the two-dimensional simulation, a periodic square domain with dimensions 2
$\mu$m $\times$ 2 $\mu$m is covered by a Cartesian grid containing 256
$\times$ 256 points with the uniform grid size $h$ in each direction. The
model dendritic spine has a diameter 1 $\mu$m and is centered within the
computational domain. A timestep $\Delta t$ = 30 $\mu$s is used for all
computations, including both the electrodiffusion and the Markov process that
controls the opening and closing of membrane channels. The electrical
potential $\varphi$, the chemical potential $\psi_{i}$ and the concentration
$c_{i}$ for $i^{\rm th}$ ionic species are defined on the cell-centered grids.
For the support of the chemical potential kernel, $w=6h$ is prescribed.

The membrane changes its permeability to Ca2+ and K+ in a voltage- and
calcium-dependent manner according to the continuous-time Markov process
described above with spatially localized ion channels. For the simulation,
five voltage-sensitive calcium channels and five calcium-activated potassium
channels are uniformly distributed on the head’s upper half-circle,
representing a postsynaptic density (PSD) with 15 nm distance between two
types of channels as a complex. These complexes are labeled indices from 1 to
5 counter-clockwise. Four ion species and background charges are considered
with initial concentrations different between extracellular and intracellular
domains, as shown in Table (1).

The membrane is stochastically permeable to calcium by the voltage-sensitive
calcium channel. The amount of influx through the channel is regulated by
controlling the chemical potential with two parameters, $A_{{\rm Ca}^{2+},{\rm
open}}$ and $A_{{\rm Ca}^{2+},{\rm closed}}$ as prescribed in Table (2).
Similarly, the membrane is stochastically permeable to potassium by the
calcium-activated potassium channel. The amount of influx through the channel
is regulated by controlling the chemical potential with two parameters,
$A_{{\rm K}^{+},{\rm open}}$ and $A_{{\rm K}^{+},{\rm closed}}$.

###  3.1 Stochastic calcium influxes by a train of back-propagating action
potential waves: non-local feedback loop between two-types of channels

The depolarization of the membrane is implemented by lowering the chemical
potential of sodium to 42.74 K${}_{\rm B}T$ at $t$=3 ms. The membrane is
repolarized by raising the sodium chemical potential to the initial level,
53.43 K${}_{\rm B}T$ at $t$=21 ms. With the membrane voltage in the
depolarization range, there comes calcium ionic influx from the extracellular
domain. The locally regulated chemical potential of Ca2+ around the area of
ion channels is shown in Fig. 1 with its dynamic changes in the amplitude.

The calcium concentration distributions in the y-section are shown when
calcium ions are flowed in, and influenced by the electrodiffusion at $t$=4.8
ms and at $t$=5.4 ms (Fig. 4a). The calcium concentration distribution on the
2D domain is also shown zoomed in at $t$=4.8 ms in Fig. 4(b). The electrical
potential distributions in the y-section and the 2D domain are shown at
$t$=7.8 ms and $t$=8.4 ms (Fig. 5). At $t$=7.8 ms, two voltage-sensitive
calcium channels are open, and the associated chemical potential is locally
lowered. Accordingly, the membrane is more depolarized locally, showing two
spikes. At $t$=8.4 ms, a calcium-activated potassium channel is open, and the
corresponding chemical potential is locally decreased. There follows a
considerable upstroke of electrical potential adjacent to the open BK channel,
and the expected membrane repolarization does appear far away _non-locally_ ,
that is an unexpected result.

The ionic current and membrane voltage in the time course of the channel
gating are presented in Fig. 6. When the membrane is depolarized, and the
channels are open, we observe the calcium currents. Interestingly even when
the membrane is repolarized, it takes a while for the calcium channels are all
closed. The stochastic calcium-activated potassium channel 1 is shortly open
around $t$=4 ms before the intracellular calcium concentration is elevated
enough. The evident response of the calcium-activated potassium channel after
the calcium inflow by the voltage-sensitive calcium channels is shown for the
complexes 3 and 5 with some persistent outflow of the potassium ionic current.
Here one interesting aspect of membrane voltage change is that local membrane
voltages are increased around the complexes 3 and 5 (red circled), and
decreased around the complexes 1 and 2 (blue circled). This shows that the
ionic current from the BK ion channel gating is influential _non-locally_ with
the fast redistribution of the electrical potential by the electrodiffusion,
as shown in Fig. 5(c) and 5(d).

###  3.2 Differential calcium handling by channel distributions

We have tested whether different distributions of the multiple channels
generate different calcium influxes, and whether the inhibition of the BK
channel drives more calcium influx. We compare the averaged intracellular
calcium concentration at $t$=30 ms for the non-uniformly distributed case with
channel complexes uniformly distributed in the whole dendritic spine head
(Fig. 7, first two datasets). The two-sample $t$-test shows a significantly
higher intracellular calcium level for the non-uniformly distributed case. It
is supposed that calcium influx induces membrane depolarization, and this
provides positive feedback to adjacent voltage-sensitive calcium channels for
them to have a higher probability of opening channels. This interaction is
thought to be stronger than repolarizing effects from adjacent calcium-
activated potassium channels. When the channel is open, the intracellular
calcium concentration around a voltage-sensitive calcium channel is elevated,
and the adjacent calcium-activated potassium channel gradually responds by
opening it (Figs. 5c and 5d). What happens here is local depolarization and
non-local repolarization. After that, this might enhance the coupled voltage-
sensitive calcium channel to stay in an open state. The only action of channel
inhibition is by calcium-sensitive inactivation in the channel itself.

We also compare the intracellular calcium concentrations with two cases of the
BK channels expressed with their non-uniform distribution or knocked out (Fig.
7, second and third datasets). The paired-sample $t$-test shows a significant
elevation of calcium influx by the BK channel inhibition or the knock-out. The
significant difference in the intracellular calcium concentration shows the
inhibitory mechanism on calcium influx through interaction between voltage-
sensitive calcium channels and calcium-activated potassium channels.

##  4 Conclusions

By applying the immersed chemical potentials on the two-dimensional domain, we
realize the regulation of membrane permeability to each ionic species and,
consequently, the ion selectivity of the membrane. With the continuous-time
Markov process in the stochastic feature of the voltage-sensitive ion channel,
the ion channel gating and the corresponding ionic current are observed
depending on the membrane voltage. With discretely placed ion channels,
individual ion channel gating is observed with spatiotemporally changing
chemical potential distributions. Compared to the point cell modeling, this
kind of approach makes it available for us to consider the spatially _non-
uniform_ distribution of ion channels and their spatial effects on the
physiology. Accordingly, we can implement the change in the ion channel
density in the sense of synaptic plasticity.

The _non-local_ interaction between voltage-sensitive calcium channels and
calcium-activated potassium channels is thought to be captured by treating the
electrophysiology with electrodiffusion in nano- to micro-scales, and
undoubtedly not feasible by drift-diffusion. At the same time, positive
feedback among _local_ channel complexes for calcium influx is also predicted
for further validation by experiments. The back-propagating action potential
was treated by the one-time event of phasic trains. In real neuronal
activities, the frequency of this kind of back-propagating action potential
spikes may induce differential responses from calcium-activated potassium
channels (Kuznetsov et al. 2006, Hage and Khaliq 2015). We expect to
reconstruct the three-dimensional electrodiffusion (Mori 2007, Cartailler and
Holcman 2019), which is straightforward due to the Cartesian grids where the
chemical potentials are immersed and the electrodiffusion is computed.

## Acknowledgments

The author appreciates careful comments from Charles Peskin.

Figure 1: Chemical potential distribution in position and in the course of ion
channel gating on 2D: it is placed along the boundary. The support of the
chemical potential kernel is a square that measures 24 mesh widths $\times$ 24
mesh widths. The chemical potential is locally and stochastically regulated at
the domain of ion channels. (a) The chemical potential for calcium is high,
i.e., five ion channels are all closed. (b) and (c) The chemical potentials
for calcium are low around 1 and 3 ion channels.

Figure 2: Markov chain of the voltage-sensitive calcium channel: The calcium
channel has four subunits. In the state of $C_{i}$, $i-1$ subunits are open.
In the state of $\rm I$, the channel is in inactivation and closed. The ion
channel is open only when it stays on the state of $\rm O$. The rate constants
for the opening and closing of each subunit are denoted by $\alpha$ and
$\beta$. Courtesy of Cox cox

.

Figure 3: Markov chain of the calcium-activated potassium channel: The channel
has four subunits, and the calcium-binding rate to each subunit is calcium-
dependent. The rate constants for the opening and closing in the five
different calcium bound states are denoted by $k_{i}$ and $k_{-i}$, each
dependent on the membrane voltage. Courtesy of Cox cox

.

$\begin{array}[]{c@{\hspace{0.1in}}c}\includegraphics[angle={0},width=433.62pt]{Figure4}\hfil\hskip
7.22743pt\end{array}$

Figure 4: Concentration distribution with ion channels open on 2D; (a) The
calcium concentration distribution in y-section at $t$=0 ms (red), $t$=4.8 ms
(black-dotted), $t$=5.4 ms (black-solid). The chemical potential distribution
is also shown (blue) at $t$=4.8 ms. (b) The calcium concentration is zoomed in
to resolve the low-level concentration in $\mu$M in the intracellular domain
at $t$=4.8 ms.

$\begin{array}[]{c@{\hspace{0.1in}}c}\includegraphics[angle={0},width=433.62pt]{Figure5}\hfil\hskip
7.22743pt\end{array}$

Figure 5: Electrical potential distribution on 2D; (a) The electrical
potential in the y-section at $t$=0.2 ms (red) and $t$=7.8 ms (black). The
chemical potential for the calcium ions is also shown in y-section at $t$=7.8
ms (blue). (b) The electrical potential on the 2D domain at $t$=7.8 ms. (c)
The electrical potential in the y-section at $t$=0.2 ms (red) and $t$=8.4 ms
(black). The chemical potential for potassium is also shown in y-section at
$t$=8.4 ms (blue). (d) The electrical potential on the 2D domain at $t$=8.4
ms.

$\begin{array}[]{c@{\hspace{0.1in}}c}\includegraphics[angle={0},width=433.62pt]{Figure6}\hfil\hskip
7.22743pt\end{array}$

Figure 6: Voltage-sensitive calcium channel and calcium-activated potassium
channel activities; The columns represent (a) calcium ionic currents by the
voltage-sensitive calcium channels, (b) potassium ionic currents by the
calcium-activated potassium channels, (c) membrane voltage, and (d)
intracellular calcium concentration adjacent to each voltage-sensitive calcium
channel. The rows represent five individual complexes of the voltage-sensitive
calcium channel and the calcium-activated potassium channel. These
electrophysiological activities are from the non-uniform distribution of the
channel complexes.

$\begin{array}[]{c@{\hspace{0.1in}}c}\includegraphics[angle={0},width=216.81pt]{channel_distribution_BK_influence}\hfil\hskip
7.22743pt\end{array}$

Figure 7: Calcium ion concentration comparison; ten samples are collected for
the comparison of the averaged intracellular calcium ion concentrations for
three cases, 1) the ion channel complexes are uniformly distributed, 2) the
ion channel complexes are non-unformed distributed, 3) the ion channels
complexes are non-uniformly distributed, and the BK channels are knocked out.
The activity of the BK channel functions as inhibiting the calcium inflow.
Non-uniform distribution of the ion channel complexes induces higher
intracellular calcium concentrations. * p $<$ 0.05 two-sample $t$-test, ** p
$<$ 0.05 paired-sample $t$-test

Table 1: Initial concentrations in all simulations (mM). X- denotes the fixed background charge and the concentrations stated for X- refer to the concentration of the background charges, not the concentration of the molecules that carry the background charges.  Ion species | extracellular concentration | intracellular concentration
---|---|---
Ca2+ | 2.0 | 0.0002
Cl- | 150 | 13
Na+ | 150 | 15
K+ | 5 | 100
X- | 9 | 102.0004
Table 2: Chemical potential (KBT) prescribed for each ionic species.  Ion species | chemical potential high | chemical potential low
---|---|---
Ca2+ | 53.43 | 8.014
Cl- | 45.41 | 45.41
Na+ | 53.43 | 42.74
K+ | 53.43 | 18.70
X- | 53.43 | 53.43

## References

  * (1) S. Balay, V. Eijkhout, W.D. Gropp, L.C. McInnes, B.F. Smith, Efficient management of parallelism in object-oriented numerical software libraries, in: E. Arge, A.M. Bruaset, H.P. Langtangen (Eds.), Modern Software Tools in Scientific Computing, Birkhauser Press, pp. 163-202 (1997)
  * (2) S. Balay, S. Abhyankar, M.-F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W.D. Gropp, D. Karpeyev, D. Kaushik, M.G. Knepley, D.-A. May, L.C. McInnes, B.F. Smith, S. Zampini, H. Zhang, PETSc Users Manual, Technical Report, ANL-95/11 Revision 3.12, Argonne National Laboratory (2019)
  * (3) S. Balay, S. Abhyankar, M.-F. Adams, J. Brown, P. Bruen, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W.D. Gropp, D. Karpeyev, D. Kaushik, M.G. Knepley, L.C. McInnes, R.T. Mills, T. Munson, K. Rupp, P. Sanan, B.F. Smith, S. Zampini, H. Zhang, PETSc, http://www.mcs.anl.gov/petsc (2019)
  * (4) H. Berkefeld, C.A. Sailer, W. Bildl, V. Rohde, J.-O. Thumfart, S. Eble, N. Klugbauer, E. Reisinger, J. Bischofberger, D. Oliver, H.-G. Knaus, U. Schulte, B. Fakler, BKCa-CaV channel complexes mediate rapid and localized Ca2+-activated K+ signaling, Science, 314, 615-620 (2006)
  * (5) B.L. Bloodgood and B.L. Sabatini, Nonlinear Regulation of Unitary Synaptic Signals by CaV2.3 Voltage-Sensitive Calcium Channels Located in Dendritic Spines, Neuron, 53, 249-260 (2007)
  * (6) B.L. Bloodgood and B.L. Sabatini, Ca2+ signalling in dendritic spines, Curr. Opin. Neurobiol., 17, 345-351 (2007)
  * (7) T. Bonhoeffer and R. Yuste, Spine motility: phenomenology, mechanisms, and function, Neuron, 35, 1019-1027 (2002)
  * (8) M. Bosch and Y. Hayashi, Structural plasticity of dendritic spines, Curr. Opin. Neurobiol., 22, 1-6 (2011)
  * (9) J. Cartailler and D. Holcman, Steady-state voltage distribution in three-dimensional cups-shaped funnels modeled by PNP, J. Math, Biol., 79, 155-185 (2019)
  * (10) D.H. Cox, Modeling a Ca2+ channel/BKCa channel complex at the single-complex level, Biophy, J., 107, 2797-2814 (2014)
  * (11) R. D. Falgout, U. M. Yang, hypre: a library of high performance preconditioners, in: P. M. A. Sloot, C. J. K. Tan, J. J. Dongarra, A. G. Hoekstra (Eds.), Computational Science - ICCS 2002 Part III, Vol. 2331 of Lecture Notes in Computer Science, Springer-Verlag, 2002, pp. 632–641, also available as LLNL Technical Report UCRL-JC-146175.
  * (12) T. Griffith and K. Tsaneva-Atanasova and J.R. Mellor, Control of Ca2+ influx and calmodulin activation by SK-channels in dendritic spines, PLoS Comput. Biol., 12, e1004949 (2016)
  * (13) J.J. Gutzmann and L. Lin and D.A. Hoffman, Functional coupling of Cav2.3 and BK potassium channels regulates action potential repolarization and short-term plasticity in the mouse hippocampus, Front. Cell. Neurosci., 34, 5261-5272 (2019)
  * (14) T.A. Hage and Z.M. Khaliq, Tonic firing rate controls dendritic Ca2+ signaling ad synaptic gain in substantia nigra dopamine neurons, J. Neurosci., 35, 5823-5836 (2015)
  * (15) C. Piochon, M. Kano, and C. Hansel, LTD-like molecular pathways in developmental synaptic pruning, Nat. Neurosci., 19, 1299-1310 (2016)
  * (16) K.M. Harris and S.B. Kater, Dendritic spines: cellular specializations imparting both stability and flexibility to synaptic function, Annu. Rev. Neurosci., 17, 341-371 (1993)
  * (17) S. He, Y.-X. Wang, R.S. Petralia, S.D. Brenowitz, Cholinergic modulation of large-conductance calcium-activated potassium channels regulates synaptic strength and spine calcium in cartwheel cells of the dorsal cochlear nucleus, J. Neurosci., 34, 5261-5272, (2014)
  * (18) D. Holcman and R. Yuste, The new nanophysiology: regulation of ionic flow in neuronal subcompartments, Nat. Rev. Neurosci., 16, 685-692 (2015)
  * (19) HYPRE: High performance preconditioners, http://www.llnl.gov/CASC/hypre.
  * (20) I.H. Kim, M.A. Rossi, D.K. Aryal, B. Racz, N. Kim, A. Uezu, F. Wang, W.C. Wetsel, R. Weinberg, H. Yin, and S.H. Soderling, Spine pruning drives antipsychotic-sensitive locomotion via circuit control of striatal dopamine, Nat. Neurosci., 18, 883-891 (2015)
  * (21) A.J. Koleske, Molecular mechanisms of dendrite stability, Nat. Rev. Neurosci., 14, 536-550 (2013)
  * (22) A.S. Kuznetsov, N.J. Kopell, and C.J. Wilson, Transient high-frequency firing in a coupled-oscillator model of the mesencephalic dopaminergic neuron, J. Neurophysiol., 95, 932-947 (2006)
  * (23) R. Lamprecht, C.R. Farb, S.M. Rodrigues, and J.E. LeDoux, Fear conditioning drives prolifin into amygdala dendritic spines, Nat. Neurosci., 9, 481-483 (2006)
  * (24) P. Lee, The immersed boundary method with advection-electrodiffusion, Ph.D. thesis, Courant Institute of Mathematical Sciences, New York University (2007)
  * (25) P. Lee and B.E. Griffith and C.S. Peskin, The immersed boundary method for advection-electrodiffusion with implicit timestepping and local mesh refinement, J. Comp. Phys., 229, 5208-5227 (2010)
  * (26) P. Lee, E.A. Sobie, C.S. Peskin, Computer simulation of voltage sensitive calcium ion channels in a dendritic spine, J. Theo. Biol., 338, 87-93 (2013)
  * (27) R. Lujan, C. Aguado, F. Ciruela, X.M. Arus, A. Martin-Belmonte, R. Alfaro-Ruiz, J. Martinez-Gomez, L. de la Ossa, M. Watanabe, J.P. Adelman, R. Shigemoto, Y. Fukazawa, SK2 channels associate with mGlu1α receptors and CaV2.1 channnels in purkinje cells, Front. Cell. Neurosci., 12, 1-16 (2018)
  * (28) I. Brunig and S. Kaech and H. Brinkhaus and T. G. Oertner and A. Matus, Influx of extracellular calcium regulates actin-dependent morphological plasticity in dendritic spines, Neuropharmacology, 47, 669-676 (2004)
  * (29) Y. Mori and J. W. Jerome and C.S. Peskin, A three-dimensional model of cellular electrical activity, Bull. Inst. Math. Acad., 2, 367-390 (2007)
  * (30) C.S. Peskin, The Immersed Boundary Method, Act. Num., 11, 479-517 (2002)
  * (31) C.S. Peskin, Mathematical Aspects of Neurophysiology, Lecture Note, Courant Institute of Mathematical Sciences, New York University (2000)
  * (32) R.A. Rice, E.E. Spangenberg, H. Yamate-Morgan, R.J. Lee, R.P.S. Arora, M.X. Hernandez, A.J. Tenner, B.L. West, K.N. Green, Elimination of microglia improves functional outcomes following extensive neuronal loss in the hippocampus, J. Neurosci., 35, 9977-9989 (2015)
  * (33) P. Sah and E.S.L. Faber, Channels underlying neuronal calcium-activated potassium currents, Progress in Neurobiology, 66, 345-353 (2002)
  * (34) M. Sheng and M. Kim, Postsynaptic Signaling and Plasticity Mechanisms, Science, 298, 776-780 (2002)
  * (35) G.M. Shepherd, The Dendritic Spine: A Multifunctional Integrative Unit, J. Neurophysiol., 75, 2197-2210 (1996)
  * (36) B.L. Sabatini and K. Svoboda, Analysis of calcium channels in single spines using optical fluctuation analysis, Nature, 408, 589-593 (2000)
  * (37) R. Yasuda and B.L. Sabatini and K. Svoboda, Plasticity of calcium channels in dendritic spines, Nat. Neurosci., 6, 948-955 (2003)
  * (38) T. Tada and M. Seung, Molecular mechanisms of dendritic spine morphogenesis, Curr. Opin. Neurobiol.,16, 95-101 (2006)


# Symmetry-group-protected microfluidics for multiplexed stress-free
manipulations

Jeremias Gonzalez  Department of Physics, University of California, Merced,
Merced, CA 95343, USA  Ajay Gopinathan  Department of Physics, University of
California, Merced, Merced, CA 95343, USA  Bin Liu  Department of Physics,
University of California, Merced, Merced, CA 95343, USA

###### Abstract

Modern micromanipulation techniques typically involve trapping using
electromagnetic, acoustic or flow fields that produce stresses on the trapped
particles thereby precluding stress-free manipulations. Here, we show that by
employing polyhedral symmetries in a multichannel microfluidic design, we can
separate the tasks of displacing and trapping a particle into two distinct
sets of flow operations, each characterized and protected by their unique
groups of symmetries. By combining only the displacing uniform flow modes to
entrain and move targeted particles in arbitrary directions, we were able to
realize symmetry-protected, stress-free micromanipulation in 3D. Furthermore,
we engineered complex, microscale paths by programming and controlling the
flow within each channel in real-time, resulting in multiple particles
simultaneously following desired paths in the absence of any supervision or
feedback. Our work therefore provides a general symmetry-group-based framework
for understanding and engineering microfluidics and a novel platform for 3D
stress-free manipulations.

##  I Introduction

Microfluidics enables the study of the behavior, control, and manipulation of
microscale flows typically using device architectures comprising channels and
junctions as their building blocks [1, 2, 3, 4]. Desired flow structures and
microfluidic functions such as transporting and mixing chemicals [5], shaping
flow profiles [6, 7], encapsulating multiphase fluids [8, 9], and trapping
microparticles [10, 11] are incorporated through the geometry of these
building blocks and the topology of their assemblies [12, 13, 14, 15]. Despite
this increased complexity, our understanding of these microfluidic functions
is mostly achieved in a case-by-case fashion relying on actual experimental
measurements or on computational fluid dynamics (CFD) studies using detailed
experimental geometries. A unified framework that allows us to understand how
these complex flow structures emerge from simpler, more fundamental flows
would prove extremely useful for the fast, modular development of rationally
designed flow structures and associated microfluidic functions.

Symmetry is a concept that has been used to understand phenomena and
structures in disciplines ranging from fundamental physics to viral structure
[16, 17, 18]. In the realm of microscale flows, symmetries are essential
features of representative flow structures, such as the helical symmetry in
swimming microorganisms [19, 20] and the reflection symmetry applicable in
many hydrodynamic interactions [21, 22]. These symmetries have profound
implications with helical symmetry giving rise to swimming motility in
bacterial species due to the breaking of kinematic reversibility while the
reflection symmetries possessed by a pair of microswimmers constrain their
synchronization and thus their collective motion. Such considerations applied
to microfluidic systems could allow us to uncover very general results
independent of system details. Recent work on multi-channel microfluidic
junctions and open-space microfluidics have shown that their main flow
structures can be obtained from the potential flows around charges
representing the channels’ ports, without consideration of detailed channel
geometries [11, 23]. This simplification has revealed intriguing flow
structures that enable rich microfluidic applications, such as flow-driven
traps for particle manipulation and dynamic confinements of fluid into
multipole-like networks. Additionally, certain combinations of rotation and
reflection symmetries have been applied both analytically and numerically to a
multi-port microfluidic junction to eliminate strain components at the
junction center, giving rise to potentially perturbation-free
micromanipulations [24]. All the above implications of symmetries suggest the
feasibility of using a more generalized symmetry-based approach, i.e., the
symmetry group, for a more fundamental understanding of the origins of these
symmetries and ways to manipulate them for specific functions.

One of the most desirable of these microfluidic functions is manipulating flow
structures that allow the displacement of individual microscale objects along
arbitrary paths [25]. Typically, micromanipulation techniques involve trapping
particles using electromagnetic, acoustic, or flow fields that produce
restoring forces in the vicinity of the trap [26, 27, 28, 29]. For example, in
cases where hydrodynamic flow alone is used for trapping, previous studies
have shown that a Stokes trap can be achieved by managing the locations of
hyperbolic points through multiple channels that intersect at a middle
junction [11]. The location of such a trap can be dynamically adjusted to
realize the direct manipulation of trapped particles, which has also been
extended to 3D flows very recently [30].

Computation-assisted feedback-control has also enabled manipulations of
particles along local flow fields away from a hyperbolic point [31, 25]. The
displacements of particles can be adjusted by aligning the local flow field
with the moving directions of particles. However, in any such
micromanipulation technique, trapping and displacement are inextricably
coupled. As illustrated in Fig. 1(a), flow fields far from these manipulated
particles are divergent, leading to inevitable stress perturbations on any
manipulated objects. It is worth noting that the stress distribution on the
particle is intrinsic to the far-field hyperbolic topology [32]. Even for a
flow field that appears uniform away from the hyperbolic point (Fig. 1(a)
inset), the stress distribution on the entrained particle is identical to that
within the trap, regardless of their relative displacement. These stressful
perturbations thus restrict such approaches to studying stress-insensitive
phenomena. Scaling such approaches with traps or local flow vector fields to
simultaneously manipulate large numbers of particles is also a challenge.
Here, we explore the use of symmetry groups to guide the design of
microfluidics that enable multiplexed and stress-free manipulation, which is
characterized by a truly uniform flow field around the manipulated particle(s)
(Fig. 1(b)). Specifically, we identify and characterize microfluidic flow
structures based on symmetry groups, which allow us to distinguish and
robustly realize fundamentally different microfluidic functions.

##  II Results

###  II.1 2D Microfluidic symmetries

We start with a cross-channel junction for its well-understood microfluidic
functions and relatively simple symmetries. Here, the four ports of identical
channels (with flow rates $f_{i}$, $i=1,2,3,4$) form a square (Fig. 2(a)), and
therefore its flow structures in the plane are restricted by the symmetry of a
square. An arbitrary flow pattern generated by a combination of $\\{f_{i}\\}$
at this junction corresponds to a set of identical flow manipulations through
certain rotations or reflections (Supplemental Material, Fig. S1). All such
possible rotations or reflections (including the trivial identity operation)
form a dihedral-4 (D4) symmetry group [16], which restricts the possible flow
patterns. This D4 group thus represents the flow manipulation functions in a
cross-channel junction.

How do the D4 group elements map to the actual flow functions? Analogous to
the subgroups contained in one symmetry group, these flow manipulation
functions are reducible into smaller classes with unique features. We note
first that there are only three independent flow modes, resulting from the
number of controls (four channels) minus the number of constraints (one
associated with the conservation of volume, i.e., $\sum_{i=1}^{4}f_{i}=0$).
Through inspection, one can readily identify these three modes, forming two
classes of flow functions. First, fluxes with opposite signs through a pair of
opposing channels give rise to displacement flows along two orthogonal
directions ($x$ and $y$), corresponding to two fundamental modes (Fig. 2(b)).
Second, fluxes with same sign through both pairs of opposing channels provide
one more fundamental elongation flow mode (Fig. 2(c)). Any flow profile
produced in this junction (with an example shown in Fig. 2(a)) must be a
combination of these three flow modes. However, none of these fundamental
modes can be generated by combining any of the rest, a consequence of their
linear independence. Importantly, the two displacement modes have a reflection
symmetry (along $45^{\circ}$ from the $x$ axis), forming a reflection group D1
or a cyclic group C2, which is known to be a subgroup of the D4 group [33].
The single elongation flow mode also represents a trivial subgroup of D4 that
contains only the identity operation.

In addition to forming different subgroups, these two classes of flows also
have distinct characteristic functions, here displacement vs. elongation. To
abstract these functions from the detailed flow geometries, we map the cross-
channel junction to a square, where the flux at each channel is mapped to an
effective charge at each lattice site ($q_{1}$, $q_{2}$, $q_{3}$, $q_{4}$),
representing a point-like fluid source (positive) or sink (negative). The flow
field is thus well-defined as a potential flow, with the flow velocity
$\mathbf{u}$ given by the gradient of the potential $\phi$, i.e.,
$\mathbf{u}=\bm{\nabla}\phi$. Since the classes of flow functions are
independent of the choice of coordinates, the candidate quantities for such
purposes are the strain rate invariants [34] that do not change under any
rotation of the microfluidic device. The first-order invariant $I_{1}$, which
is the trace of the rate of strain tensor $\dot{\gamma}$, is trivial and
vanishes for incompressible fluids. We therefore consider the first non-
trivial invariant $I_{2}=-\frac{1}{2}\mathrm{Tr}(\dot{\gamma}^{2})$. For a 2D
potential flow with a total number of $N=4$ sources and sinks (with magnitudes
{$q_{i}$}, located at {$\mathbf{r}_{i}$}), the scalar potential can be
expressed as
$\phi(\mathbf{r})=\sum_{i=1}^{4}q_{i}\ln\left(\left|\mathbf{r}-\mathbf{r_{i}}\right|\right)$.
Since the rate of strain is linearly related to the magnitudes of sources (or
sinks), we can show (Supplemental Material, Equations S8-9) that the invariant
$I_{2}$ at the center of the junction has a quadratic form with respect to
$q_{i}$, i.e., $I_{2}=\sum_{i,j=1}^{4}L_{ij}q_{i}q_{j}$, where $L_{ij}$ is the
element of a coupling tensor, depending on the distance between the pair of
$q_{i}$ and $q_{j}$. On a square, there are only two possible types of pairs:
pairs along the lateral sides or along the diagonals. With these
considerations, we can show (Supplemental Material, Equation S15) that the
$\mathbf{L}$ tensor has the form

$\mathbf{L}=\left[\begin{array}[]{cccc}0&a&b&a\\\ a&0&a&b\\\ b&a&0&a\\\
a&b&a&0\end{array}\right].$ (1)

Elements “a” and “b” are symbolic representations of the above lateral and
diagonal couplings, respectively, forming a symmetric circulant matrix for
$\mathbf{L}$. A circulant matrix has all the same rows with each row being
rotated one element to the right relative to the row above. One mutual
eigenvector for any circulant matrix corresponds to the trivial mode where all
four charges are identical, i.e., $(q_{1},q_{2},q_{3},q_{4})$ = $(1,1,1,1)$.
This mode cannot be accessed by a fluid system due to the invalidation of the
continuity condition. This leaves $N-1=3$ nontrivial eigenmodes (all
satisfying $\sum_{i=1}^{4}q_{i}=0$), in agreement with the degrees-of-freedom
argument. Every nontrivial mode exhibits a one-to-one map to the fluidic
manipulation function as previously classified through subgroups of the D4
group (Fig. 2(d)). More specifically, the degeneracy of two orthogonal dipoles
corresponds to a common displacement function of the flow along orthogonal
directions (Fig. 2(e)). The remaining quadrupole-like mode is nondegenerate
and maps to the elongation function (Fig. 2(f)).

###  II.2 3D Microfluidic symmetries

As demonstrated above, the eigenvalue analysis of the stress invariant
($I_{2}$) provides an automatic strategy of classifying microfluidic
functions, which we now extend to 3D structures to incorporate more
sophisticated symmetries. The potential flow now adopts its 3D form with $N$
sources,
$\phi(\mathbf{r})=\sum_{i=1}^{N}-\frac{q_{i}}{\left|\mathbf{r}-\mathbf{r}_{i}\right|}$.

For relatively more practical applications in microfluidics, we consider here
the 3D symmetries contained in structures with fewer vertices, namely the
tetrahedron ($N=4$) and the octahedron ($N=6$).

A tetrahedron has a permutation S4 symmetry [33], corresponding to an
unaltered structure (identical through rotation and reflection) by permuting
all four vertices (Fig. 3(a)). The eigenvalue analysis shows that all three
eigenmodes for this tetrahedron are degenerate (Fig. 3(b)), due to equivalent
neighbors for every vertex. These three modes together represent a threefold
rotation symmetry group (C3) (achieved by $120^{\circ}$ rotation about any of
the face norms), which is a subgroup of S4. Notably, the potential flow of
such a mode exhibits a mixture of dipole and quadrupole moments (upper Fig.
3(b)), leading to a microfluidic function corresponding to simultaneously
displacing and elongating (transverse to the displacement axis) the fluid
within the junction. To obtain a purely displacing flow in 3D, we therefore
look to the next available polyhedral symmetry, i.e., the octahedral symmetry.

An octahedron (Fig. 3(c)) can be formed by the middle points of the six edges
of a tetrahedron (Fig. 3(a)). Its symmetry group is isomorphic to
S4$\times$C2, with the additional symmetry arising from an extra 2-fold
rotation symmetry that is absent for a tetrahedron (Fig. 3(a)) [33]. A similar
eigenvalue analysis leads to two groups of degenerate modes, corresponding
respectively to dipoles and quadrupoles. The three degenerate dipole modes
($p_{x}$, $p_{y}$, $p_{z}$) are orthogonal in direction, corresponding to a C3
group (through $120^{\circ}$ rotation about any of the face norms), one
subgroup of the octahedral symmetry group [33]. In such dipole modes, a pair
of effective charges along the diagonal (e.g., $q_{1}$ and $q_{2}$ for
$p_{x}$) have opposite signs (with the rest of the charges being neutral),
corresponding to the activation of a pair of a source and a sink. The two
degenerate quadrupoles are distributed in two perpendicular planes (here,
$x$-$y$ and $x$-$z$ planes), forming a D1 subgroup (switchable by a reflection
along the middle plane between $x$-$y$ and $x$-$z$ plane). In such quadrupole
modes, a pair of sources and a pair of sinks along two diagonals are
activated. It should be noted that an equivalent quadrupole in the $y$-$z$
plane can be achieved by superimposing the above two quadrupoles, and is thus
not one of the eigenmodes. Analogous to the square symmetry in the 2D case,
the octahedral symmetry gives rise to completely separate displacement and
elongation flow functions (Fig. 3(d) and (e)), classified by distinct
symmetries. The pure-dipole-like potential in the displacements also ensures
that the rate-of-strain invariant $I_{2}$ is zero in the middle of the
junction. Thus, a superposition of these three displacement flow modes
generates 3D omnidirectional flows while preserving the stress-free condition
(at least in the middle of the junction).

###  II.3 Realizing 3D Stress-free microfluidic manipulations

To realize such a stress-free microfluidic junction while enabling microscope
observation, we rotate the octahedron so that one of its faces (e.g.,
$(1,1,1)$ axis) is aligned along the visualization axis ($z^{\prime}$). This
rotated geometry fits well into a double-layer microfluidic design [24]: six
channels ($0.5$ mm wide and $0.3$ mm deep in cross-sections) at two elevations
intersecting a middle cylindrical chamber (with radius $R=1.0$ mm and height
$H=1.6$ mm) along the radial directions, with the locations of intersections
matching that of the octahedral vertices (Fig. 4(a)). This rotation is also
determined by the coordinates of the laboratory frame of reference
$(x^{\prime},y^{\prime},z^{\prime})$, with a radial channel (here, $f_{1}$)
aligned along the $x^{\prime}$ axis. We fabricated the microfluidic channel
part of this device by fusing multiple sheets of laser-etched glass
(Citrogene, see Methods), which was then mounted on a customized adapter for
flow control and microscopy.

To ensure that a displacement flow can be robustly created in all possible
directions, we first considered flows along the three orthogonal axes
($x^{\prime}$, $y^{\prime}$, $z^{\prime}$) of the laboratory coordinates
forming the basis of the velocity space. The potential representations
(isosurfaces) of these orthogonal flows (along $x^{\prime}$, $y^{\prime}$, and
$z^{\prime}$) are shown in Fig. 4(b), which are linear combinations of the
above three degenerate dipole modes with their coefficients
$(c_{1},c_{2},c_{3})$ being $\frac{1}{\sqrt{6}}(2,-1,-1)$,
$\frac{1}{\sqrt{2}}(0,1,-1)$, and $\frac{1}{\sqrt{3}}(1,1,1)$, respectively.
In experiments, these effective charges are realized by independently
offsetting the fluid pressure on the corresponding channel using two 4-channel
piezoelectric regulators (Elveflow OB1 MK3+), with a positive offset for a
positive charge (or source) and a negative offset for a negative one (or
sink). Custom software was programmed to convert the above flow coefficients
into actual pressure settings for flow generation (Methods). The three
orthogonal flows ($x^{\prime}$, $y^{\prime}$, $z^{\prime}$) were successfully
generated in our microfluidic device, within the same volume near the center
(Fig. 4(c)-(e)), as demonstrated by the 3D traces of the seeding particles
(Methods). Despite their different geometries in flow dipoles (Fig. 4(b)),
these flows are all uniform over almost the entire 3D volume ($\approx$ $200$
$\mu$m $\times$ $200$ $\mu$m $\times$ $150$ $\mu$m) captured by the
microscope, demonstrating the experimental realization of a robust
“stressless” condition that is guaranteed by symmetry.

To further investigate the experimental capacity of such stress-free
microfluidics, we incorporated temporal dependence in the flow by dynamically
varying the pressures on all channels. To facilitate this time-dependent
control, we represented each possible configuration of the displacement flow
by its coefficients $(c_{1},c_{2},c_{3})$ when expressed in terms of three
degenerate dipole modes, which essentially form a 3D phase space (Fig. 5(a)).
Any time-dependent flow manipulation can thus be generated by a series of
points in this phase space. Noting that an axial ($z^{\prime}$) flow is
represented by a vector along the $\frac{1}{\sqrt{3}}(1,1,1)$ direction in the
phase space, all orthogonal flows must satisfy $c_{1}+c_{2}+c_{3}=0$, which
restricts all horizontal (or in-plane) displacement flows to a plane in the
phase space. Hopping in this plane with equal angular separations and
distances with respect to the origin $(0,0,0)$ gives rise to polygon-shaped
flow patterns, for instance, a triangle, a square and a circle in the
continuum limit (Fig. 5(b)-(d), Movies S1-S3). Such patterns are almost
identical for all seeding particles, with the fluid in the bulk translating
like a piece of solid.

By incorporating phase spaces out of the $c_{1}+c_{2}+c_{3}=0$ plane, we also
generated time-dependent flows in 3D. For instance, imposing an oscillatory
$z^{\prime}$ motion (with doubled frequency) to the above circular mode gives
rise to a 3D flow that represents a Lissajous curve (Fig. 6(a), Movie S4).
Invariably, individual seeding particles within the volume of observation
trace out the same pattern as desired. demonstrating the robustness of the
stressless condition even under 3D dynamic control.

To assess the flow uniformity, we center the trajectory of every seeding
particle within the entire volume and show them side by side. It is clear that
these trajectories overlap with each other very well over a full period. In
perfectly stress-free flows, all trajectories must be identical, leaving no
deviation from the mean path. We thus defined a characteristic strain based on
the deviation from the mean trajectory, i.e.,
$\epsilon(t)=\textrm{deformation}/\textrm{size}=\Delta\tilde{r}(t)/\Delta{r}(t)$,
where
$\Delta{r}\equiv\sqrt{\left<(\mathbf{r}-\left<\mathbf{r}\right>)\cdot(\mathbf{r}-\left<\mathbf{r}\right>)\right>}$
and
$\Delta\tilde{r}(t)\equiv\sqrt{\left<(\tilde{\mathbf{r}}-\left<\tilde{\mathbf{r}}\right>)\cdot(\tilde{\mathbf{r}}-\left<\tilde{\mathbf{r}}\right>)\right>}$
are the fluctuations of displacements of all particles ($\mathbf{r}$) before
(Fig. 6(a)) and that ($\tilde{\mathbf{r}}$) after (Fig. 6(b)) centering, with
operators $\left<\cdot\right>$ corresponding to spatial averages. As shown in
Fig. 6(c), such characteristic strains are of the order of $\epsilon\sim
10^{-2}$ for typical 3D manipulations, illustrating the nearly stress-free
condition.

The agile flow responses to the controlling pattern in the phase space
indicate a direct map of the phase space to the real velocity space. The
previously mentioned control in phase space can thus be generalized to
arbitrary flow motions to realize more sophisticated manipulation
capabilities. To demonstrate this concept, we combine both continuous motions
and discrete hops in the phase space (or the velocities) to manipulate, in a
stress-free manner, individual particles for “printing” discrete letters
(here, “UCM”) onto the focal plane. Each letter was traced out by combining
motions (along smooth curves) and hops (at the corners of the letter) in the
phase space. The discrete gaps between adjacent letters were achieved by
abruptly offsetting the fluid in the axial direction so that the previously
focused particles were moved away from the focal plane before their
reappearance for “printing” the next letter. As shown in Fig. 7, this
“printing” task can indeed be achieved in our stress-free microfluidic
channel, which further demonstrates its versatile manipulation capabilities
(Movie S5).

##  III Discussion

Our results show a substantial role for symmetry groups in classifying
microfluidic functions. These symmetry-group-protected functions are
insensitive to detailed geometries, which we exploited to realize 3D stress-
free flows in a microfluidic device. Such flows are qualitatively different
from all available micromanipulation approaches, where physical traps of
certain forms must be present, resulting in stresses on the manipulated
objects. Both the stress-free modes and the stressful modes can be
characterized by distinct subgroups of a group of polyhedral symmetry
(possessed by the microfluidic device), suggesting independent transport and
trapping in microfluidic manipulations.

From our experimental observations, these stress-free flows are demonstrated
by extremely parallel trajectories of the manipulated particles, regardless of
the detailed curvatures of the trajectories. This enables us to achieve truly
multiplexed stress-free manipulations, a significant challenge for trap-based
approaches. It is worth noting that our trajectories inevitably deviate from
the ideal geometries due to the absence of feedback controls. Meanwhile, any
potential differences (e.g., different flow resistances) among the six
channels lead to channel-sensitive responses to applied pressures, which
ultimately modifies our time-dependent flow patterns. Remarkably, these
“defects” in manipulations do not alter the uniformity in the flows (with all
particle trajectories remaining parallel), suggesting the robustness in the
symmetry-protected flows.

Our work therefore opens up new avenues of experiments on microscopic
phenomena that occur in truly stress-free flows. Combinations of the above
subgroups easily lead to other subgroups with higher orders (e.g., the C4
subgroup of the D4 group), associated with flow characteristics more than the
strain invariants at the center. Flow patterns belonging to these subgroups
maybe coordinated in time to form a mutual impact on the manipulated flow,
e.g., creating a “mean” vortex that is absent from steady potential flows.
These advanced symmetry features and flow manipulation functions will be
explored in our future work.

## Acknowledgements

This work was supported by National Science Foundation Grant CBET-2046822,
CBET-1706511, NSF-CREST: Center for Cellular and Bio-molecular Machines (CCBM)
at UC Merced (HRD-1547848), and Department of Defense DURIP ARO No.73839-MA-
RIP.

##  IV Appendix: Methods

Microfluidic fabrication and its assembly for microscope observations. The
microfluidic channel was fabricated as a stack of glass sheets (0.13 - 1 mm
thick), with each containing a special pattern of channels. We made the design
of these sheets through a computer assisted design (CAD) software (Autodesk
Fusion 360), which were then sent to Citrogene for fabrications of these
customized sheets in borosilicate glass and their integration into one single
multilayer microfluidic chip (using optically clear adhesives). The overall
exterior dimension of the microfluidic chip is compact, here, 62 mm $\times$
22 mm $\times$ 2.7 mm in $x$, $y$ and $z$, comparable to the footprint of a
typical glass slide. The chip was also designed to be self-enclosed, leaving
only 6 small pores (1 mm in diameter) open for accessing the microfluidic
flows (see Supplemental Material). A 0.13 mm thick glass sheet at the bottom
layer sealed the microfluidic channels and enabled high optical quality for
conventional microscope observations. Before being connected to a multi-
channel pump, the glass chip was mounted to a customized adapter for better
sealing results.

Multiple-channel flow controls. Flows within the multilayer microfluidic
channel were generated by two 4-channel microfluidic flow controllers
(Elveflow OB1 MK3+) with each channel’s pressure regulated independently
between -0.9 and 1 bar relative to the atmospheric pressure. A custom program
(written in Python) was used to modulate the pressures on each channel and
record the microscope image in real time. To accommodate the finite resolution
of our pressure controllers, here 120 $\mu$bar, we used a water-glycerol
mixture to increase the viscosity of the fluid in the microfluidic channel,
which maintained the flow speed within a reasonable range for microscope
observations while using a decent fraction of the pressure range.

Microscope imaging and three-dimensional flow reconstructions. The
microfluidic flows were visualized by mixing the fluid with polystyrene beads
($2$ $\mu$m in diameter) as seeding particles, imaged under an inverted
microscope (Nikon Eclipse Ti2) at a $60\times$ magnification, operated in its
phase-contrast mode. These phase-contrast images were recorded by a USB
Scientific CMOS (sCMOS) video camera (Andor Zyla 4.2) in full resolution
($2048\times 2048$ pixels) at 50 fps. The diffraction pattern of each seeding
particle was visible within a range of $160$ $\mu$m along the optical axis,
which were calibrated to restore the axial positions of all seeding particles
in the view [35], forming a sizable (at least $200$ $\mu$m $\times$ $200$
$\mu$m $\times$ $100$ $\mu$m) visible zone of the 3D flow field.

## References

## References

  * [1] Squires, T. M. & Quake, S. R.  Microfluidics: Fluid physics at the nanoliter scale.  _Reviews of Modern Physics_ 77, 977–1026 (2005).
  * [2] Whitesides, G. M.  The origins and the future of microfluidics.  _Nature_ 442, 368–373 (2006).
  * [3] Mark, D., Haeberle, S., Roth, G., von Stetten, F. & Zengerle, R.  Microfluidic lab-on-a-chip platforms: requirements, characteristics and applications.  _Chemical Society Reviews_ 39, 1153 (2010).
  * [4] Battat, S., Weitz, D. A. & Whitesides, G. M.  Nonlinear Phenomena in Microfluidics.  _Chemical Reviews_ 122, 6921–6937 (2022).
  * [5] Manz, A. _et al._ Planar chips technology for miniaturization and integration of separation techniques into monitoring systems.  _Journal of Chromatography A_ 593, 253–258 (1992).
  * [6] Mao, X., Waldeisen, J. R. & Huang, T. J.  “microfluidic drifting”–implementing three-dimensional hydrodynamic focusing with a single-layer planar microfluidic device.  _Lab on a Chip_ 7, 1260 (2007).
  * [7] Drescher, K., Shen, Y., Bassler, B. L. & Stone, H. A.  Biofilm streamers cause catastrophic disruption of flow with consequences for environmental and medical systems.  _Proceedings of the National Academy of Sciences of the United States of America_ 110, 4345–4350 (2013).
  * [8] Garstecki, P. _et al._ Formation of monodisperse bubbles in a microfluidic flow-focusing device.  _Applied Physics Letters_ 85, 2649–2651 (2004).
  * [9] Amstad, E., Datta, S. S. & Weitz, D. A.  The microfluidic post-array device: high throughput production of single emulsion drops.  _Lab on a Chip_ 14, 705–709 (2014).
  * [10] Perkins, T. T., Smith, D. E. & Chu, S.  Single Polymer Dynamics in an Elongational Flow.  _Science_ 276, 2016–2021 (1997).
  * [11] Shenoy, A., Rao, C. V. & Schroeder, C. M.  Stokes trap for multiplexed particle manipulation and assembly using fluidics.  _Proceedings of the National Academy of Sciences of the United States of America_ 113, 3976–3981 (2016).
  * [12] Kaigala, G. V., Lovchik, R. D. & Delamarche, E.  Microfluidics in the “open space” for performing localized chemistry on biological interfaces.  _Angewandte Chemie International Edition_ 51, 11224–11240 (2012).
  * [13] Bhargava, K. C., Thompson, B. & Malmstadt, N.  Discrete Elements for 3D Microfluidics.  _Proceedings of the National Academy of Sciences of the United States of America_ 111, 15013–15018 (2014).
  * [14] Au, A. K., Huynh, W., Horowitz, L. F. & Folch, A.  3D-Printed Microfluidics.  _Angewandte Chemie International Edition_ 55, 3862–3881 (2016).
  * [15] Takahashi, K., Kamiya, S., Takao, H., Shimokawa, F. & Terao, K.  Stainless microfluidic probe with 2D-array microapertures.  _AIP Advances_ 11, 015331 (2021).
  * [16] Coxeter, H. S. M.  _Regular polytopes_ (Dover Publications, New York, 1973), 3d ed edn.
  * [17] Gross, D.  The role of symmetry in fundamental physics.  _Proceedings of the National Academy of Sciences of the United States of America_ 93, 14256–14259 (1996).
  * [18] Zandi, R., Reguera, D., Bruinsma, R. F., Gelbart, W. M. & Rudnick, J.  Origin of icosahedral symmetry in viruses.  _Proceedings of the National Academy of Sciences of the United States of America_ 101, 15556–15560 (2004).
  * [19] Purcell, E. M.  Life at low Reynolds number.  _American Journal of Physics_ 45, 3–11 (1977).
  * [20] Lauga, E. & Powers, T. R.  The hydrodynamics of swimming microorganisms.  _Reports on Progress in Physics_ 72, 096601 (2009).
  * [21] Berke, A. P., Turner, L., Berg, H. C. & Lauga, E.  Hydrodynamic Attraction of Swimming Microorganisms by Surfaces.  _Physical Review Letters_ 101, 038102 (2008).
  * [22] Elfring, G. J. & Lauga, E.  Hydrodynamic Phase Locking of Swimming Microorganisms.  _Physical Review Letters_ 103, 088101 (2009).
  * [23] Goyette, P.-A. _et al._ Microfluidic multipoles theory and applications.  _Nature Communications_ 10, 1781 (2019).
  * [24] Gonzalez, J. & Liu, B.  Symmetry-based nonperturbative micromanipulation in a three-dimensional microfluidic device.  _Physical Review Fluids_ 5, 044202 (2020).
  * [25] Fang, W., Xiong, T., Pak, O. S. & Zhu, L.  Data-Driven Intelligent Manipulation of Particles in Microfluidics.  _Advanced Science_ 2205382 (2022).
  * [26] Ashkin, A., Dziedzic, J. M., Bjorkholm, J. E. & Chu, S.  Observation of a single-beam gradient force optical trap for dielectric particles.  _Optics Letters_ 11, 288 (1986).
  * [27] Wu, J.  Acoustical tweezers.  _The Journal of the Acoustical Society of America_ 89, 2140–2143 (1991).
  * [28] Grier, D. G.  Optical tweezers in colloid and interface science.  _Current Opinion in Colloid & Interface Science_ 2, 264–270 (1997).
  * [29] Luan, Q., Macaraniag, C., Zhou, J. & Papautsky, I.  Microfluidic systems for hydrodynamic trapping of cells and clusters.  _Biomicrofluidics_ 14, 031502 (2020).
  * [30] Tu, M. Q., Nguyen, H. V., Foley, E., Jacobs, M. I. & Schroeder, C. M.  3D manipulation and dynamics of soft materials in 3D flows.  _Journal of Rheology_ 67, 877 (2023).
  * [31] Shenoy, A., Kumar, D., Hilgenfeldt, S. & Schroeder, C. M.  Flow Topology During Multiplexed Particle Manipulation Using a Stokes Trap.  _Physical Review Applied_ 12, 054010 (2019).
  * [32] The stress distribution is intrinsic to the far-field hyperbolic topology because any offset from the hyperbolic point, such as with a minimal planar extensional flow [36], is equivalent to adding a mean flow that imposes zero stress. For instance, taking that example flow velocity $\mathbf{v}\propto x\hat{x}-y\hat{y}$, the strain rate is found to be $\dot{\gamma}\propto\hat{x}\hat{x}-\hat{y}\hat{y}$. A displaced field, i.e. $\mathbf{v}\propto(x-a_{x})\hat{x}-(y-a_{y})\hat{y}$ for some physical displacement of the center $(a_{x},a_{y})$, is trivially seen to separate into the original field with an added mean flow, and the strain rate is identical and likewise invariant.
  * [33] Pal, P. B.  _A physicist’s introduction to algebraic structures: vector spaces, groups, topological spaces and more_ (Cambridge University Press, Cambridge, United Kingdom, 2019).
  * [34] Spencer, A. J. M.  _Continuum mechanics_.  Longman mathematical texts (Longman, London ; New York, 1980).
  * [35] Taute, K., Gude, S., Tans, S. & Shimizu, T.  High-throughput 3D tracking of bacteria on a standard phase contrast microscope.  _Nature Communications_ 6, 8776 (2015).
  * [36] Batchelor, G. K.  _An introduction to fluid dynamics_ (Foundation Books, Cambridge, U.K., 2007).  OCLC: 928892326.

Figure 1: Trap-associated v.s. trap-free manipulations of microscale particle
through surrounding flows. (a) In trap-associated manipulations, a manipulated
particle is led by a fluid “trap” or guided by the local flow director
(zoomed-in 2D views (i) and (ii) in the inset, respectively) along the desired
pathway (shown in the thick curve), surrounded by overall divergent flows
(with speeds of fluid tracers shown in colors) due to its trapping nature. The
stress distribution around the particle (with dashed circles showing target
positions) is intrinsic to the far-field flow topology (inset), with
compression (in red) and extension (in blue) along orthogonal axes,
independent of its relative position to the trap. (b) An ideal perturbation-
free manipulation is achieved by entraining the particle in a uniform flow
field. All surrounding fluid flows in parallel paths at identical velocities
and leaves the entrained particle stress-free (in blank color).  Figure 2:
Functional classification of a cross-channel microfluidic junction. (a) A 2D
CFD simulation of the flow pattern (with streamlines in solid curves and speed
$|\mathbf{u}|$ in colors) from an arbitrary flux setting ($f_{1}$ \- $f_{4}$
at four ports) can be decomposed into three linearly independent flow modes:
(b) two orthogonal displacement flow modes, and (c) one elongation flow mode.
(d) The flow pattern in (a) is reproduced by a flow potential $\phi$
associated with fluid sources (or sinks) at the four corners of a square
($q_{1}$ \- $q_{4}$). The square symmetry gives rise to (e) two degenerate
dipole modes equivalent to two displacement modes in orthogonal directions,
and (f) one quadrupole mode equivalent to elongation flow. Figure 3:
Microfluidic devices with 3D symmetries. (a) Fluid sources (or sinks)
following a tetrahedral symmetry ($q_{1}$ \- $q_{4}$) can be classified into
three degenerate functions, (b) a mix of dipole $p$ and quadrupole $Q$ in each
of the three orthogonal directions, representing simultaneous displacement and
elongation of the fluid. (c) Fluid sources (or sinks) following an octahedral
symmetry ($q_{1}$ \- $q_{6}$) can be classified into two subgroups: (d) three
dipoles along orthogonal directions ($x$, $y$ and $z$) and (e) two quadrupoles
in orthogonal planes (e.g., $x$-$y$ and $x$-$z$), contributing to separate
displacement and elongation functions. Figure 4: Realization of the octahedral
symmetry in 3D microfluidics. (a) An octahedron viewed along its original
$(1,1,1)$ axis (inset) is converted into a microfluidic junction with all
sources (or sinks) distributed at two planes perpendicular to the imaging axis
$z^{\prime}$. (b) The dipole flows in microscope coordinates
$(x^{\prime},y^{\prime},z^{\prime})$ are obtained by linear combination of the
three orthogonal dipole modes (Fig. 3(d)), with their experimental
realizations shown in (c), (d), and (e), respectively. The axial $z^{\prime}$
positions of the seeding particles are obtained by correlating their phase-
contrast images (projected in the bottoms of volumes) with a calibrated
library of diffraction patterns.  Figure 5: Dynamic stress-free displacement
flows in 2D. (a) Stress-free displacement flows were controlled in a phase
space $(c_{1},c_{2},c_{3})$, as the components of three orthogonal dipole
modes. Traces in the $c_{1}+c_{2}+c_{3}=0$ plane of this phase space
correspond to in-plane flows (in the $x^{\prime}$-$y^{\prime}$ plane) in the
real space. The dashed lines correspond to instantaneous transitions between
discrete states (dots) of the flow with the arrow showing the transition
directions, forming a triangle, a square, and a circle in the continuous
limit. The corresponding motions of the flow are shown in (b), (c), and (d),
respectively, with individual seeding particles shown in different colors. The
periods of these patterns are $3$ s, $4$ s, and $4$ s, respectively. Figure 6:
3D stress-free and multiplexed manipulation. (a) A design of the 3D Lissajous
pattern (inset, with 3D positions in color and 2D projections in gray) was
realized in the flow by incorporating additional axial motion into the
circular pattern [Fig. 5(d)], with all seeding particles in the view tracing
out the desired 3D pattern (in $8$ s). (b) Visualized trajectories (gray
curves) of seeding particles (11 particles shown) within the 3D view (200 $\mu
m$ $\times$ 200 $\mu m$ $\times$ 150 $\mu m$ in $x^{\prime}\times
y^{\prime}\times z^{\prime}$) were overlaid and exhibited a fair degree of
uniformity. (c) A characteristic strain $\epsilon$ computed along the
manipulation path further confirms the uniformity (with $\epsilon$ bounded by
$2\%$) of the manipulation flows.  Figure 7: More complex 3D flow
manipulations were demonstrated by rapidly imprinting the “UCM” (with its
design shown in the inset) onto the focal plane by seeding particles, which
are color-coded in time (over $23$ s). Some of those particles initially out
of the focal plane imprint the connecting lines (blue segments in the inset)
between letters. The scale bar shows 20 $\mu m$.


# Image storage on synthetic DNA using autoencoders

###### Abstract

Over the past years, the ever-growing trend on data storage demand, more
specifically for ”cold” data (rarely accessed data), has motivated research
for alternative systems of data storage. Because of its biochemical
characteristics, synthetic DNA molecules are now considered as serious
candidates for this new kind of storage. This paper presents some results on
lossy image compression methods based on convolutional autoencoders adapted to
DNA data storage.

The model architectures presented here have been designed to efficiently
compress images, encode them into a quaternary code, and finally store them
into synthetic DNA molecules. This work also aims at making the compression
models better fit the problematics that we encounter when storing data into
DNA, namely the fact that the DNA writing, storing and reading methods are
error prone processes. The main take away of this kind of compressive
autoencoder is our quantization and the robustness to substitution errors
thanks to the noise model that we use during training.

Index Terms—  image, compression, DNA, autoencoder, denoising

##  1 Introduction

Autoencoders has become a focus in the field of image compression. Quite
early, some autoencoders [1],[2] and recurrent neural networks[3] were already
showing better performances than the more conventional compression algorithms
like JPEG or JPEG2000. Some research [4], [5] aimed at having variable bit-
rate compression autoencoders, other work put their focus on the
differentiability problem of the quantization [1], [6] or the entropy
evaluation [1],[7]. Our goal is to adapt compression autoencoders to DNA data
storage. For that, we have based our proposed network on the compression
autoencoder designed by Lucas Theis et al. in [1], and introduced a DNA coder
for the autoencoder’s latent space tensor.

##  2 Context

###  2.1 DNA data storage

The memory of humanity relies on our ability to manage increasingly large
amounts of data, over periods of time ranging from a few years to several
centuries. Current tools are no longer sufficient and it is necessary to
consider game-changing solutions that can become operational quickly. One of
the most promising solutions is to store information in the form of DNA, just
like the genome by living beings. Indeed, DNA provides a very stable storage
medium over very long periods with simple implementation conditions [8].

The first step in the encoding workflow is the construction of a dictionary of
codewords composed by the symbols A, G, C and T also called nucleotides. The
DNA coded information stream must respect some biochemical constraints on the
combinations of bases that form a DNA fragment: homopolymers, high/low GC
content and repeated patterns should be avoided. One must also take into
account that the process involves some biochemical procedures which can
corrupt the encoded data. Synthesis, sequencing, storage and the manipulation
of DNA (mainly PCR amplification) may introduce errors by introducing
substitutions or indels (insertions or deletions of nucleotides), and may
jeopardize the integrity of the stored content [9].

However, to ensure better adaptation to the characteristics of the storage
medium, i.e., DNA, and possibly achieve higher storage efficiencies, it is
better to design coding algorithms specific for DNA storage. Indeed, since DNA
synthesis cost is relatively high, it is also important to take full advantage
of an optimal compression that can be achieved before synthesizing the
sequence into DNA. As an example, among other relevant works, [10] proposed
Discrete Wavelet Transform (DWT) image decomposition where the DWT
coefficients are scalar quantized and then encoded using a quaternary code. On
the contrary, following the example of what is done today in the context of
JPEG standardization for image coding with the development of JPEG
AI111https://jpeg.org/jpegai/index.html, we propose in this work to use neural
networks such as autoencoders for image coding on quaternary codes.

###  2.2 Compressive autoencoders for DNA image storage

Fig. 1: General scheme of an autoencoder used for image storage in synthetic
DNA. $Z$ is the latent vector in the latent space. For compression purpose,
$Z$ is quantized and encoded in a quaternary code for DNA storage. Then the
quaternary code is used for DNA synthesis. The image reconstruction process
needs sequencing before decoding with the autoencoder.

The general idea of our proposed encoding process is depicted in the figure
(1) and can be very roughly described by the following steps. Firstly, the
input image has to be compressed using a lossy/near lossless image coder.
Here, we propose to use a compressive autoencoder to learn the image
characteristics as well as the biochemical noisy process. Convolutionnal
autoencoders have shown good properties for image denoising and could then
appear as a good solution when it comes to deal with the noise introduced by
the DNA data storage process. Secondly, a quantization operation is introduced
to quantize the latent vector $Z$ in the latent space. The reason for that is
that the latent space doesn’t suffice as a compressed representation because
it uses floating numbers. The quantization goal is to compress it and reduce
its cost on the one hand and to facilitate its encoding on the other hand.
Thirdly, in the case of DNA storage, the latent vector is encoded using a
quaternary code using the $\\{A,C,T,G\\}$ alphabet. Decoding is ensured thanks
to the decoder. The quantization operation is one of the major problems when
training compressive autoencoders. Indeed, to be able to train a model, all
the operations have to be differentiable, which is not the case for the
quantization (most of the times it is equivalent to a rounding function).
However, this can be tackled by the use of a linear approximation function as
in [1].

During training, the compressive autoencoder has to minimize two quantities:

  * •

The distortion which corresponds to the squared difference between the
original image and the image reconstructed after compression and
decompression.

  * •

The entropy of the quantized latent space computed in a quaternary basis. The
entropy being closely linked to the rate, we will use indifferently both terms
in the rest of the paper.

The structure of the compressive autoencoder model we propose in this paper is
closely linked the one proposed by Theis et al. in [1]. It is described in the
following section 3. The introduction of the biochemichal noise in the model
is described in the section 4.

##  3 Non-noisy latent space

###  3.1 Neural Network

We designed our proposed autoencoder to obtain a latent space with higher
dimensionality than the one proposed in [1], the objective of such a model was
to obtain higher bit-rates. Because we did not reduce the size of the latent
space as much as in [1], our new model is shallower. This means that our new
model now has a reduced number of parameters, which could cause some
overfitting problems during the training process. To compensate such a risk,
we decided to add new residual blocks to the model with skipped connections.

When encoding data, the common practice is to operate with integers, not with
floating point numbers. For that reason, the output $Z$ of an autoencoder
group of layers in the latent space cannot be encoded using quaternary code as
it is and a quantization operator must then be introduced in the process.
Given an input image $I$, the output image (decoded image) $\hat{I}$ of the
autoencoder can be expressed as a combination of several functions such as
follows:

$\hat{I}=g(\beta_{DNA}(\alpha_{DNA}(Q(f(I))))$ (1)

where, $f$ represents the encoding part of the autoencoder and $g$ the
decoding, $Q$ represents the quantizer that rounds the components z of the
latent vector $Z=f(I)$ into integer values and is defined as:

$Q(z)=q\times\left\lfloor\frac{z}{q}+\frac{1}{2}\right\rfloor$ (2)

Finally, $\alpha_{DNA}$ and $\beta_{DNA}$ represent the DNA-encoding and
decoding algorithms that encode a sequence of symbols into a quaternary stream
composed by the letters of the alphabet $\\{A,C,T,G\\}$. When no noise, it is
clear that $\beta_{DNA}(\alpha_{DNA}(k))=k\ \forall k\in\mathbb{Z}$. In this
work we have used the DNA fixed-length code proposed in [10].

###  3.2 Loss function

We used in our work the following classical loss function:

$L=||I-\hat{I}||^{2}+\lambda H(Q(Z)),$ (3)

where $I$ and $\hat{I}$ are respectively the input and output images and
$H(Q(Z))$ corresponds to the entropy of the quantized latent space computed in
base $4$ as follows:

$H(Q(Z))=-\sum_{i=1}^{n}Pr\\{Q(z_{i})\\}log_{4}Pr\\{Q(z_{i})\\}$ (4)

with $Q(Z)=(Q(z_{1}),Q(z_{2}),...,Q(z_{n}))$ and $Pr\\{Q(z_{i})\\}$ the
probability of a quantized component. The entropy $H$ is expressed in
nucleotides per component.

To encode the quantized values, we use a fixed-length coding system. Every
value is coded with a DNA codeword of the same length $n$. The set of all
different codewords of length $n$ is a finite set. In the case of quaternary
code, there are $4^{n}$ different codewords, but only a subpart are DNA
codewords that respect biochemical constraints (mainly no codes containing
homopolymers runs such as AAAA or TTTT for example can be generated). We call
a codebook the set of constrained DNA codewords of length $n$ that can be used
for this fixed-length encoding. Ensuring encodability of any compressed image
means that the number of different possible values output by the quantizer has
to be smaller or equal to the number of codewords available.

The compression model output was bounded using a hyperbolic tangent function.
The uniform quantization was then applied to the output of that bounded
function, giving a finite number of different possible values. The number of
possible values can be adjusted with the quantization step.

##  4 Introducing biochemical noise

As described previously, one must also take into account the process involves
some biochemical procedures which can corrupt the encoded data. Synthesis,
sequencing, storage and the manipulation of DNA may introduce errors by
introducing substitutions or indels (insertions or deletions of nucleotides
$A,T,C$ or $G$), and may jeopardize the integrity of the stored content [9].
DNA storage can then be viewed as a naturally noisy channel for which
appropriately resilient encoding solutions need to be defined.

In order to adapt the autoencoders to the noisy channel of DNA data storage,
one needs to introduce a noise model between the encoding and decoding parts.
Since the substitution noise is prevalent in the noise of the DNA storage
channel, we decided in this work to focus on it. A substitution, as mentioned
earlier, is the phenomenon of one nucleotide being changed into another one
(see figure 2).

Furthermore, since we are using the fixed length codes of [10] for encoding
the quantized values, a substitution error affects only one code and has no
effect on the decoding of its neighbors as shown in figure 2. This means that
the model for the substitution noise can be established and applied at the
quantized tensor level (coefficients in the latent space) and not necessarily
at the nucleotide level.
Let’s call $\eta$ the substitution noise introduced by the biochemical
process. One can rewrite the autoencoder input/output function given in
equation (1) as:

$\hat{I}=g(\overline{\beta_{DNA}}(\alpha_{DNA}(Q(f(I)))+\eta)),$ (5)

where the noise is introduced at the level of the latent space and the
quaternary decoder $\overline{\beta_{DNA}}$ is designed to decode a noisy
code. Because of the biochemical constraints, the DNA code is not a quaternary
complete code and thus, a noisy code might not be decodable. To ensure
decodability, when a noisy codeword is not decodable, we replace it with the
closest valid possible code in the codebook, in terms of Hamming distance.

Fig. 2: Example of a substitution error on a DNA code. It only impacts one
decoded symbol and not its neighbors.

In this work we assumed that this noise can be modeled by a i.i.d Gaussian
process. The optimization of the autoencoder during training is still
performed by minimizing the loss function given by equation (3).

##  5 Experimental results

###  5.1 Implementation of the training

We trained the model with the 30k Flickr image dataset. During the training
step, we used batches of 32 random crops of size 96x96 from the Flickr images.
The training process has been separated into two steps: the first with a
learning rate of 1e-4 during 200 epochs and the second one of 1e-5 that would
be used for 500 epochs. Models have been trained independently for each
quantization step.

The model was then evaluated using the Kodak dataset. For each compression
rate chosen, we computed the performances for each image of the dataset, and
also the average on the whole dataset. In figure 4, the gray curve (called
avg) represents the average performance of the model and it’s resistance to
different levels of noise.

(a) Without noise optimization

(b) Optimized for noise

Fig. 3: Visual results on kodim18 and kodim21 compressed at 6 bits/nt with
noise levels of 5%. Left: Model trained without any substitution noise
considerations, right: training adapted to noise

###  5.2 Training process without substitution noise

In this section, we describe the performance of the networks trained as
described in section 3 using the model of formula (1), and how their
performance is maintained when noise is introduced into the channel at the
encoding step and not during the training. The experiments were conducted on
different autoencoder models, each one trained to a given compression rate
(given quantization step $q$ as defined in formula (2)).

The result are presented in the figure 4 where the PSNR is evaluated for
different noise levels and averaged on all the images of the Kodak dataset.
The results are provided for two different rates, 6 and 4 bits/nucleotide. The
reconstruction maintains a good visual quality until the substitution noise
level reaches around 5%, where a lot of artifacts start to affect drastically
the reconstruction, only maintaining an approximation of the general image
features as shown in the figure 3(a). What we can quickly understand from
those results is that the substitution noise has a big influence on the
performances of our models. After a few percentages of error, the images
technically become unusable. This justifies the development for compression
methods adapted to noise.

###  5.3 Training process including substitution noise

In this section, we evaluate the performance of our autoencoder models trained
with a noisy latent space as proposed in section 4 and formula (5). For each
quantization step $q$ (or equivalently each rate in nucleotides per component)
we trained a specific autoencoder.

Here, the most important parameter is the noise level. It increases throughout
the training from 0 (meaning that at the beginning of the training, the model
is learning with a non-noisy data) and a value $max\\_level$, which is the
maximum noise level, used during the last epochs of the training session. Note
that a more complex alternative would be to train several models each one for
different noise levels, instead of training a unique model for all the
possible noise levels.

Fig. 4: Comparison of robustness to noise at different coding rates for models
trained without considering noise substitution (avg) and model trained by
considering the substitution noise (avg-noise).

What can be taken out from the results shown in the figure 4 is that by
training our models to adapt to a substitution noise we managed to obtain a
gain between 0.5 and 1.5 dB on the PSNR depending on the level of noise. On
the other hand, the models adapted to noise seem to underperform when no noise
is introduced in the latent space. Furthermore, the visual results provided by
the optimized autoencoder remain very interesting for a substitution noise
level around 5% (see figure 3(b)) showing a strong robustness to high
substitution noise levels and the good performance of the proposed solution.

##  6 Conclusion

In this work, we have developed a compression solution for image storage on
synthetic DNA, robust to substitution noise. The proposed approach is based on
compressive autoencoder optimized for DNA fixed-length encoding technologies.
A noise model was developed and introduced in the autoencoder optimization to
analyze its effects on the quality of reconstruction of the decoded image.
Training the compression neural network including the noise model showed
improvements over the network trained without the noise model (figure 4).

In future works, experimenting with entropy-based DNA coding systems instead
of fixed-length encoding might show some improvements since the compression
network minimizes entropy. Introducing new types of noise (insertions and
deletions), and solutions to minimize their effects is also another field of
improvement.

## References

  * [1] Lucas Theis, Whenzhe Shi, Andrew Cunningham, and Ferenc Huszár,  “Lossy image compression with compressive autoencoders,”  International Conference on Learning Representations, 2017.
  * [2] Ballé Johannes, David Minnen, Singh Saurabh, Sun Jin Hwang, and Nick Johnston,  “Variational image compression with a scale hyperprior,”  International Conference on Learning Representations, 2018.
  * [3] George Toderici, Damien Vincent, Nick Johnston, Sun Jin Hwang, David Minnen, Joel Shor, and Michelle Covell,  “Full resolution image compression with recurrent neural networks,”  IEEE Conference on Computer Vision and Pattern Recognition, 2017\.
  * [4] Fei Yang, Luis Herranz, Joost van de Weijer, José A. Iglesias Guitián, Antonio M. López, and Mikhail G. Mozerov,  “Variable rate deep image compression with modulated autoencoder,”  IEEE Signal Processing Letters, vol. 27, pp. 331–335, 2020.
  * [5] Yoojin Choi, Mostafa El-Khamy, and Jungwon Lee,  “Variable rate deep image compression with a conditional autoencoder,”  IEEE/CVF International Conference on Computer Vision (ICCV), 2019\.
  * [6] Thierry Dumas, Aline Roumy, and Christine Guillemot,  “Autoencoder based image compression: can the learning be quantization independent?,”  2018\.
  * [7] V Oliveira, T Oberlin, M. Chabert, Charly Poulliat, Mickaël Bruno, C Latry, M Carlavan, S Henrot, F Falzon, and Roberto Camarero,  “Simplified entropy model for reduced-complexity end-to-end variational autoencoder with application to on-board satellite image compression,”  09 2020.
  * [8] S. M. Hossein Tabatabaei Yazdi, Ryan Gabrys, and Olgica Milenkovic,  “Portable and error-free dna-based data storage,”  Nature, 2017.
  * [9] N. Goldman, P. Bertone, and S. Chen,  “Towards practical, high-capacity, low-maintenance information storage in synthesized dna,”  Nature, 2013.
  * [10] Melpomeni Dimopoulou, Marc Antonini, Pascal Barbry, and Raja Appuswamy,  “A biologically constrained encoding solution for long-term storage of images onto synthetic dna,”  European Signal Processing Conference (EUSIPCO), 2019.