# Mathematical Model for Chemical Reactions in Electrolyte Applied to
Cytochrome $c$ Oxidase: an Electro-osmotic Approach

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

###### Abstract

A mathematical model for chemical reactions in electrolytes is developed using
an Energy variational method consistent with classical thermodynamics.
Electrostatics and chemical reactions are included in properly defined
energetic and dissipative functionals. The energy variation method is
generalized to deal with open systems with inputs of charge, mass, and energy.
The open systems can transform input energy of one type into output energy of
another type. The generalized method is used to analyze the conversion of
electrical current into proton flow by cytochrome c oxidase, an important
enzyme in mitochondria that helps generate the ’energetic currency of life’
ATP. The structure of the oxidase guides flows of current and mass that
interact according to the energetic and dissipative functionals of the
generalized theory. Kirchhoff’s current law provides important coupling when
the enzyme is in its natural setting in the mitochondrial membrane. The
natural function of the enzyme is the result. Electron flows are converted
into proton flows and gradients.

##  1 Introduction

History seems to have separated much of chemistry [25, 36] from the classical
theory of fields [42, 78, 2]. Chemical reactions are found throughout the
ionic solutions of biology and chemistry but they are usually described in a
language apparently disjoint from that of classical field theory even though
the reactants and products of chemical reactions are almost always charged and
carry significant electrical current. The reactants, catalysts and enzymes of
chemistry and biology depend on charge interactions for much of their
function.

Field theory has much to offer chemistry, particularly in the study of charged
systems as admirably reviewed in [115]. The Maxwell equations are as universal
and precise as any theory and apply to chemical reactions in the plasmas of
gases and ionic solutions, and liquids in general. Indeed, the Maxwell
equations can be written without any adjustable parameters, implying the
universal conservation of total current [91, 97]. A theory of the
electromechanical response of charge to the electric field is needed in that
case to make a complete description of the charged system, i.e. an
electromechanical theory of polarization phenomena [109]. But some properties
of the electromagnetic field (e.g., conservation of total current) are true
entirely independent of the electromechanical response. The special systems
that are almost completely described by conservation of current (without
specification of charges) include the electronic circuits of our computers and
many properties of the action potential of nerve and muscle fibers. Both
systems are almost entirely specified by Kirchhoff’s current law. The question
then arises how do we fit chemical reactions into the framework of Kirchhoff’s
law and conservation of total current that is so remarkably simpler to
implement than a full accounting of all the charges in a chemical system?

Here we show one way to describe chemical reactions in ionic solutions with an
extension of classical field theory that does not violate the traditions of
either chemistry or electrodynamic field theory. In this work, we take
advantage of the energy variation method [90, 103, 113] that treats ionic
solutions as complex fluids, with interactions, internal stored energy, flow,
and dissipation like other complex fluids [70, 96, 101, 99, 85]. It begins
with defining two functionals for the total energy and dissipation of the
system, and introducing the kinematic equations based on physical laws of
conservation. The specific forms of the flux and stress functions in the
kinematic equations are obtained by taking the time derivative of the total
energy functional and comparing with the defined dissipation functional. More
details of this method can be found in [103]. We use energy variation methods
to link the electric field and reaction dynamics as Wang, et al, [104, 115]
have for reactions that do not involve charge or electrodynamics.

The generalization to charged systems allows us to study systems of some
importance. We study active transporters of biological membranes. In a sense,
we extend the electrical treatment of the membrane proteins called channels to
deal with transporters. Conservation of current (in the form of Kirchhoff’s
law) provides the crucial coupling between the properties of disjoint sodium
and potassium channel proteins that act independently in the atomic and
molecular sense because they are well screened. Channels are many Debye
lengths apart, shielded from each other by the ions, water dipoles (and
quadrupoles), that also form the ionic atmosphere of proteins and lipid
bilayer. The atomic scale function of these proteins is crucial to their
biological function. Their coupling is just as important to their function,
but the coupling of these channel proteins is not chemical. The coupling is
provided by the cable equation [24] as biophysicists have called the
telegrapher equation version of the Maxwell equations (and Kirchhoff’s law)
used by Kelvin to design the Atlantic cable [1] well before Maxwell wrote his
equations.

Active transport is one of the most important processes in life. It maintains
the concentration gradients and membrane potentials that allow biological
cells to function. Indeed, without active transport animal cells swell, burst,
and die. Active transport powers the generation of ATP in both animals (in
mitochondria, by oxidative phosphorylation) and plants (in chloroplasts, by
photosynthesis). ATP is the currency of chemical energy in all plants and
animals, storing in a small organic compound the energy from photosynthesis or
oxidation. When hydrolyzed to ADP, the chemical energy is available for the
myriad of dissipative processes essential to life. Nearly all of them use ATP
as their energy source. Life is complex with many facets. The energy source of
life is not.

The coupling of flows in transporters in the inner membrane of mitochondria
allows one substance to move uphill (against its gradient of electrochemical
potential) using the energy from the downhill movement of another substance.
Coupling is inherently about the relationship of flows, yet most analysis and
simulations of active transporters do not explicitly include a variable for
flow. Most do not use the electrodynamic equations for flow, whether
conservation of total current or Kirchhoff’s law or a nonlinear version of
Ohm’s law. Indeed, most analysis and simulations use methods that are derived
assuming zero flow and do not include changes in potential or free energy
associated with flow. The problems with this approach become apparent if one
tries to analyze electron flow through a resistor or semiconductor diode or
rectifier in that tradition.

We analyze an active transporter by studying the currents through it, as a
specific example of our general field theory for ionic flows with chemical
reactions. We combine Kirchhoff’s law and chemical reaction energetics with
diffusion, migration, and convection of ions and water to make an ’electro-
osmotic’ model. The name is chosen to emphasize the important role of
electricity in this system, implying the need to deal with electrical flows by
methods used to deal with electrical flows in other systems, like electrical
and electronic circuits. This approach seems sensible for the systems like
oxidative phosphorylation and photosynthesis where electron flows are
involved. The analysis of electron flows has been well established in physics
and engineering for more than a century. The analysis of ionic flows has been
well established in membrane biophysics for some seventy years. We combine
them here with chemical reactions hoping to construct a useful electro-osmotic
theory of cytochrome $c$ oxidase.

We do not have to deal with the myriad of charges involved in the transport of
current or the incalculable number of interactions of those charges. Analysis
of current flow is all we use in this conservative coupling approach,
following the practice of circuit analysis. We do not have to assume
equilibrium or zero flow. We do not have to deal explicitly with the charges
involved in the transport. Analysis of current flow is all we use. In
particular, we do not assume equilibrium or near equilibrium flows, as has
been done in previous analysis. Note that near equilibrium analysis (using the
Green-Kubo formalism, for example) is inappropriate for devices that function
with large flows. These devices are not near equilibrium. The electronic
devices of our digital technology function far from equilibrium and they are
not analyzed by assuming nearly zero flow. Indeed, they usually use power
supplies to maintain spatially nonuniform potentials, often described by
inhomogeneous Dirichlet boundary conditions. Traditionally, electronic devices
are analyzed by studying small changes around a nonequilibrium operating
point, which we hasten to add is maintained by large—not small— flows from the
power supplies. But even that linearization is not necessary nowadays because
the full flow nonequilibrium problems can be solved conveniently with readily
available software.

We are certainly not the first to exploit the simplification provided by the
analysis of current instead of charge. Circuit designers have used this
approach ’forever’ [4, 5]. Charges are hardly mentioned when circuits are
designed as a glance at textbooks shows [77, 80, 40, 100, 49, 54, 102],
perhaps most eloquently in symbolic circuit design [73]. Analysis of
current—not charge—characterizes the study of ion channels since they were
discovered as conductances by Hodgkin and Huxley some seventy years ago [32,
8, 14, 47, 6, 7]. Analysis of current is not a prominent feature of the study
of active transporters, however, even in the work of Hodgkin and collaborators
that started the physiological analysis of active transport in cell membranes
(not whole epithelia or biochemical preparations), at much the same time as
their work defining channels [10, 9, 15].

Most models of active transport include conformational changes of the protein
that require a model to compute the spatial distribution of mass in the
protein, as it changes during active transport [52, 57]. The conformational
changes usually provide alternating access to an occluded state that is not
connected (or accessible) to either side of the protein. The occluded state
blocks the conduction path (and incidentally often traps ions in the ’middle’
of the transporter) and thus prevents backflux. Alternating access mechanisms
create flux across the transporter protein without allowing backflux that
would seriously degrade the efficiency of the transporter. Transporters were
first thought to use quite different mechanisms from channels [29, 28].
However, recent work [31] shows otherwise. Alternating access is apparently
created by correlated motions of gates that account for activation and
inactivation in classical voltage activated sodium channels [52, 64]. The
physical basis of the gates is not discussed in the classical literature.
Later in the paper, we speculate that the switches that provide alternating
access might be like the switches of bipolar transistors. .

In this work, we use an electro-osmotic approach to describe cytochrome $c$
oxidase (or Complex IV). We choose cytochrome $c$ oxidase because experimental
work and simulations of the highest quality have shown that “cytochrome $c$
oxidase is a remarkable energy transducer [i.e., coupled transporter of
electrons and protons] that seems to work almost purely by Coulombic
principles without the need for significant protein conformational changes”
[89]. Cytochrome $c$ oxidase depends on an "occluded state" containing the
reaction center(s) to prevent backflow as do other active transporters, but it
uses some type of the ’water-gate’ mechanism proposed in [50, 89]. The
alternating access in cytochrome $c$ oxidase occurs without conformation
change (of the spatial distribution of mass), in marked contrast to the usual
alternating access models of transporters. Perhaps the gate in the oxidase is
like the switch in a semiconductor (diode) rectifier. The switches might be
rectifiers produced by spatial distributions of permanent charge, of opposite
signs, as rectification (and switching) is produced in PN diodes, and bipolar
transistors. Diode rectifiers depend on changes in shape (i.e., conformation)
of the electric field, not changes in the distribution (i.e., conformation) of
mass. This idea is outlined in the Discussion Section following [38, 37, 59].

The rest of the paper is organized as follows. In section 2, we derive the
general three dimensional field equations for an ionic system with reaction
and use it to create a general framework of an electro-osmotic model. In
section 3, we propose a specific, simplified model for cytochrome $c$ oxidase.
In section 4, we carry out computational studies of our cytochrome $c$ oxidase
model and explore the effects of various conditions on the transport of
protons across the mitochondria membrane. In section 5, we conclude our paper
with a discussion of our cytochrome $c$ oxidase model and future directions.

##  2 Derivation of Electro-osmotic Model

We mainly focus on a mathematical model of elementary reactions

$\alpha_{1}C_{1}^{z_{1}}+\alpha_{2}C_{2}^{z_{2}}+\alpha_{3}C_{3}^{z_{3}}\overset{k_{f}}{\underset{k_{r}}{\rightleftharpoons}}\alpha_{4}C_{4}^{z_{4}},$
(1)

where $k_{f}$ and $k_{r}$ are two constants for forward and reverse
directions, $[C_{i}]$ is the concentration of $i^{th}$ species, respectively.
Here $\alpha_{i}$ is the stoichiometric coefficient, $z_{i}$ is the valence of
$i^{th}$ species and together they satisfy

$\sum_{i=1}^{3}\alpha_{i}z_{i}=\alpha_{4}z_{4}.$ (2)

In particular, we have in mind a case where an active transporter (’pump’)
uses the energy supplied by a chemical reaction to pump molecules. Later, we
will focus on the reaction for cytochrome $c$ oxidase, i.e., for Complex IV of
the respiratory chain

$2H^{+}+\frac{1}{2}O_{2}+2e^{-}\overset{k_{f}}{\underset{k_{r}}{\rightleftharpoons}}H_{2}O.$
(3)

According to the conservation laws, we have the following conservation of
chemical elements (like sodium, potassium and chloride). Note that this
conservation is in addition to the conservation of mass, because nuclear
reactions that change one element into another are prohibited in our
treatment, as in laboratories and most of life.

$\displaystyle\frac{d}{dt}(\alpha_{4}[C_{1}]+\alpha_{1}[C_{4}])=0,$ (4a)
$\displaystyle\frac{d}{dt}(\alpha_{4}[C_{2}]+\alpha_{2}[C_{4}])=0,$ (4b)
$\displaystyle\frac{d}{dt}(\alpha_{4}[C_{3}]+\alpha_{3}[C_{4}])=0.$ (4c)

In order to derive a thermal dynamical consistent model, the Energy Variation
Method [103] is used. Based on the laws of conservation of elements and
Maxwell equations, we have the following kinematic system

$\left\\{\begin{array}[]{l}\frac{d[C_{1}]}{dt}=-\nabla\cdot\bm{j}_{1}-\nabla\cdot\bm{j}_{p}-\alpha_{1}\mathcal{R},\\\
\frac{d[C_{2}]}{dt}=-\nabla\cdot\bm{j}_{2}-\alpha_{2}\mathcal{R},\\\
\frac{d[C_{3}]}{dt}=-\nabla\cdot\bm{j}_{3}-\alpha_{3}\mathcal{R},\\\
\frac{d[C_{4}]}{dt}=-\nabla\cdot\bm{j}_{4}+\alpha_{4}\mathcal{R},\\\
\nabla\cdot(\bm{D})=\sum_{i=1}^{4}z_{i}[C_{i}]F,\\\
\nabla\times\bm{E}=\bm{0},\end{array}\right.$ (5)

where $\bm{j}_{l},l=1,2,3,4$ are the passive fluxes and $\bm{j}_{p}$ is the
pump flux, $\mathcal{R}$ is reaction rate function. All these variables are
unknown and will be derived by using the Energy Variational method.

$\bm{j}_{ex}$ is the flux of electrons supplied from an external source. In
mitochondiral membranes this will include special pathways linking one enzyme
and one complex to another by the movement of lipid soluble or water soluble
electron donors and water soluble electron donors or acceptors like the
quinones.

$\bm{D}$ is Maxwell’s electrical displacement field and
$\bm{D}=\varepsilon_{0}\varepsilon_{r}\bm{E}$ with electric field $\bm{E}$,
dielectric constant $\varepsilon_{0}$ and relative dielectric constant
$\varepsilon_{r}$. The equation $\nabla\times\bm{E}=\bm{0}$ implies that there
exists a $\phi$ such that $\bm{E}=-\nabla\phi$. We consider a system with
structure and boundary conditions defined on that structure.

The structures are given to us by structural biologists. The structures are
decorated with molecules (proteins and lipids for the most part) that use
particular atomic arrangements to channel physical forces into physiological
function. We describe a constant flux for one species, say $C_{3}$, that serve
as the input of the system. In cytochrome $c$ oxidase the input is electrons
carried on the heme groups of cytochrome oxidase.

$\left\\{\begin{array}[]{ll}\bm{j}_{i}\cdot\bm{n}=j_{i,extra},i=1\cdots,4,&\mbox{on~{}}\partial\Omega,\\\
\bm{D}\cdot\bm{n}=0,&\mbox{on~{}}\partial\Omega.\end{array}\right.$ (6)

######  Remark 2.1.

By multiplying $z_{i}e$ on both sides of the first three equations and $-e$ on
both sides of the fourth equation, we have

$\displaystyle\frac{d}{dt}(\nabla\cdot\bm{D})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{4}z_{i}F\frac{d[C_{i}]}{dt}$ (7) $\displaystyle=$
$\displaystyle-\sum_{i=1}^{4}\nabla\cdot(z_{i}F\bm{j}_{i})-z_{i}F\nabla\cdot\bm{j}_{p}-(z_{1}\alpha_{1}+z_{2}\alpha_{2}+z_{3}\alpha_{3}-z_{4}\alpha_{4})F\mathcal{R}$
(8) $\displaystyle=$
$\displaystyle-\sum_{i=1}^{4}\nabla\cdot(z_{i}F\bm{j}_{i})-z_{i}F\nabla\cdot\bm{j}_{p},$
(9)

which is consistent with the electrostatic Maxwell equations. Treatment of
transient problems, involving displacement currents is needed to deal with
some important experimental work [56, 51, 55, 72, 87, 63].

The total energetic functional is defined as the summation of entropies of
mixing, [22] internal energy and electrical static energy.

$\displaystyle E_{tot}$ $\displaystyle=$ $\displaystyle
E_{ent}+E_{int}+E_{ele}$ (10) $\displaystyle=$
$\displaystyle\sum_{i=1}^{4}\int_{\Omega}RT\left\\{[C_{i}]\left(\ln{\left(\frac{[C_{i}]}{c_{0}}\right)}-1\right)\right\\}dx+\int_{\Omega}\sum_{i=1}^{4}[C_{i}]U_{i}dx+\int_{\Omega}\frac{\bm{D}\cdot\bm{E}}{2}dx.$

Then the chemical potentials can be calculated from the variation of total
energy

$\tilde{\mu}_{l}=\frac{\delta
E_{tot}}{\delta[C_{i}]}=RT\ln\frac{[C_{i}]}{c_{0}}+U_{i}+z_{l}\phi
e,l=1,\cdots,4.$ (11)

It is assumed in the present work that dissipation of the system energy is due
to passive diffusion, chemical reaction and the pump. Accordingly, the total
dissipation functional $\Delta$ is defined as follows

$\displaystyle\Delta=\int_{\Omega}\left\\{\sum_{j=1}^{4}|\bm{j}_{i}|^{2}+RT\mathcal{R}\ln\left(\frac{\mathcal{R}}{k_{r}\left(\frac{[C_{4}]}{c_{0}}\right)^{\alpha_{4}}}+1\right)\right\\}dx-\int_{\Omega}f_{p}dx,$
(12)

where $f_{p}=f_{p}(\mathcal{R},\mu,x)\geq 0$ is the term from the pump.

Open systems in which some fluxes flow in or out, entering or leaving the
system altogether, have distinctive energy dissipation laws that differ from
those of closed systems. The natural mitochondrion is an UNclamped system, in
which the electrical potential assumes whatever value satisfies the field
equations. The sum of all currents across the membrane of the natural
mitochondrion is zero (including the capacitive displacement current) as it is
in small biological cells. Many experiments are done in voltage clamped
systems. In these the sum of the currents does not equal zero just as the sum
of currents in the classical Hodgkin Huxley experiments was not zero. Of
course, the ratio of fluxes will be different in the clamped and unclamped
cases, as we document at length later in this paper.

In the natural Unclamped mitochondrion, we have the following generalized
energy dissipation law

$\frac{dE_{tot}}{dt}=J_{E,\partial\Omega}-\Delta.$ (13)

Here $J_{E,\partial\Omega}$ is the rate of boundary energy absorption or
release that measures the energy of flows that enter or leave the system
altogether through the boundary. Recall that the chemical potential of a
species is the energy that can be absorbed or released due to a change of the
number of particles of the given species and $J_{i}\cdot n$ is the total
number of $i^{th}$ particles passing through the boundary, per area per unit
time. We define $J_{E,\partial\Omega}$ as follows

$J_{E,\partial\Omega}=\int_{\partial\Omega}\sum_{i=1}^{4}\tilde{\mu}_{i}\bm{j}_{i}\cdot\bm{n}dS.$
(14)

In general, different types of boundary conditions can be written in the
following general format

$\bm{j}_{i}\cdot\bm{n}=g_{i}(f([C_{i}])-f([C_{i}]_{extra})),$ (15)

where $g_{i}$ is the conductance of $i^{th}$ species on the boundary,
$[C_{i}]_{extra}$ is the fix reservoir’s concentration of $i^{th}$ species and
$f$ is some specific function. Then the rate of boundary energy absorption or
release is

$J_{E,\partial\Omega}=\int_{\partial\Omega}\sum_{i=1}^{4}g_{i}\tilde{\mu}_{i}(f([C_{i}])-f([C_{i}]_{extra})).$
(16)

In this case energy can change both because of the flux across the boundary
and also because of the change in dissipation. and

$\frac{dE}{dt}-J_{E,\partial\Omega}=-\Delta.$ (17)

######  Remark 2.2.

Boundary Conditions, Structure, Evolution, and Engineers

These boundary conditions serve as the link between general field equations
and structures that serve as devices. Structures are chosen and devices
designed (by evolution or engineers) so these boundary conditions are
satisfied. The boundary conditions are chosen so devices have almost the same
properties no matter where they are placed in a network. The structures and
boundary conditions on those structures are not automatic properties of
nature. The structures are decorated with (i.e., include) specific
substructures (like power transistors) that exploit arrangements of atoms
(like doping charges) to create properties that are useful. The properties are
summarized by boundary conditions located on the structures provided by
evolution and engineers. These boundary conditions help make the idea of a
component useful. They help ensure that a component in one part of a system
does what it does in another part of the system and so can be described by a
’transfer function’ independent of its location in the system.

It is clear that channels and transporters in biological systems behave as
components. Indeed, most of classical physiology and biophysics is devoted to
identifying such components, on a wide variety of length scales from atoms to
organisms, and studying how they interact in the hierarchy of structures that
make animals and plants [61, 83, 79, 66, 93, 75, 35, 30].

######  Remark 2.3.

  1. 1.

A closed system allows no flux across the boundaries. It has the following no-
flux boundary conditions

$\left\\{\begin{array}[]{ll}\bm{j}_{i}\cdot\bm{n}=0~{}i=1,2,3,4,&\mbox{on~{}}\partial\Omega,\\\
\bm{D}\cdot\bm{n}=0,&\mbox{on~{}}\partial\Omega.\end{array}\right.$ (18)

In a closed system, $J_{E,\partial\Omega}=0$ and the energy dissipation law is

$\frac{dE_{tot}}{dt}=-\Delta.$

In a closed system, the energy changes into dissipation. That is the only way
energy can change in a closed system.

  2. 2.

An open system has flow across the boundaries. An open system might have
constant inflow/outflow $\bm{j}_{i}\cdot\bm{n}=J_{i,extra}$. In that case,
$g_{i}=\frac{J_{i,extra}}{f([C_{i}])-f([C_{i}]_{extra})}$, and

$J_{E,\partial\Omega}=\int_{\partial\Omega}\sum_{i}^{4}\tilde{\mu}_{i}J_{i,extra}dS.$

  3. 3.

For the Dirichlet boundary condition $[C_{i}]=[C_{i}]_{extra}$ on
$\partial\Omega$, the flux $\bm{j}_{i}\cdots\bm{n}$ is unknown and part of the
solution. In this case, $J_{E,\partial\Omega}$ is an unknown flux needed to
ensure that the Dirichlet condition $[C_{i}]=[C_{i}]_{extra}$ is obeyed on
$\partial\Omega$.

It is very important to understand this requirement. In reality, i.e., in
experiments and their models, supplying the unknown flux requires specialized
instrumentation, for example, a patch clamp amplifier in a voltage clamp
setup. Almost always, that flux is supplied at one location in space. In that
way a classical voltage clamp can be established. However, if one wishes to
"clamp" a field, one must control the potential at many locations. Each
location requires a different flux and thus a different amplifier and
different electrodes to supply that flux. Without such complicated apparatus,
it is almost impossible to maintain a constant field in space [34]. Indeed, it
is nearly impossible to maintain any pre-specified field because it is
practically impossible to apply different fluxes at different locations. If
one assumes a constant field in a theory, without such apparatus in an
experiment, one is in effect introducing flux into the calculation and model
that is not present in the experimental setup. One is introducing an
artifactual flux likely to produce artifactual conclusions that are not
relevant to the original experiment. [69, 41].

.

By taking the time derivative of total energy function (10), we have

$\displaystyle\frac{dE_{tot}}{dt}$ $\displaystyle=$
$\displaystyle\int_{\Omega}\sum_{i=1}^{4}\left\\{\mu_{i}\frac{d[C_{i}]}{dt}\right\\}dx+\int_{\Omega}\bm{E}\cdot\frac{d\bm{D}}{dt}dx$
(19) $\displaystyle=$
$\displaystyle\int_{\Omega}\sum_{i=1}^{4}\left\\{\mu_{i}\frac{d[C_{i}]}{dt}\right\\}dx-\int_{\Omega}\nabla\phi\cdot\frac{d\bm{D}}{dt}dx$
$\displaystyle=$
$\displaystyle\int_{\Omega}\sum_{i=1}^{4}\left\\{\mu_{i}\frac{d[C_{i}]}{dt}\right\\}dx+\int_{\Omega}\phi\nabla\cdot\left(\frac{d\bm{D}}{dt}\right)dx$
$\displaystyle=$
$\displaystyle\int_{\Omega}\sum_{i=1}^{4}\left\\{\mu_{i}\frac{d[C_{i}]}{dt}\right\\}dx+\int_{\Omega}\phi
F\sum_{i=1}^{4}\left\\{z_{i}\frac{d[C_{i}]}{dt}\right\\}dx$ $\displaystyle=$
$\displaystyle\int_{\Omega}\sum_{i=1}^{4}\left\\{\tilde{\mu}_{i}\frac{d[C_{i}]}{dt}\right\\}dx$
$\displaystyle=$
$\displaystyle-\int_{\Omega}\sum_{i=1}^{4}\left\\{\tilde{\mu}_{i}\nabla\cdot\bm{j}_{i}\right\\}dx-\int_{\Omega}\tilde{\mu}_{1}\nabla\cdot\bm{j}_{p}dx-\int_{\Omega}\mathcal{R}(\sum_{i=1}^{3}\alpha_{i}\tilde{\mu}_{i}-\alpha_{4}\tilde{\mu}_{4})dx$
$\displaystyle=$
$\displaystyle\int_{\Omega}\sum_{i=1}^{4}\left\\{\nabla\tilde{\mu}_{i}\cdot\bm{j}_{i}\right\\}dx+\int_{\Omega}\nabla\tilde{\mu}_{1}\cdot\bm{j}_{p}dx-\int_{\Omega}\mathcal{R}(\sum_{i=1}^{3}\alpha_{i}\mu_{i}-\alpha_{4}\mu_{4})dx+\int_{\partial\Omega}\sum_{i=1}^{4}\mu_{i}\bm{j}_{i}\cdot\bm{n}dS,$

where $\mu_{i}=RT\ln\frac{[C_{i}]}{c_{0}}+U_{i}$ and Eq. (2) is used.

By comparing with the dissipation functional, we have

$\displaystyle\bm{j}_{i}=-\frac{D_{i}}{RT}[C_{i}]\nabla\tilde{\mu}_{i},~{}i=1,2,3,$
(20a) $\displaystyle
RT\ln\left(\frac{\mathcal{R}}{k_{r}\left(\frac{[C_{4}]}{c_{0}}\right)^{\alpha_{4}}}+1\right)=\sum_{i=1}^{3}\alpha_{i}\mu_{i}-\alpha_{4}\mu_{4}.$
(20b)

And the corresponding energy influx rate is

$J_{E}=\sum_{i=1}^{4}\int_{\partial\Omega}\tilde{\mu}_{i}j_{i,extra}.$ (21)

For the pump flux, if we assume flux is only along the $z$ direction, then,

$\bm{j}_{p}=(0,0,\frac{f_{p}}{\partial_{z}\mu_{1}}).$ (22)

At equilibrium, we have

$\left\\{\begin{array}[]{l}\bm{j}_{i}=\nabla[C_{i}]_{eq}+\frac{z_{i}F}{RT}[C_{i}]_{eq}\nabla\phi_{eq}=0,\\\
\sum_{i=1}^{3}\alpha_{i}\mu_{i}([C_{i}]_{eq})-\alpha_{4}\
\mu_{4}([C_{4}]_{eq})=0,\end{array}\right.$

The last equation means

$\displaystyle
0=RT\ln\left(\frac{\Pi^{3}_{i=1}(\frac{[C_{i}]_{eq}}{c_{0}})^{\alpha_{i}}}{(\frac{[C_{4}]_{eq}}{c_{0}})^{\alpha_{4}}}\right)+(\sum_{i=1}^{3}\alpha_{i}U_{i}-\alpha_{4}U_{4}),$
(23)

According to the definition of equilibrium constant $k_{eq}$,

$k_{eq}=\frac{\Pi^{3}_{i=1}(\frac{[C_{i}]_{eq}}{c_{0}})^{\alpha_{i}}}{(\frac{[C_{4}]_{eq}}{c_{0}})^{\alpha_{4}}}$
(24)

Eq. (23) yields

$\ln{k_{eq}}=e^{-\frac{\Delta U}{RT}},$ (25)

with $\Delta U=\sum_{i=1}^{3}\alpha_{i}U_{i}-\alpha_{4}U_{4}$ .

Then combining Eqs. (20b) and (25) yields

$\displaystyle\ln\left(\frac{\mathcal{R}}{k_{r}\left(\frac{[C_{4}]}{c_{0}}\right)^{\alpha_{4}}}+1\right)$
$\displaystyle=$
$\displaystyle\ln\left(\frac{\Pi_{i=1}\left(\frac{[C_{i}]}{c_{0}}\right)^{\alpha_{i}}}{\left(\frac{[C_{4}]}{c_{0}}\right)^{\alpha_{4}}}\right),$
(26)

which implies

$\mathcal{R}=k_{f}\left(\frac{[C_{1}]}{c_{0}}\right)^{\alpha_{1}}\left(\frac{[C_{2}]}{c_{0}}\right)^{\alpha_{2}}\left(\frac{[C_{3}]}{c_{0}}\right)^{\alpha_{3}}-k_{r}\left(\frac{[C_{4}]}{c_{0}}\right)^{\alpha_{4}},$

where $k_{f}=\frac{k_{r}}{k_{eq}}$ [104].

######  Remark 2.4.

Here $k_{eq}$ is dimensionless. $k_{r}$ and $k_{f}$ are with unit $s^{-1}$
[ozcan2022equilibrium].

Then the whole system is as follows

$\left\\{\begin{array}[]{l}\frac{d[C_{1}]}{dt}=\nabla\cdot(D_{1}\nabla[C_{1}]+D_{1}\frac{z_{1}F}{RT}[C_{1}]\nabla\phi)-\partial_{z}j_{p}-\alpha_{1}\mathcal{R},\\\
\frac{d[C_{2}]}{dt}=\nabla\cdot(D_{2}\nabla[C_{2}]+D_{2}\frac{z_{2}F}{RT}[C_{2}]\nabla\phi)-\alpha_{2}\mathcal{R},\\\
\frac{d[C_{3}]}{dt}=\nabla\cdot(D_{3}\nabla[C_{3}]+D_{3}\frac{z_{3}F}{RT}[C_{3}]\nabla\phi)-\alpha_{3}\mathcal{R},\\\
\frac{d[C_{4}]}{dt}=\nabla\cdot(D_{4}\nabla[C_{4}]+D_{4}\frac{z_{4}F}{RT}[C_{4}]\nabla\phi)+\alpha_{4}\mathcal{R},\\\
-\nabla\cdot(\varepsilon_{0}\varepsilon_{r}\nabla\phi)=\sum_{i=1}^{4}z_{i}F[C_{i}],\end{array}\right.$
(27)

with

$\displaystyle\mathcal{R}=k_{f}\left(\frac{[C_{1}]}{c_{0}}\right)^{\alpha_{1}}\left(\frac{[C_{2}]}{c_{0}}\right)^{\alpha_{2}}\left(\frac{[C_{3}]}{c_{0}}\right)^{\alpha_{3}}-k_{r}\left(\frac{[C_{4}]}{c_{0}}\right)^{\alpha_{4}}.$
(28)

and boundary conditions

$\left\\{\begin{array}[]{ll}\bm{j}_{i}\cdot\bm{n}=j_{extra},i=1\cdots
4,&\mbox{on~{}}\partial\Omega,\\\
\bm{D}\cdot\bm{n}=0,&\mbox{on~{}}\partial\Omega.\end{array}\right.$ (29)

######  Remark 2.5.

If we assume that one of the reactants is an electron, for instance $C_{3}$,
and supplied by an thin electrode along z direction, the density of electron
$[C_{3}]=\rho_{e}=\rho(z)\delta(x_{0},y_{0})$. Then the model is changed to

$\left\\{\begin{array}[]{l}\frac{d[C_{1}]}{dt}=\nabla\cdot(D_{1}\nabla[C_{1}]+D_{1}\frac{z_{1}e}{RT}[C_{1}]\nabla\phi)-\partial_{z}(j_{p})-\alpha_{1}\mathcal{R}\delta(x_{0},y_{0}),\\\
\frac{d[C_{2}]}{dt}=\nabla\cdot(D_{2}\nabla[C_{2}]+D_{2}\frac{z_{2}e}{RT}[C_{2}]\nabla\phi)-\alpha_{2}\mathcal{R}\delta(x_{0},y_{0}),\\\
\frac{d[C_{4}]}{dt}=\nabla\cdot(D_{4}\nabla[C_{4}]+D_{4}\frac{z_{4}e}{RT}[C_{4}]\nabla\phi)+\alpha_{4}\mathcal{R}\delta(x_{0},y_{0}),\\\
\frac{d\rho(z)}{dt}=-\partial_{z}\bm{j}_{e}-\alpha_{3}\mathcal{R}\delta(x_{0},y_{0}),\\\
-\nabla\cdot(\varepsilon_{0}\varepsilon_{r}\nabla\phi)=\sum_{i=1,2,4}z_{i}e[C_{i}]-F\rho(z)\delta(x_{0},y_{0}).\end{array}\right.$
(30)

######  Remark 2.6.

When the reaction and ions are in an electrolyte, the fluid effect needs to be
taken into consideration. In this case, the energy functional is changed to be

$\displaystyle E_{tot}$ $\displaystyle=$ $\displaystyle
E_{kin}+E_{ent}+E_{int}+E_{ele}$ (31) $\displaystyle=$
$\displaystyle\int_{\Omega}\frac{\rho|\mathbf{u}|^{2}}{2}dx+\sum_{i=1}^{4}\int_{\Omega}RT\left\\{[C_{i}]\left(\ln{\left(\frac{[C_{i}]}{c_{0}}\right)}-1\right)\right\\}dx+\int_{\Omega}\sum_{i=1}^{4}[C_{i}]U_{i}dx$

and the dissipation functional is changed to

$\displaystyle\Delta=\int_{\Omega}2\eta|\bm{D}_{\eta}|^{2}dx+\int_{\Omega}\left\\{\sum_{j=1}^{4}|\bm{j}_{i}|^{2}+RT\mathcal{R}\ln\left(\frac{\mathcal{R}}{k_{r}\left(\frac{[C_{4}]}{c_{0}}\right)^{\alpha_{4}}}+1\right)\right\\}dx-\int_{\Omega}f_{p}dx,$
(32)

where $\bm{D}_{\eta}=\frac{\nabla\bm{u}+(\nabla\bm{u})^{T}}{2}$ and $\bm{u}$
is the velocity. We can use the Energy variation method to get the diffusion-
reaction-convection model as follows

$\left\\{\begin{array}[]{l}\frac{\partial[C_{1}]}{\partial
t}+\nabla\cdot([C_{1}]\bm{u})=\nabla\cdot(D_{1}\nabla[C_{1}]+D_{1}\frac{z_{1}F}{RT}[C_{1}]\nabla\phi)-\partial_{z}j_{p}-\alpha_{1}\mathcal{R},\\\
\frac{\partial[C_{2}]}{\partial
t}+\nabla\cdot([C_{2}]\bm{u})=\nabla\cdot(D_{2}\nabla[C_{2}]+D_{2}\frac{z_{2}F}{RT}[C_{2}]\nabla\phi)-\alpha_{2}\mathcal{R},\\\
\frac{\partial[C_{3}]}{\partial
t}+\nabla\cdot([C_{3}]\bm{u})=\nabla\cdot(D_{3}\nabla[C_{3}]+D_{3}\frac{z_{3}F}{RT}[C_{3}]\nabla\phi)-\alpha_{3}\mathcal{R},\\\
\frac{\partial[C_{4}]}{\partial
t}+\nabla\cdot([C_{4}]\bm{u})=\nabla\cdot(D_{4}\nabla[C_{4}]+D_{4}\frac{z_{4}F}{RT}[C_{4}]\nabla\phi)+\alpha_{4}\mathcal{R},\\\
-\nabla\cdot(\varepsilon_{0}\varepsilon_{r}\nabla\phi)=\sum_{i=1}^{4}z_{i}F[C_{i}],\\\
\rho(\frac{\partial\bm{u}}{\partial t}+(\bm{u}\nabla)\cdot\bm{u}+\nabla
p=\nabla(\eta(\nabla\bm{u}+(\nabla\bm{u})^{T}))-(\sum_{i=1}^{4}z_{i}F[C_{i}])\nabla\phi\\\
\nabla\cdot\bm{u}=0.\end{array}\right.$ (33)

Note we are not here considering transient problems in which charge is stored
in polarization fields. These will be studied separately so we can deal with
the important experiments reported in [56, 51, 55, 72, 87, 63]. Transient
problems are obviously important if reactions are studied on the atomic scale
of distance and time (angstroms and femtoseconds) because the polarization
currents are large. Dealing with those currents requires use of a universal
form of the Maxwell equations combined with an appropriate model of the stress
strain relation of charge in a viscoelastic structure, commonly called
polarization. Speaking loosely, the transient problems can be dealt with in
circuits by a generalization of Kirchhoff’s law [17, 91] to describe the
actual transient currents that flow through an ideal resistor [88].

It is important to realize that currents (and fluxes) cannot be computed by
methods that assume the currents and fluxes are zero. Electrostatics cannot
compute currents because currents and fluxes involve time and electrostatics
does not [62]. Electrostatics does not include Maxwell’s Ampere law that is
the universal coupler of current to electric and magentic fields. In the
context of cytochrome $c$ oxidase these issues come to the fore. Models
without electron or proton current as variables do not describe the ’transfer
function’ of the transporter being studied. Models cannot calculate Ohm’s law
(for system with large and small currents and electrical potentials) if the
models assume currents are zero.

In fact, using a formulation of electrodynamics that explicitly involves
current is straightforward, as engineers have known for a very long time,
going back to Heaviside [11, 27] and are worked out in practical detail in
[17]. Kirchhoff’s current law allows analysis of systems of great importance,
without dealing with charges explicitly. That is why analysis of electronic
circuits does not need to use distributions of charges but rather uses
Kirchhoff’s current law or its generalization, conservation of total current.
Kirchhoff’s current law is an exact corollary of the Maxwell equations
themselves, if current includes the displacement current [17, 88, 91]. It
might seem that another corollary of the Maxwell equations, the continuity
equation, can be used instead of Kirchhoff’s law for total current. And it is
certainly true that the continuity equation of electrodynamics contains the
same information as conservation of (total) current, all conjoined with the
Maxwell equations. But that information is not useful when enormous numbers of
charges are involved as in cytochrome $c$ oxidase or in other macroscopic
scale systems like the electronic circuits of our computers. The information
implicit in the flux of charges is only usable when written as total current
that is conserved perfectly whenever the Maxwell equations are valid. This
formulation using conservation of total current does not require explicit
treatment of charges. The continuity equation does require the explicit
treatment of charges, and their significant interactions, whether involving
two charge interactions, three charge, …. or the interactions of an entire
cluster expansion. The significant interactions of charge are difficult to
understand or even enumerate and more difficult to compute [84, 67].
Kirchhoff’s current law is easy to understand and trivial to compute.

##  3 An Electro-osmotic Model of cytochrome $c$ oxidase

Here we propose a specific model of cytochrome $c$ oxidase (or Complex IV) as
an example so our approach can be seen in action. The schematic structure of
cytochrome $c$ is shown in Fig.1 a, where both two channels from mitochondria
matrix (inside), D and K, are taken into consideration. Here, $E$ denotes the
end of $D$ channel. And the end of K channel is asumped to be the binuclear
center (BNC) denoted by $B$ where the chemical reaction (36) occurs. The
protons accumulated in E are transported to the BNC and the proton loading
site (PLS), denoted by X. A pump is located between E and PLS. The pump
provides energy that comes from concentration gradients, namely gradients of
chemical potential at BNC. Then finally, the proton is pushed out from PLS to
the inter membrane space, outside the mitochondrion.

It is clear that this model is incomplete at best, and in some sense wrong, at
worst. We depend on our experimental colleagues to help us correct and improve
the model, for example, by including mechanisms we have over simplified.
Enormous detail of the chemical reactions is described in the literature, with
more intermediates being reported frequently. We do not include these
intermediates.

Let $\mathcal{\rho}_{e}=\rho_{0}\delta(x_{0},y_{0},z_{0},t)$. Integrating the
diffusion-reaction equation

$\frac{d[C]_{i}}{dt}=-\nabla\cdot\bm{j}_{i}-\alpha_{i}\mathcal{R},$ (34)

in the complex IV compartment yields

$\eta\frac{d\bar{C}_{i}}{dt}=J_{i}^{in}-J_{i}^{out}-\alpha_{i}\mathcal{R},$
(35)

where $\eta_{mat}$, $\eta_{ims}$ and $\eta$ are the volumes of mitochondrial
compartment, inter membrane space and reaction compartment, respectively.

The chemical reaction in the cytochrome $c$ oxidase Complex IV is

$2H^{+}+\frac{1}{2}O_{2}+2e^{-}\overset{k_{f}}{\underset{k_{r}}{\rightleftharpoons}}H_{2}O.$
(36)

Figure 1: Schematic. (a) Complex IV structure; (b) Circuit diagram of model
(37) - (38); (c) Circuit diagram when the rectifier is located between the
protein loading site PLS and P side. The drawing of (a) is taken from [68]
according to the policies of the Proceedings of the National Academy of
Sciences https://www.pnas.org/author-center/publication-charges#author-rights-
and-permissions examined on June 18, 2022. We thank the authors for providing
all of us such a helpful figure.

For simplicity, we follow the Hodgkin Huxley tradition and fix the proton
concentrations in the mitochondrial matrix (inside) and the inter membrane
space outside the mitochondria so they do not vary with time or flow. More
general treatments in which concentrations are changed with time by flow are
possible as have been done in even more complex structures, Such analysis has
been done in a bi-domain model of the lens of the eye and a tri-domain model
of the optic nerve and glia [95, 105, 105, 110, 111].

Here for simplicity, we assume that the concentration of oxygen at the B site
is a constant. If the oxygen varies with time, an additional equation can be
used to describe the dynamics of oxygen. The properties of that term will be
determined either directly by experimentation or by a higher resolution model
as in [76]. We do not expect the extra term to introduce significant
mathematical, numerical, or computational difficulties.

Many variables are needed to keep track of all the potentials and
concentrations in the variouis regions of our model. The concentrations and
potentials at E242; at BNC; and at the proton loading site (PLS), are
different inside and outside the mitochondria, as is the electron
concentration (see Fig.(1) (a)). They are described by the variables
$[H]_{E}$, $[H]_{B}$, $[H]_{x}$,$\phi_{E}$, $\phi_{B}$, $\phi_{X}$,
$\phi_{N}$, $\phi_{P}$ and $\rho_{e}$, respectively.

$\displaystyle\frac{d[H]_{E}}{dt}=\frac{S_{v}}{F}(I_{N2E}-I_{E2X}-I_{E2B}),$
(37a)
$\displaystyle\frac{d[H]_{B}}{dt}=\frac{S_{v}}{F}(I_{E2B}+I_{N2B})-2\mathcal{R},$
(37b) $\displaystyle\frac{d[H]_{X}}{dt}=\frac{S_{v}}{F}(I_{E2X}-I_{X2P}),$
(37c) $\displaystyle\frac{d\rho_{e}}{dt}=\frac{-S_{v}}{F}I_{e}-2\mathcal{R},$
(37d) $\displaystyle
C_{E}\frac{d(\phi_{E}-\phi_{N})}{dt}=(I_{N2E}-I_{E2X}-I_{E2B}),$ (37e)
$\displaystyle C_{B}\frac{d(\phi_{B}-\phi_{N})}{dt}=I_{E2B}+I_{N2B}+I_{e},$
(37f) $\displaystyle C_{X}\frac{d(\phi_{X}-\phi_{P})}{dt}=(I_{E2X}-I_{X2P}),$
(37g) $\displaystyle
C_{m}\frac{d(\phi_{N}-\phi_{P})}{dt}+I_{leak}+I_{X2p}+I_{e}=0,$ (37h)

with currents

$\displaystyle
I_{N2B}=g_{K}(\phi_{N}-\phi_{B}-\frac{RT}{F}\ln\frac{[H]_{B}}{[H]_{N}})=\frac{g_{K}}{F}(\mu_{N}-\mu_{B}),$
(38a) $\displaystyle
I_{E2B}=g_{B}(\phi_{E}-\phi_{B}-\frac{RT}{F}\ln\frac{[H]_{B}}{[H]_{E}})=\frac{g_{B}}{F}(\mu_{E}-\mu_{B}),$
(38b) $\displaystyle I_{E2X}=I_{pump}+I_{xleak},$ (38c) $\displaystyle
I_{leak}=g_{m}(\phi_{N}-\phi_{P}-E_{other}),$ (38d) $\displaystyle
I_{e}=-FJ_{e},$ (38e) $\displaystyle I_{xleak}=-g_{E}(\mu_{X}-\mu_{E}),$ (38f)
$\displaystyle
I_{pump}=\left\\{\begin{array}[]{cc}g_{pump}max(R_{c},0)(\mu_{X}-\mu_{E}),\mu_{X}-\mu_{E}<\delta_{th},\\\
g_{pump}max(R_{c},0)\delta_{th}\exp{\left(-\frac{(\mu_{X}-\mu_{E})}{\varepsilon}\right)},\mu_{X}-\mu_{E}\geq\delta_{th},\end{array}\right.$
(38i)
$\displaystyle\mathcal{R}=k_{f}[H^{+}]^{2}[O_{2}]^{1/2}\rho_{e}^{2}-k_{r}[H_{2}O].$
(38j)

We follow the review of Wikström [89] and implement switching functions
without invoking conformation changes of the distribution of mass. We treat
cytochrome c oxidase as a Coulombic system and use rectifiers to implement the
switching functions that provide alternating access of an occluded state. Here
we discuss two cases. In one case, a rectifier between $N$ and $E$ blocks the
proton flows. In the other case, a rectifier between $X$ and $P$ blocks the
backward proton flows. Then, the currents $I_{N2E}$ and $I_{X2P}$ are modelled
in the follows two cases.

  * •

Case 1: the rectifier is between $N$ and $E$ as shown in Fig. 1b

$\displaystyle
I_{N2E}=max\left(g_{D}\left(\phi_{N}-\phi_{E}-\frac{RT}{F}\ln\frac{[H]_{E}}{[H]_{D}}\right),-SW_{0}\right)=max\left(\frac{g_{D}}{F}(\mu_{N}-\mu_{E}),-SW_{0}\right),$
(39a) $\displaystyle
I_{X2P}=g_{X}(\phi_{X}-\phi_{P}-\frac{RT}{F}\ln\frac{[H]_{P}}{[H]_{X}})=\frac{g_{X}}{F}(\mu_{X}-\mu_{P}),$
(39b)

  * •

Case 2: the rectifier is between $X$ and outside as shown in Fig. 1c

$\displaystyle
I_{N2E}=g_{D}\left(\phi_{N}-\phi_{E}-\frac{RT}{F}\ln\frac{[H]_{E}}{[H]_{D}}\right),$
(40a) $\displaystyle
I_{X2P}=max\left(g_{X}(\phi_{X}-\phi_{P}-\frac{RT}{F}\ln\frac{[H]_{P}}{[H]_{X}}),-SW_{0}\right)=max\left(\frac{g_{X}}{F}(\mu_{X}-\mu_{P}),-SW_{0}\right),$
(40b)

where $SW_{0}$ is the threshold for the turn-off of the rectifier and
$SW_{0}=0$ stands for perfect rectifier. We reiterate that pn junctions are
used to rectify the movement of the pseudo-ions holes and electrons throughout
our digital circuitry. Analogous distributions of permanent charge provided by
acid and base side chains of proteins produce rectification of charge movement
in ionic systems.

We use Kirchhoff’s law and the conductance formulation of Hodgkin and Huxley.
Complex properties are hidden by a nonlinear time dependent version of Ohm’s
law and modelled by a conductance, as did Hodgkin and Huxley [32, 8, 14, 47,
6, 7]. Alternating access (with its implied occluded state) is described by a
switching function for the D channel using equation (39a). This is a classical
rectifier function and (when $SW_{0}=0$), allows current only to flow from D
to E: no backward flow is allowed.

Many properties of the model depend on the pump current between E242 and the
Proton Loading Site PLS. We assume that in the ordinary situation the pump
strength depends on the reaction rate and the chemical reaction difference
between two sites. However, in the less ordinary situation, when the
difference is too large, the pump may not be able to overcome the barrier. A
turn-off threshold is assigned to the pump for that reason. We assume that
when the difference in chemical potential $\mu_{x}-\mu_{E}$ is greater than
the threshold, the pump current decreases exponentially to zero as it turns
off. More realistic, and complex properties of the pump will undoubtedly be
needed to explain some functions of cytochrome c oxidase. They can easily be
incorporated into our model, as these properties are measured and modelled.

Our ‘electro-osmotic’ model is a ‘master equation’ approach building on the
work of Hummer and Kim, [74, 58, 65] but showing how to exploit conservation
of current in the form of Kirchhoff’s current law. This approach is used
throughout electrical and electronic engineering to design semiconductor
devices, as textbooks document ($op.cit.$) perhaps most eloquently in the
modern automated circuit design literature built on Kirchhoff’s law [73].
Currents are sufficient for such automated design. Charges are not needed
except in switched-capacitor networks (p. 64 of [73]) .

We extend the classical use of Kirchhoff’s law that forms the foundation of
circuit design to include chemical reactions. We must include chemical
reactions to drive currents of electrons, protons and other ions because that
is how cytochrome $c$ oxidase functions. The essential function of cytochrome
$c$ oxidase is to convert a flow of electrons to a flow of protons from inside
the mitochondrion to outside it. The electrons that are inputs to the
cytochrome $c$ oxidase are presented to the enzyme attached to the heme group
of cytochrome $c$ itself.

The existing literature analyzes these systems without explicitly dealing with
currents, making the task either impossible (by using a theory that assumes a
zero value for the fluxes being studied) or very difficult (by involving a
staggering number of charges). Using currents instead of charges avoids these
difficulties and has the added advantage of automatically satisfying Maxwell’s
equations, if total current is used in Kirchhoff’s law.

This approach is incomplete because it does not deal with all the charge in
the circuit formed by cytochrome $c$ oxidase But those details of charge are
not needed in the design of electronic circuits. That simple fact can be
verified by examining textbooks of circuit design (as already cited).

In circuit analysis of this type, some questions about charges need not be
asked: for example, the atomic mechanism of current flow (particularly
electron flow) can be ignored. The function of the circuit is independent of
the details of the components of current in wires, for example, with only a
few exceptions [17]. Thus, we de-emphasize the atomic details of the various
pathways that provide electron flow to the main reaction centers. For us,
these pathways are wires. The atomic and chemical details of electron flow in
these wires are known in breathtaking detail and we are sorry that we do not
seem to need to use these magnificent results, but, at this resolution, we do
not.

A key biological result is that some of the coupling so important to
understand the electro-osmotic properties of a mitochondrion depends on the
macroscopic conservation of current, i.e., Kirchhoff’s law applied to the
entire mitochondrion. The application of Kirchhoff’s law to mitochondrial
transport, and active transport in general, is not common in the literature.
But Kirchhoff’s law has been used in another branch of biophysics for a long
time, for some eighty years. Kirchhoff’s law is the keystone of the analysis
of ion channels. Kirchhoff’s law is the keystone that supports the structure
of the Hodgkin approach to the action potential by balancing the various
components of current, summing them to zero in the appropriate (finite)
geometries, like those of mitochondria, the way the keystone of an arch sums
mechanical forces.

Conservation of current provides the coupling in other biophysical
applications, e.g., generation and propagation of the action potential,
linking atomic scale properties of ion channels of one type to properties
(e.g., opening) of another type. In a classical action potential, the opening
of sodium channels is coupled to the opening of other sodium channels, and to
the closing of potassium channels by the electric field, not by anything else.
There is no steric or chemical interaction between the channels. The coupling
is essential to the function of the nerve cell, but that coupling is described
by a version of the Maxwell equations (called the cable or telegrapher’s
equation) not by equations of chemical kinetics. The ion channels of the
action potential act independently in the chemical sense because they are so
far apart, without opportunity for short range or chemical interactions. The
ion channels are not independent, in the physiological or physical sense,
however. Rather they are coupled by the electric field. The electrical field
is that which satisfies the Maxwell equations, or their equivalent,
Kirchhoff’s current law.

What we propose here is in the tradition of Hodgkin’s treatment of ionic
channels, but we include the chemical reactions that are the essence of
oxidative phosphorylation and the life and function of mitochondria. We are
more than aware that a detailed analysis of alternating access, the occluded
states, and the switching function is needed to understand cytochrome $c$
oxidase. That analysis needs the currents flowing to be analyzed along with
the atomic detail of the water-gate switch [50, 89], in our view. The switches
act on currents and of course satisfy conservation of current. An analysis of
charges cannot easily guarantee conservation of current, and classical
chemical analysis precludes large currents because assuming equilibrium or
near equilibrium conditions is clearly inappropriate for a system like
Reaction Center IV designed for the efficient handling of large flows of
electrons and protons. Here we describe alternating access with the classical
equation of a rectifier to highlight the possibility that occluded states and
alternating access are the biochemical names for what engineers call
rectification. It is important to realize that rectification is an automatic
unavoidable consequence of the distribution of doping in semiconductors, for
example in the classical PN diode. This rectification occurs with no change in
the spatial distribution of mass (i.e., with no change in what is usually
called conformation) and so it is compatible with the view cited above that
cytochrome c oxidase functions without changes in the spatial distribution of
mass, i.e., without what is classically called conformation change. In the
rectification mechanism, the switching (rectification) occurs because of a
radical change in the distribution of electrical potential, which in turn
allows current flow in one direction and not another. The distribution of
potential depends on the distribution of mobile electrons which have almost no
mass. The conformation of the potential profile and thus the electric field
creates a barrier for current flow in one direction but not in another,
because of the effects of doping (permanent charge) and mobile charge combined
in the Maxwell Gauss law, or the Poisson equation. This system is rather
complex, although completely understood and used in literally billions of
places in each of our computers, The system involves diffusive and electrical
movement of electrons (and holes) driven by the gradients of chemical
potential (e.g., concentration) and electrical potential. As the electrical
potential changes sign, diffusive and electrical flow changes. As
concentrations change, diffusive and electrical flow change in other ways. All
interact through the changing fields of electrical and chemical potential. A
few pages (not just a few words or sentences) are needed to explain how each
kind of movement (diffusive, electrical, holes and electrons) contribute to
rectification. See textbooks of semiconductor circuit design, e.g, [53, 39].
It is also important to consult research articles [44, 45] to understand the
oversimplifications of the textbook discussions and to validate them. More
elaborate patterns of doping, starting with the PNPN designs (thyristors,
Silicon Controlled Rectifiers = SCR) are used in power transistors. Analogous
spatial distributions of permanent charge (and acid and base side chains)
might be used to implement switches in cytochrome c oxidase.

We note rectification arising from the distribution of permanent charge in a
protein was proposed by Mauro a very long time ago [13, 12]. Such rectifiers
of ionic current have even more complex properties than semiconductor
rectifiers because concentrations of current carriers in biological solutions
can be changed independently of electrical potential, which is not often the
case in analogous semiconductor systems. Ionic rectifiers were built a long
time ago using a biological protein as a template [60] and are now used
routinely in the ionic channels of nanotechnology [108] even in a hybrid chip
that can enable a scalable integrated ionic circuit platform for micro-total-
analytical systems [82].

The switch of Reaction Complex IV is likely to involve both the distribution
of permanent charge (mostly acid and base side chains), and chemical
interactions as described in the water-gate model [50, 89], perhaps also
involving spatial distribution of dielectric properties as well [106]. It
seems premature to attack this problem here, as important as it is for the
function of cytochrome $c$ oxidase, and all alternating access transporters,
for that matter. Here, we simply describe the rectification without further
analysis of how it arises from the distribution of permanent charge and other
properties of the transporter structure.

Of course, other possibilities exist. Alternating access might arise, for
example, from bubbles in the conduction pathway, as we are studying in other
work [114]. Two bubbles might act as coupled activation and inactivation
gates, correlated to provide alternating access to an occluded state, for
example.

##  4 Computational Studies

In this section, we carry out several computational studies to explore the
effects of various conditions on proton transport efficiency. The initial
values and default parameters are listed in Table 1-2.

Variable | Notations | Values (with Unit)
---|---|---
$E_{242}$ site $H^{+}$ concentration | $[H]_{E}$ | 0.01196 $\mu M$
BNC site $H^{+}$ concentration | $[H]_{B}$ | 0.01682 $\mu M$
PLS site $H^{+}$ concentration | $[H]_{X}$ | 0.01441 $\mu M$
BNC site electric density | $\rho_{e}$ | 0.01166 $\mu M$
$E_{242}$ site electric potential | $\phi_{E}$ | -5 $mV$
BNC site electric potential | $\phi_{B}$ | -14.1562 $mv$
PLS site electric potential | $\phi_{X}$ | 200 $mv$
N site electric potential | $\phi_{N}$ | 0 $mv$
P site electric potential | $\phi_{P}$ | 160 $mv$
Table 1: Default Initial Values Variable | Notations | Values (with Unit)
---|---|---
$E_{242}$ site effective capacitance | $C_{D}$ | 1E-1 $fAms/mV/(\mu m)^{2}$
BNC site effective capacitance | $C_{B}$ | 1E-1 $fAms/mV/(\mu m)^{2}$
PLS site effective capacitance | $C_{X}$ | 1E-1 $fAms/mV/(\mu m)^{2}$
Membrane capacitance | $C_{X}$ | 7.5E-2 $fAms/mV/(\mu m)^{2}$
D channel conductance for $H^{+}$ | $g_{D}$ | 3.75E-3$pS/(mum)^{2}$
K channel conductance for $H^{+}$ | $g_{K}$ | 1E-3 $pS/(\mu m)^{2}$
E2B channel conductance for $H^{+}$ | $g_{B}$ | 5E-2 $pS/(\mu m)^{2}$
E2X channel conductance for $H^{+}$ | $g_{E}$ | 1E-3 $pS/(\mu m)^{2}$
E2X Pump rate for $H^{+}$ | $g_{P}$ | 369 $pSms/(\mu m)^{2}\mu M$
X2P channel conductance for $H^{+}$ | $g_{X}$ | 9.8E-4 $pS/(\mu m)^{2}$
Membrane conductance for leak | $g_{m}$ | 1E-5$pS/(\mu m)^{2}$
Mito. matrix $H^{+}$ concentration | $[H]_{mat}$ | 0.01 $\mu M$
Mito. inner membrane space $H^{+}$ concentration | $[H]_{ims}$ | 0.063 $\mu M$
Nernst Potential due to other Ions | $E_{Other}$ | $-160mV$
Reaction site $[O_{2}]$ concentration | $[O_{2}]$ | 0.0028 $\mu M$
Reaction site $[H_{2}O]$ concentration | $[H_{2}O]$ | 0 $\mu M$
Electron current | $I_{e}$ | -5.24 $fA$
Forward reaction rate coefficient | $k_{f}$ | 1333
Backward reaction rate coefficient | $k_{r}$ | 0.005
surface volume ratio | $S_{v}$ | 1000
Potential Threshold | $\delta_{th}$ | 210 mv
Decay rate | $\varepsilon$ | 1 $(ms)^{-1}$
Table 2: Parameters

###  4.1 Effect of electron current: Input to Output Relations

We first check Figs.2 the effect of electron current $I_{e}$ on the efficiency
of Complex IV. The case 1 (inside to E rectifier) results are represented by
blue circle and the case 2 (X to outside rectifier) results are represented by
red square.

The ratios between currents and the supplied electron current are measures of
the transfer function or ’gain’ of cytochrome c oxidase. According to the
previous study [89], the ratios are $\frac{I_{X2P}}{I_{e}}=-1$,
$\frac{I_{E2X}}{I_{e}}=-1$ and $\frac{I_{E2P}}{I_{N2E}+I_{N2B}}=-0.5$ at the
normal state. These ratios mean that (nominally) each input electron will
bring 2 protons from the N side. One of the protons is used for the chemical
reaction and the other one is pumped to the P side becoming an output in that
way.

Fig. 2 (d)-(h) confirm that when the electron supply is sufficient
($|I_{e}|\geq 5.24fA$) these ratios can be maintained. However, if the input
electron current decreases (in magnitude), the reaction rate decreases
linearly (see Fig. 2 (a)) since $\mathcal{R}=\frac{-S_{v}}{2F}I_{e}$ at
equilibrium according to Eq. (37d). The pump strength depends on reaction rate
(Eq.(38i)), so the pump current $I_{Pump}$ decreases hand in hand with
reaction rate.

Beyond a threshold, the total current between E242 and PLS $I_{E2x}$ becomes
negative, $provided$ the rectifier is located $between$ the inside and E site
(blue lines with circles). The the protons leak back form PLS to E242. This
‘leak back’ can be seen in Fig. 2 (d), where the positive ratio means means
that the $I_{E2X}$ is negative (because $I_{e}$ is negative, with our sign
conventions). Proton back flow from outside to the proton loading site PLS
occurs in this case. The accumulated protons in E242 increase the chemical
potential $\mu_{E}$ to be greater than $\mu_{N}$ and $\mu_{E}$, which leads to
more current from E242 to the reaction cite B (see Fig. 2 (f)) and back flow
from reaction cite to N side (see Fig.2 (h)). The rectifier blocks the direct
back flow from E242 to the N side. The urrent $I_{N2E}$ becomes zeros (see
Fig. 2 (g)). In this case, the proton flow pattern is shown in Fig. 3 (b).

The location of the rectifier is important. Behavior is different when the
rectifier is moved. When the rectifier is between the protein loading site PLS
and outside (red lines with squares), the backward flow from outside to PLS is
blocked to be zero. Then the current $I_{E2X}$ is also zero at the
equilibrium, according to Eq. (37c), which means $\mu_{E}=\mu_{X}$ as shown in
Fig. 2 (b). The protons are still transported from inside to E242 then to BNC.
$\frac{I_{N2E}}{I_{e}}=-0.68$ through D channel and directly to BNC with ratio
$\frac{I_{N2B}}{I_{e}}=-0.32$ through the K channel. Fig. 3(c) shows the
proton flow pattern in this situation.

We suspect that the rectifier between the protein loading site PLS and outside
is closer to the real biology setup, because it blocks backward flow. For that
reason, we mainly present the results with the rectifier between PLS and
outside. Our approach can of course handle almost any location or properties
of the rectifier/switch once they are specified by experiment or higher
resolution models.

Figure 2: Pump efficiency at equilibrium states with different electron
current and different threshold. (a) Reaction rate; (b) $\mu_{X}-\mu_{E}$; (c)
$I_{Pump}$;(d) $I_{E2X}/I_{e}$; (e) $I_{X2P}/I_{e}$; (f) $I_{E2B}/I_{e}$; (g)
$I_{N2E}/I_{e}$; (h) $I_{N2B}/I_{e}$. Red line: Switch between X and outside;
Blue dash line: Switch between inside and E.

Figure 3:  Schematic for the proton flow pattern. Equation numbers that define
arrows are shown, e.g., eq. (38b) for $I_{X2P}$. (a) Normal state; (b)
Backward flow state with Perfect switch between N and E242; (c) Flow state
with perfect switch between the proton Loading Site PLS and outside.

###  4.2 Effect of Proton Concentration

In this section, we study the effect of proton concentrations in the
intermembrane space (outside) by increasing from the default value $0.06\mu M$
to $0.15\mu M$.

Figs. 5 and 4 show the equilibrium states of concentrations and pump
efficiency at different proton concentrations with difference leak conductance
$g_{m}$.

First, Fig. 4 a. illustrates the reaction rate with different $[H]_{P}$ keeps
a constant since $\mathcal{R}=\frac{-S_{v}}{2F}I_{e}$ at equilibrium. When the
leak conductance is zero, the complex IV efficiency does not change, i.e.
$\frac{I_{X2P}}{I_{e}}=1$ and the flow pattern is shown in Fig.3a. When the
shunt conductance $g_{m}$ is larger than zero, the pump resistance increases
with the outside proton concentration. This produces a decrease of the complex
IV efficiency all the way to zero after the threshold. Then the proton pattern
is the same as in Fig.3c, where all the protons pumped from inside through D
and K channels are consumed by the reaction at BNC.

The concentrations of electrons and protons at $E$ and $B$ sites are small
perturbations in all cases. The concentration in the PLS is almost a constant
with different $[H]_{P}$ when the leak conductance is large. However, it
increases tremendously when the leak conductance is small. Large leak
conductance simulates voltage clamp conditions, which do not describe the
normal functional state of the mitochondrion. Small leak conductance
presumably corresponds to the natural state in which the sum of all currents
across the mitochondrion is ‘clamped’ to zero, by Kirchhoff’s current law,
because there is nowhere else for the current to flow.

Figure 4: Pump efficiency at equilibrium states with different proton
concentration and different leak conductance. (a) Reaction rate; (b)
$I_{N2E}/I_{e}$; (c) $I_{N2B}/I_{e}$; (d) $I_{E2X}/I_{e}$; (e)
$I_{E2B}/I_{e}$; (f) $I_{X2P}/I_{e}$. Black dash line: $g_{m}=0$; Red line
with circle: $g_{m}=10^{-5}$; Red line with square: $g_{m}=10^{-3}$.

Figure 5: Concentration at equilibrium states with different proton
concentration and different leak conductance. (a) Electron concentration
$\rho_{e}$; (b) $[H]_{E}$; (c) $[H]_{B}$;(d) $[H]_{X}$. Black dash line:
$g_{m}=0$; Red line with circle: $g_{m}=10^{-5}$; Red line with square:
$g_{m}=10^{-3}$.

###  4.3 Kirchhoff clamp

Most of this paper describes cytochrome c oxidase embedded in a mitochondrion
approximating a preparation without other members of the respiratory chain,
but with otherwise normal properties. The mitochondrion is a small cell, as it
were, in which the interior potential is unlikely to vary substantially with
macroscopic location, on the micron scale, because the cell is much smaller
than the length constant of cable theory. In such a system, Kirchhoff’s
current law ensures that the sum of all currents across the mitochondrial
membrane is zero. The currents are necessarily coupled by electrodynamics,
whether or not they are coupled by chemistry. If one current increases, the
sum of the others must decrease. A graph of one current against another will
give a definite ratio, a coupling ratio, if you will, when other variables are
held constant. The various currents are coupled by the electric potential. The
electrical potential contributes to the forces that drive the currents.
Chemical reactions may contribute as well. But even without chemical
reactions, coupling can occur through the electric field. While this may seem
strange in context of classical transport biophysics, it is well precedent ed
and understood in channel biophysics, that presumably follows the same laws of
physics. The coupling of sodium and potassium conductances that allow the
action potential to propagate are an example of coupling by the electric
field, without chemical interaction of the underlying protein molecules, as we
have discussed previously in this paper.

Coupling occurs because the electric field adopts the values that conserve
current. That is easy to prove from the Maxwell equations [92, 98]. In fact,
the conservation of current is a form of Kirchhoff’s law, so currents are
clamped to one another (i.e., coupled) in a “Kirchhoff clamp” if we want to
coin a phrase for what is really just the UNclamped, natural situation.

If the electrical potential is controlled, so it is not free to adopt the
value that conserves current, a different situation occurs altogether. This
situation is called a voltage clamp in electrophysiology, and was invented by
Cole and used by Hodgkin and Huxley to understand the mechanism of the action
potential. The voltage clamp loosens the Kirchhoff clamp because it has an
amplifier (that is outside the biological system) to supply current and
energy. Indeed, the Kirchhoff clamp of the natural mitochondrion is entirely
removed by the currents supplied by the voltage clamp amplifier.

In the voltage clamp, one current is not coupled to another current by the
voltage. They cannot be, because the voltage does not vary with the current.
The result is that coupling and flux ratios reflect chemical coupling, not
voltage coupling, in the voltage clamp setup. The result is that nearly every
experimental result is different in a voltage clamp and the natural UNclamped
situation.

The voltage clamp was invented to gain experimental control of currents so
they can be studied as Cole made abundantly clear, followed by Hodgkin and
Huxley. But the voltage clamp is not natural. It removes a natural form of
flux coupling. Flux coupling by the electric field is absent. Flux coupling by
the electric field is natural, just as natural in the mitochondrion as in the
nerve, just as natural in the generation of ATP as it is in the generation of
the nerve signal.

It is difficult to voltage clamp mitochondria, and preparations reconstituted
into bilayers (that can be voltage clamped) have other difficulties that
experimentalists often wish to avoid. Other methods are used to simulate a
voltage clamp, quite well, as it turns out.

In work on mitochondria, voltage clamp is usually produced indirectly, by
artificially increasing the leak conductance. An effective carrier of
potassium current like valinomycin is often added to solutions. When
valinomycin is present in large enough concentrations, it partitions into the
mitochnondrial membrane, and the leak conductance dominates. The potential
across the mitochondrial membrane is set by the equilibrium potential of the
leak conductance. If valinomycin is used to increase the leak conductance, the
potential is in fact close to the potassium equilibrium potential, independent
of current because valinomycin is remarkably selective for potassium ions.
Valinomycin clamps the potential to the potassium equilibrium potential.

Figs. 6-8 illustrate the effect of Nernst potential $E_{other}$ on
concentrations, electric potentials and currents under different leak
conductance. $E_{other}$ is an approximation to the potassium equilibrium
potential The red lines with circles, squares and triangles are denotes the
different leak conductance $g_{m}=10^{-6},10^{-5},10^{-3}$, correspondingly.
The black dash lines are the results by setting zero leakage, i.e. $g_{m}=0$.
And the blue dash lines denotes voltage clamp results where the electric
potential at inside $\phi_{N}$ is set to be zero and electric potential at
outside $\phi_{P}=-E_{other}$ in system (37).

First, consider the natural case, when shunt conductance $g_{m}=0$. Cytochrome
$c$ oxidase is not affected. The efficiency of complex IV is not changed by
the Nernst potential because $I_{leak}$ is always zero in this case. When
$g_{m}>0$, as the Nernst potential become more negative, the resistance for
proton pumping increases which leads to the decrease of proton pump
efficiency. Of course, when the leak conductance is large enough, the system
is nearly voltage clamped to the equilibrium potential for the leak. The
Kirchhoff clamp (red lines with triangles) is unlocked, removed by the large
leak conductance. The results are the same as the results of the voltage clamp
(blue dash lines). Fig. 6(d), Fig. 7(c)-(d), and Fig. (8) show the difference
in various quantities between voltage and UNclamped natural situations.

The dramatic effect of antibiotics on membrane properties was studied by [23].
Valinomycin has been studied in many other papers, which include [19, 16, 26,
18, 20, 94].

Figure 6: Concentration at equilibrium states with different $E_{other}$and
$g_{m}$. (a) Electron concentration $\rho_{e}$; (b) $[H]_{E}$; (c)
$[H]_{B}$;(d) $[H]_{X}$. Black dash line: $g_{m}=0$; Red line with circles:
$g_{m}=10^{-6}$; Red line with squares: $g_{m}=10^{-5}$; Red line with
triangles: $g_{m}=10^{-3}$; Blue dash lines: Voltage clamp.

Figure 7: Electric potential at equilibrium states with different
$E_{other}$and $g_{m}$. (a) $\phi_{E}$; (b) $\phi_{B}$; (c) $\phi_{X}$;(d)
$\phi_{P}$. Red line with squares: $g_{m}=10^{-5}$; Red line with triangles:
$g_{m}=10^{-3}$; Blue dash lines: Voltage clamp.

Figure 8: Evolution of current at equilibrium states with different
$E_{other}$and $g_{m}$. (a) $I_{N2E}/I_{e}$; (b) $I_{N2B}/I_{e}$; (c)
$I_{E2B}/I_{e}$; (d) $I_{E2X}/I_{e}$. (e) $I_{X2P}/I_{e}$. Red line with
squares: $g_{m}=10^{-5}$; Red line with triangles: $g_{m}=10^{-3}$; Blue dash
lines: Voltage clamp.

##  5 Discussion and Conclusion

As history separated chemistry and field theory, so it separated chemical
theory from devices. Almost all of chemical theory devalues the significance
of boundary conditions and flows. Almost none of chemical theory allows flows
from boundary to boundary. These bald statements are easy to confirm. Most
monographs and texts of chemical theory barely mention flow from boundaries,
and most deal with equilibrium zero flow systems, using those results to
discuss what happens when flow is not zero with certain inherent difficulties
(given the obvious inconsistencies involved). References are hard to find that
discuss spatially nonuniform boundary conditions or flows from boundary to
boundary driven by external sources. Power supplies are essential for most
engineering devices and they require different locations on boundaries to have
different potential. That is to say, they require spatially nonuniform
boundary conditions for potential with nonzero flows on the boundary.

Devices are important. Our entire electronic technology is built from devices
that function more or less the same way no matter where they are located
(within reasonable limits, it goes without saying. Nothing in engineering or
technology is true in general. Everything exists and functions only within
reasonable limits). It would obviously be useful if chemical systems could be
easily and routinely shaped into devices.

Devices depend on spatially complex boundary conditions. A device has inputs
and outputs with different locations and different boundary equations. If
inputs and outputs are at the same location, and have the same properties,
there is no device! Most devices have power supplies as well as inputs and
outputs. These supply flows of energy that allow the device to have well
defined input output relations that are robust, quite independent of what is
connected to the input or output of the device. A transfer function relates
input and output when they are related by a constant coefficient causal
ordinary differential equation in time.

Devices maintain these properties almost entirely by using electricity and
energy from power supplies. They are fundamentally nonequilibrium systems with
spatially nonuniform Dirichlet boundary conditions for the electrical
potential.

Electricity is used to make engineering devices for good reason. Electrical
potentials, currents, and electrical energy are described by the Maxwell
equations with greater precision over a wider range of conditions than almost
any other physical phenomenon.

The generality of the Maxwell equations in classical form is obscured because
they embody an outdated and seriously inadequate representation of dielectrics
and polarization. We say this with no disrespect for the enormous
contributions of Maxwell and Heaviside, et al. But a representation that
described dielectric measurements slower than some 0.05 sec (in the 1890’s)
cannot be expected to describe systems studied in the 2020’s that function on
atomic time scales of femtoseconds, let alone the much faster time scales of
visible and ultraviolet light ($3x10^{12}$ \- $7x10^{12}$) to $10^{16}Hz$ and
even higher energy radiation like x-rays ($10^{16}$ \- $10^{20}Hz$) and gamma
rays($>10^{20}Hz$).

The classical Maxwell equations use a single real number to describe how
charge moves when an electric field is changed. Charge moves (and is said to
polarize) on an enormous range of times scales, when electric fields are
applied to matter, that cannot usefully be described by the single dielectric
constant. The classical formulation of the Maxwell equations embody—as well as
depend on— the classical, but crude representation of dielectrics. They use a
dielectric constant in the very definition of key variables. Those variables
in fact depend on an out of date constitutive model. If time dependent,
complicated charge movements are present, as is always the case in liquids,
and in solids over the time scale of modern technology, the Maxwell equations
need to be rewritten to isolate the movements of material charge (with mass)
from other types of current (e.g., from the displacement current found
throughout space $\varepsilon_{0}\partial E/\partial t$). The rewritten “Core
Maxwell Equations” must then be joined to a description of polarization
showing how charge moves when electric and magnetic forces are applied [109].
That description is not very different from the stress strain description of
how mass moves when a force is applied. Mathematicians and physicists have
dealt with such stress strain relations in solids and complex fluids with many
types of flow (migration, convection, diffusion) and those methods (chiefly of
the theory of complex fluids, in its energetic variational flavor) can be
applied to the polarization phenomena of charge.

The Core Maxwell Equations have a special property not found in many other
field equations because of the displacement current term
$\varepsilon_{0}\partial E/\partial t$ that is universal, present everywhere
including inside atoms and in empty space, wherever the Maxwell equations are
valid.

The Core Maxwell Equations are of use even when the spatial and temporal
location of charge is not known because the displacement current term
$\varepsilon_{0}\partial E/\partial t$ is universal, It is present everywhere,
inside atoms and between stars. For this reason, the Core Maxwell Equations
(surprisingly) are very useful for any description of polarization and for any
type of charge movement, whether created or driven by electrodynamics,
diffusion, convection, sunlight (in solar cells), or even heat flow even
though they themselves do not describe polarization phenomena.

The Core Maxwell Equations have this special property because of Maxwell’s
Ampere law. This law has no counterpart in the field theories of mechanics.
This law allows electric and magnetic fields to flow in perfect vacuum so they
can create propagated waves we call ‘the light of the sun’, even though the
vacuum contains no charge (with mass).

The Maxwell Ampere law has a corollary: the divergence of total current is
always zero, everywhere to the same accuracy and in the same domain that the
Maxwell equations are valid. The Maxwell Ampere law implies a unique
definition of total current as the entire source (i.e., right hand side) of
curl B in Maxwell’s Ampere law no matter how material charge moves. In fact,
Maxwell’s Ampere law defines total current everywhere, including in a perfect
vacuum where there is no mass or charge with mass or movement of charge with
mass whatsoever. In a total vacuum, total current is $\varepsilon_{0}\partial
E/\partial t$). In the presence of matter total current is $J_{total-
current}=\ J+\ \varepsilon_{0}\partial E/\partial t$) where $J$ described the
movement of charge with mass no matter how fast, transient, or small it is.
This fundamental property of electricity was well known to Maxwell and his
followers, as is made clear on p.155 and p. 511-512 of [3]. and is central to
the discussion of current flow in the classic text of [17] and the work of
Landauer [43, 33]. The idea of total current is defined to focus attention on
these issues and is discussed in many papers cited in [107].

If total current is confined to a single one dimensional path, conservation of
total current becomes equality of total current everywhere (see Fig. 2 of [81]
for an extensive physical discussion with examples). If that current is
confined to a circuit, conservation of total current becomes Kirchhoff’s
current law [112]. In a particular system, it is easy to verify whether the
one dimensional approximation for current flow is accurate (enough). Just
measure the currents and see if they are equal (in unbranched systems) or
whether they sum to zero as Kirchhoff’s law requires at nodes in circuits in
general as Kirchhoff’s law requires. It would be interesting to learn to
specify both the necessary and sufficient conditions in abstract mathematical
terms (that can be evaluated before a particular system is specified) under
which Kirchhoff’s law is accurate (enough).

The Kirchhoff’s law just described is not quite the Kirchhoff’s law of
textbooks of electrodynamics or engineering. The textbook law is derived and
presented as valid for long times many orders of magnitude longer than the
time scales of electronic devices, let alone atomic motion. There should be no
misunderstanding of this crucial point. The literature of circuit design shows
that Kirchhoff’s law is used as an essential design tool or analog and digital
circuits and integrated circuits [80, 46, 100, 46, 21, 49, 71, 48] that can
even function on the $10^{-12}$ second time scale [86]. It seems clear to us
that the generalization of Kirchhoff’s law to total current either solves this
problem, or makes it moot, as you wish [91, 88], following the practice of
Maxwell himself, and his followers, according to p.155 and p. 511-512 of [3].

This paper combines Kirchhoff’s law for total current with a quite general
description of chemical reactions with a general description of ion and water
flow, using the EnVarA (energy variational) approach in the tradition of the
theory of complex fluids. We use this electro-osmotic framework to analyze a
coarse grained description of cytochrome $c$ oxidase in the tradition of
‘master equations’. Our model is built on the carefully constructed and well
analyzed models of many others, but here we analyze the master equations using
currents defined as in Kirchhoff’s law for total current.

Our analysis is in the tradition of engineering. It does not depend on details
of electron current flow. It does not compute properties of charges and their
interactions, except implicitly as defined by the Maxwell equations,
particularly Maxwell’s Ampere’s law $curl\ B=\mu_{0}\ J_{total-
current}=\mu_{0}\ J+\ \mu_{0}\varepsilon_{0}\partial E/\partial t$

It seems obvious to us that each of the systems of oxidative phosphorylation
and photosynthesis require a circuit analysis embedded in the theory of
complex fluids. Evolution has built structures that conduct electron currents,
as in our electronic technology and so should be analyzed by the
extraordinarily successful methods of electronic circuit analysis. Evolution
has used ion current flow and chemical reactions in addition to electron flow
so the circuit approach is embedded here using the theory of complex fluids.
Complex fluid theory is designed to combine a wide variety of flows and the
forces and energies that drive them in a mathematically consistent way. The
structures evolution uses to control these flows form the geometry—or anatomy
or histology, depending on the length scale—of the system. The channels and
transporters (and electron transport pathways) built by evolution form
boundary conditions that decorate (i.e., are located on) these structures and
thus describe how they work.

Cytochrome $c$ oxidase is called Complex IV for a reason. It is embedded in a
lipid bilayer, connected to electron pathways we approximate as wires (in the
engineering tradition), and is surrounded by electrolytes that store energy in
their electric and chemical potential fields that form a complex fluid.
Complex IV includes the cytochrome $c$ oxidase enzyme, electron pathways,
channels for protons and potassium ions, pathways for oxygen diffusion, and
the membrane that encloses it and ions that surround it. Each subsystem stores
energy and responds to energy gradients with different types of flows. Complex
fluid mathematics is designed to handle systems of this complexity, although
biological applications involve more preset structural complexity than in many
physical systems of fluids. Complex fluid theory treats all fields, flows and
boundary conditions—including spatially nonuniform conditions that power the
system as they power electronic devices—consistently. The electro-osmotic
extension of the theory hopefully joins the forces, flows and energies of
chemical reactions into this formulation, while preserving the mathematical
consistency of the original theory, without violating the traditions of
chemistry.

Our work is significantly incomplete and limited. We over-approximate several
important biophysical mechanisms, including the water-gate switch, and the
oxygen reduction mechanism. We are more than aware of the need for higher
resolution in later work, with specific atomic scale models that compute the
electric field and flows from underlying structures and chemical reactions on
time scales of displacement (capacitive) currents that have been so well
resolved in experiments of great difficulty. These currents are important in
understanding the switches and mechanisms by which cytochrome c oxidase
couples electron flow, oxidative chemical reactions, and proton flow to make
oxidative phosphorylation possible in mitochondria.

## Acknowledgement

This work was partially supported by the National Natural Science Foundation
of China no. 12071190 and Natural Sciences and Engineering Research Council of
Canada (NSERC). We also would like to thank American Institute of Mathematics
where this project initiated.

## References

  * [1] William Thompson and L Kelvin  “On the theory of the electric telegraph”  In _Proc. Royal Soc. London_ 7, 1855, pp. 382–399
  * [2] C Chree  “The Mathematical Theory of Electricity and Magnetism”  In _Nature_ 78.2031 Nature Publishing Group, 1908, pp. 537–538
  * [3] James Hopwood Jeans  “The mathematical theory of electricity and magnetism”  Cambridge University Press, 1908
  * [4] Vannevar Bush  “Operational circuit analysis”  In _New York_ , 1929, pp. 197
  * [5] E.A. Guillemin  “Communications Networks Vol. 1 The Classical Theory of Lumped Constant Networks”  John Wiley, 1931
  * [6] Kenneth Stewart Cole  “Dynamic electrical characteristics of the squid axon membrane”  In _Archives des sciences physiologiques_ 3.2, 1949, pp. 253–258
  * [7] Alan Hodgkin, Andrew Huxley and Bernhard Katz  “Ionic Currents underlying activity in the giant axon of the squid”  In _Arch. Sci. physiol._ 3, 1949, pp. 129–150
  * [8] A.L. Hodgkin, A.F. Huxley and B. Katz  “Measurement of current- voltage relations in the membrane of the giant axon of Loligo”  In _J. Physiol. (London)_ 116, 1952, pp. 424–448
  * [9] P.. Caldwell, A.. Hodgkin, R.. Keynes and T.. Shaw  “The effects of injecting ‘energy-rich’ phosphate compounds on the active transport of ions in the giant axons of Loligo”  In _The Journal of Physiology_ 152.3, 1960, pp. 561–590  URL: http://jp.physoc.org/content/152/3/561.short
  * [10] PC Caldwell, AL Hodgkin, RD Keynes and TI Shaw  “Partial inhibition of the active transport of cations in the giant axons of Loligo”  In _The Journal of Physiology_ 152.3, 1960, pp. 591
  * [11] Vitold Belevitch  “Summary of the history of circuit theory”  In _Proceedings of the IRE_ 50.5 IEEE, 1962, pp. 848–855
  * [12] Alexander Mauro  “Space Charge Regions in Fixed Charge Membranes and the Associated Property of Capacitance”  In _Biophysical Journal_ 2.2, Part 1, 1962, pp. 179–198  DOI: http://dx.doi.org/10.1016/S0006-3495(62)86848-9
  * [13] Alan Finkelstein and Alexander Mauro  “Equivalent Circuits as Related to Ionic Systems”  In _Biophysical Journal_ 3.3, 1963, pp. 215–237  DOI: http://dx.doi.org/10.1016/S0006-3495(63)86817-4
  * [14] Andrew F Huxley  “The quantitative analysis of excitation and conduction in nerve” 1963, Les Prix Nobel, 1963, pp. 242–260
  * [15] P.. Caldwell, A.. Hodgkin, R.. Keynes and T.. Shaw  “The Rate of Formation and Turnover of Phosphorus Compounds in Squid Giant Axons”  In _J Physiol_ 171, 1964, pp. 119–31  URL: http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=Citation&list_uids=14170137
  * [16] B.. Pressman  “Induced active transport of ions in mitochondria”  In _Proceedings of the National Academy of Sciences_ 53.5, 1965, pp. 1076–1083  DOI: 10.1073/pnas.53.5.1076
  * [17] S. Ramo, J.R. Whinnery and T. Van Duzer  “Fields and Waves in Communication Electronics”  J. Wiley, 1965  URL: https://books.google.com/books?id=ABpRAAAAMAAJ
  * [18] Thomas E Andreoli, M Tieffenberg and Daniel C Tosteson  “The effect of valinomycin on the ionic permeability of thin lipid membranes”  In _The Journal of general physiology_ 50.11, 1967, pp. 2527–2545
  * [19] Eric J Harris, Graham Catlin and Berton C Pressman  “Effect of transport-inducing antibiotics and other agents on potassium flux in mitochondria”  In _Biochemistry_ 6.5, 1967, pp. 1360–1370
  * [20] Daniel C Tosteson, Pl Cook, Thomas Andreoli and M Tieffenberg  “The effect of valinomycin on potassium and sodium permeability of HK and LK sheep red cells”  In _The Journal of general physiology_ 50.11, 1967, pp. 2513–2525
  * [21] Mohammed Shuaib Ghausi and John Joseph Kelly  “Introduction to distributed-parameter networks: with application to integrated circuits”  Holt, RinehartWinston, 1968
  * [22] Clifford Truesdell  “Rational thermodynamics: a course of lectures on selected topics”  McGraw-Hill, 1969
  * [23] J. Leung and R.. Eisenberg  “The effects of the antibiotics gramicidin A, amphotericin B, and nystatin on the electrical properties of frog skeletal muscle”  In _Biochim Biophys Acta_ 298.3, 1973, pp. 718–23  URL: http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=Citation&list_uids=4541500
  * [24] Robert S Eisenberg  “Electrophysiology: Electric Current Flow in Excitable Cells. JJB Jack, D. Noble, and RW Tsien. Clarendon (Oxford University Press)”  In _Science_ 190.4219 American Association for the Advancement of Science, 1975, pp. 1087–1087
  * [25] Stephen G Brush  “The kind of motion we call heat”  North-Holland Amsterdam, 1976
  * [26] Berton C Pressman  “Biological applications of ionophores”  In _Annual review of biochemistry_ 45.1, 1976, pp. 501–530
  * [27] Sidney Darlington  “A history of network synthesis and filter theory for circuits composed of resistors, inductors, and capacitors”  In _IEEE transactions on circuits and systems_ 31.1 IEEE, 1984, pp. 3–13
  * [28] B. Hille  “Transport Across Cell Membranes: Carrier Mechanisms, Chapter 2”  In _Textbook of Physiology_ 1 Philadelphia: Saunders, 1989, pp. 24–47.
  * [29] Daniel Tosteson  “Membrane Transport: People and Ideas”, 1989, pp. 414
  * [30] James Darnell, Harvey Lodish and David Baltimore  “Molecular Cell Biology”  New York: Scientific American Books, 1990
  * [31] R.S. Eisenberg  “Channels as enzymes: Oxymoron and Tautology”  In _Journal of Membrane Biology_ 115, 1990, pp. 1–12. Available on arXiv as http://arxiv.org/abs/1112.2363
  * [32] A.L. Hodgkin  “Chance and Design”  New York: Cambridge University Press, 1992, pp. 401
  * [33] Rolf Landauer  “Conductance from transmission: common sense points”  In _Physica Scripta_ 1992.T42, 1992, pp. 110  URL: http://stacks.iop.org/1402-4896/1992/i=T42/a=020
  * [34] Siyuan Han, J Lapointe and JE Lukens  “SUPERCONDUCTING QUANTUM INTERFERENCE DEVICES”  In _Activated Barrier Crossing: Applications in Physics, Chemistry and Biology_ 4 World Scientific, 1993, pp. 241
  * [35] B. Alberts et al.  “Molecular Biology of the Cell”  New York: Garland, 1994, pp. 1294
  * [36] Elizabeth Garber, Stephen G Brush and CWF Everitt  “Maxwell on Heat and Statistical Mechanics: On" Avoiding All Personal Enquiries" of Molecules”  Lehigh University Press, 1995
  * [37] R.S. Eisenberg  “Atomic Biology, Electrostatics and Ionic Channels.”  In _New Developments and Theoretical Studies of Proteins_ 7 Philadelphia: World Scientific, 1996, pp. 269–357. Published in the Physics ArXiv as arXiv:0807.0715
  * [38] R.S. Eisenberg  “Computing the field in proteins and channels.”  In _Journal of Membrane Biology_ 150, 1996, pp. 1–25. Preprint available on physics arXiv as document 1009.2857
  * [39] Robert F Pierret  “Semiconductor device fundamentals”  Pearson Education India, 1996
  * [40] Roger T. Howe and Charles G. Sodini  “Microelectronics: an integrated approach”  Upper Saddle River, NJ USA: Prentice Hall, 1997, pp. 908
  * [41] Bob Eisenberg  “Ionic channels in biological membranes-electrostatic analysis of a natural nanotube”  In _Contemporary Physics_ 39.6 Taylor & Francis, 1998, pp. 447–466
  * [42] Thomas K Simpson  “Maxwell on the Electromagnetic Field: A Guided Study.”  American Association of Physics Teachers, 1998
  * [43] Yoseph Imry and Rolf Landauer  “Conductance viewed as transmission”  In _Reviews of Modern Physics_ 71.2, 1999, pp. S306
  * [44] Steven E Laux and Karl Hess  “Revisiting the analytic theory of pn junction impedance: Improvements guided by computer simulation leading to a new equivalent circuit”  In _IEEE Transactions on Electron Devices_ 46.2 IEEE, 1999, pp. 396–412
  * [45] Amr Haggag and Karl Hess  “Analytical theory of semiconductor pn junctions and the transition between depletion and quasineutral region”  In _IEEE Transactions on Electron Devices_ 47.8 IEEE, 2000, pp. 1624–1629
  * [46] P.. Gray, P.. Hurst, S.. Lewis and R.. Meyer  “Analysis and Design of Analog Integrated Circuits”  New York: John Wiley, 2001, pp. 875
  * [47] A.F. Huxley  “From overshoot to voltage clamp”  In _Trends in Neurosciences_ 25.11, 2002, pp. 553–558
  * [48] Howard W Johnson and Martin Graham  “High-speed signal propagation: advanced black magic”  Prentice Hall Professional, 2003
  * [49] R.S. Muller, M. Chan and T.I. Kamins  “Device Electronics For Integrated Circuits, 3rd Ed”  Wiley India Pvt. Limited, 2003  URL: https://books.google.com/books?id=X9Q-QQAACAAJ
  * [50] M. Wikstrom, M.. Verkhovsky and G. Hummer  “Water-gated mechanism of proton translocation by cytochrome c oxidase”  In _Biochim Biophys Acta_ 1604.2, 2003, pp. 61–5  URL: http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=Citation&list_uids=12765763
  * [51] Dmitry Bloch et al.  “The catalytic cycle of cytochrome <em>c</em> oxidase is not the sum of its two halves”  In _Proceedings of the National Academy of Sciences of the United States of America_ 101.2, 2004, pp. 529–533  DOI: 10.1073/pnas.0306036101
  * [52] D.. Gadsby  “Ion transport: spot the difference”  In _Nature_ 427.6977, 2004, pp. 795–7  DOI: 10.1038/427795a
  * [53] J-P Colinge and Cynthia A Colinge  “Physics of semiconductor devices”  Springer Science & Business Media, 2005
  * [54] Paul Scherz and S. Monk  “Practical electronics for inventors”  McGraw-Hill, Inc., 2006, pp. 1056
  * [55] M.. Verkhovsky, I. Belevich, D.. Bloch and M. Wikstrom  “Elementary steps of proton translocation in the catalytic cycle of cytochrome oxidase”  In _Biochim Biophys Acta_ 1757.5-6, 2006, pp. 401–7  URL: http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=Citation&list_uids=16829227
  * [56] Ilya Belevich et al.  “Exploring the proton pump mechanism of cytochrome c oxidase in real time”  In _Proceedings of the National Academy of Sciences_ 104.8, 2007, pp. 2685–2690
  * [57] D.. Gadsby  “Structural biology: ion pumps made crystal clear”  In _Nature_ 450.7172, 2007, pp. 957–9  DOI: 10.1038/450957a
  * [58] Young C Kim, Mårten Wikström and Gerhard Hummer  “Kinetic models of redox-coupled proton pumping”  In _Proceedings of the National Academy of Sciences_ 104.7 National Acad Sciences, 2007, pp. 2169–2174
  * [59] H. Miedema et al.  “A Biological Porin Engineered into a Molecular, Nanofluidic Diode”  In _Nano Lett._ 7.9, 2007, pp. 2886–2891  URL: http://pubs3.acs.org/acs/journals/doilookup?in_doi=10.1021/nl0716808
  * [60] Henk Miedema et al.  “A biological porin engineered into a molecular, nanofluidic diode”  In _Nano letters_ 7.9 ACS Publications, 2007, pp. 2886–2891
  * [61] Walter Boron and Emile Boulpaep  “Medical Physiology”  New York: Saunders, 2008, pp. 1352
  * [62] Ryogo Sugitani, Emile S Medvedev and Alexei A Stuchebrukhov  “Theoretical and computational analysis of the membrane potential generated by cytochrome c oxidase upon single electron injection into the enzyme”  In _Biochimica et Biophysica Acta (BBA)-Bioenergetics_ 1777.9 Elsevier, 2008, pp. 1129–1139
  * [63] Marina L Verkhovskaya et al.  “Real-time electron transfer in respiratory complex I”  In _Proceedings of the National Academy of Sciences_ 105.10, 2008, pp. 3763–3767
  * [64] D.. Gadsby  “Ion channels versus ion pumps: the principal difference, in principle”  In _Nat Rev Mol Cell Biol_ 10.5, 2009, pp. 344–52  DOI: 10.1038/nrm2668
  * [65] Young C Kim, Mårten Wikström and Gerhard Hummer  “Kinetic gating of the proton pump in cytochrome c oxidase”  In _Proceedings of the National Academy of Sciences_ 106.33 National Acad Sciences, 2009, pp. 13707–13712
  * [66] Clifford. Prosser, Brian A. Curtis and Esmail Meisami  “A History of Nerve, Muscle and Synapse Physiology”  Stipes Public License, 2009, pp. 572
  * [67] Isaiah Shavitt and Rodney J Bartlett  “Many-body methods in chemistry and physics: MBPT and coupled-cluster theory”  Cambridge university press, 2009
  * [68] Ilya Belevich et al.  “Initiation of the proton pump of cytochrome c oxidase”  In _Proceedings of the National Academy of Sciences_ 107.43 National Acad Sciences, 2010, pp. 18469–18474
  * [69] Bob Eisenberg  “Computing the field in proteins and channels”  In _arXiv preprint arXiv:1009.2857_ , 2010
  * [70] Bob Eisenberg, Yunkyong Hyon and Chun Liu  “Energy variational analysis of ions in water and channels: Field theory for primitive models of complex ionic fluids”  In _The Journal of Chemical Physics_ 133.10 American Institute of Physics, 2010, pp. 104104
  * [71] Elya B. Joffe and Kai-Sang Lock  “Grounds for Grounding”  NY: Wiley-IEEE Press, 2010, pp. 1088
  * [72] Margareta R.. Blomberg and Per E.. Siegbahn  “The mechanism for proton pumping in cytochrome c oxidase from an electrostatic and quantum chemical perspective”  In _Biochimica et Biophysica Acta (BBA) - Bioenergetics_ 1817.4, 2012, pp. 495–505  DOI: 10.1016/j.bbabio.2011.09.014
  * [73] Georges Gielen and Willy MC Sansen  “Symbolic analysis for automated design of analog integrated circuits”  Springer Science & Business Media, 2012
  * [74] Young C Kim and Gerhard Hummer  “Proton-pumping mechanism of cytochrome c oxidase: A kinetic master-equation approach”  In _Biochimica et Biophysica Acta (BBA)-Bioenergetics_ 1817.4 Elsevier, 2012, pp. 526–536
  * [75] N. Sperelakis and N. Sperelakis  “Cell Physiology Source Book: Essentials of Membrane Biophysics”  Elsevier Science, 2012  URL: https://books.google.com/books?id=n7kxScnqQsQC
  * [76] Takefumi Yamashita and Gregory A Voth  “Insights into the mechanism of proton transport in cytochrome c oxidase”  In _Journal of the American Chemical Society_ 134.2 ACS Publications, 2012, pp. 1147–1152
  * [77] R.L. Boylestad and L. Nashelsky  “Electronic Devices and Circuit Theory: Pearson New International Edition PDF eBook”  Pearson Education, 2013  URL: https://books.google.com/books?id=zyypBwAAQBAJ
  * [78] Andrew Zangwill  “Modern electrodynamics”  Cambridge University Press, 2013
  * [79] J. Keener and J. Sneyd  “Mathematical Physiology: I: Cellular Physiology”  Springer New York, 2014  URL: https://books.google.com/books?id=Mf7HsgEACAAJ
  * [80] Paul Horowitz and Winfield Hill  “The Art of Electronics”  Cambridge University Press, 2015, pp. 1224
  * [81] Robert S. Eisenberg  “Mass Action and Conservation of Current”  In _Hungarian Journal of Industry and Chemistry Posted on arXiv.org with paper ID arXiv:1502.07251_ 44.1, 2016, pp. 1–28  DOI: 10.1515/hjic-2016-0001
  * [82] Gongchen Sun, Satyajyoti Senapati and Hsueh-Chia Chang  “High-flux ionic diodes, ionic transistors and ionic amplifiers based on external ion concentration polarization by an ion exchange membrane: a new scalable ionic circuit platform”  In _Lab on a Chip_ 16.7, 2016, pp. 1171–1177  DOI: 10.1039/c6lc00026f
  * [83] Joseph J Feher  “Quantitative human physiology: an introduction”  Academic press, 2017
  * [84] Sacha Friedli and Yvan Velenik  “Statistical mechanics of lattice systems: a concrete mathematical introduction”  Cambridge University Press, 2017
  * [85] Mi-Ho Giga, Arkadz Kirshtein and Chun Liu  “Variational modeling and complex fluids”  In _Handbook of mathematical analysis in mechanics of viscous fluids_ Springer International Publishing, Cham, 2017, pp. 1–41
  * [86] Lianfeng Zou, Shulabh Gupta and Christophe Caloz  “A Simple Picosecond Pulse Generator Based on a Pair of Step Recovery Diodes”  In _IEEE Microwave and Wireless Components Letters_ 27.5, 2017, pp. 467–469  DOI: 10.1109/lmwc.2017.2690880
  * [87] Xiuhong Cai et al.  “Network analysis of a proposed exit pathway for protons to the P-side of cytochrome c oxidase”  In _Biochimica et Biophysica Acta (BBA) - Bioenergetics_ 1859.10, 2018, pp. 997–1005  DOI: 10.1016/j.bbabio.2018.05.010
  * [88] Bob Eisenberg, Nathan Gold, Zilong Song and Huaxiong Huang  “What current flows through a resistor?”  In _arXiv preprint arXiv:1805.04814_ , 2018
  * [89] Mårten Wikström and Vivek Sharma  “Proton pumping by cytochrome c oxidase–a 40 year anniversary”  In _Biochimica et Biophysica Acta (BBA)-Bioenergetics_ 1859.9, 2018, pp. 692–698
  * [90] Shixin Xu, Bob Eisenberg, Zilong Song and Huaxiong Huang  “Osmosis through a semi-permeable membrane: a consistent approach to interactions”  In _arXiv preprint arXiv:1806.00646_ , 2018
  * [91] Robert S Eisenberg  “Kirchhoff’s Law Can Be Exact”  In _arXiv preprint arXiv:1905.13574_ , 2019
  * [92] Robert S Eisenberg  “Updating Maxwell with Electrons, Charge, and More Realistic Polarization”  In _arXiv preprint available at https://arxiv.org/abs/1904.09695_ , 2019
  * [93] D.U. Silverthorn et al.  “Human Physiology: An Integrated Approach”  Pearson Education, Incorporated, 2019  URL: https://books.google.com/books?id=JOPStAEACAAJ
  * [94] Zhangfei Su et al.  “How Valinomycin Ionophores Enter and Transport K+ across Model Lipid Bilayer Membranes”  In _Langmuir_ 35.51, 2019, pp. 16935–16943  DOI: 10.1021/acs.langmuir.9b03064
  * [95] Yi Zhu, Shixin Xu, Robert S. Eisenberg and Huaxiong Huang  “A Bidomain Model for Lens Microcirculation”  In _Biophysical Journal_ 116.6, 2019, pp. 1171–1184 Preprint available at https://arxiv.org/abs/1810.04162  DOI: https://doi.org/10.1016/j.bpj.2019.02.007
  * [96] Ren-Chuen Chen et al.  “Differential Capacitance of Electric Double Layers: A Poisson-Bikerman Formula”  In _arXiv preprint arXiv:2012.13141_ , 2020
  * [97] Robert S Eisenberg  “Electrodynamics Correlates Knock-on and Knock-off: Current is Spatially Uniform in Ion Channels”  In _arXiv preprint arXiv:2002.09012_ , 2020
  * [98] Robert s. Eisenberg  “Maxwell Equations for Material Systems”  In _doi: 10.20944/preprints202011.0201.v1_ , 2020  DOI: 10.20944/preprints202011.0201.v1
  * [99] Chin-Lung Li and Jinn-Liang Liu  “Generalized Debye–Hückel Equation From Poisson–Bikerman Theory”  In _SIAM Journal on Applied Mathematics_ 80.5 SIAM, 2020, pp. 2003–2023
  * [100] Jens Lienig and Juergen Scheible  “Fundamentals of layout design for electronic circuits”  Springer Nature, 2020
  * [101] Jinn-Liang Liu and Bob Eisenberg  “Molecular mean-field theory of ionic solutions: A Poisson-Nernst-Planck-Bikerman model”  In _Entropy_ 22.5 Multidisciplinary Digital Publishing Institute, 2020, pp. 550
  * [102] A.S. Sedra et al.  “Microelectronic Circuits”  Oxford University Press, Incorporated, 2020  URL: https://books.google.com/books?id=qGfRzQEACAAJ
  * [103] Lingyue Shen et al.  “An energy stable C0 finite element scheme for a quasi-incompressible phase-field model of moving contact line with variable density”  In _Journal of Computational Physics_ 405 Elsevier, 2020, pp. 109179
  * [104] Yiwei Wang, Chun Liu, Pei Liu and Bob Eisenberg  “Field theory of reaction-diffusion: Law of mass action with an energetic variational approach”  In _Physical Review E_ 102.6 APS, 2020, pp. 062147
  * [105] Yi Zhu, Shixin Xu, Robert S. Eisenberg and Huaxiong Huang  “A Tridomain Model for Potassium Clearance in Optic Nerve”  In _arXiv:2012.03303_ , 2020  DOI: arxiv:2012.03303
  * [106] Luigi Catacuzzeno, Luigi Sforna, Fabio Franciolini and Robert S. Eisenberg  “Multiscale modeling shows that dielectric differences make NaV channels faster than KV channels”  In _Journal of General Physiology_ 153 DOI: 10.1085/jgp.202012706.2, 2021  DOI: 10.1085/jgp.202012706
  * [107] Robert S Eisenberg  “Maxwell Equations Without a Polarization Field, Using a Paradigm from Biophysics”  In _Entropy_ 23.2 Multidisciplinary Digital Publishing Institute, 2021, pp. 172
  * [108] Budi Riza Putra et al.  “Microscale Ionic Diodes: An Overview”  In _Electroanalysis_ 33.6, 2021, pp. 1398–1418  DOI: 10.1002/elan.202060614
  * [109] Yiwei Wang, Chun Liu and Bob Eisenberg  “On variational principles for polarization in electromechanical systems”  In _arXiv e-prints_ , 2021, pp. arXiv–2108
  * [110] Yi Zhu, Shixin Xu, Robert S Eisenberg and Huaxiong Huang  “Membranes in Optic Nerve Models”  In _arXiv preprint arXiv:2105.14411_ , 2021
  * [111] Yi Zhu, Shixin Xu, Robert S Eisenberg and Huaxiong Huang  “Optic nerve microcirculation: Fluid flow and electrodiffusion”  In _Physics of Fluids_ 33.4, 2021, pp. 041906
  * [112] Robert Eisenberg  “A Necessary Addition to Kirchhoff’s Current Law of CircuitsVersion”  In _Engineering Archive EngArXiv_ https://doi.org/10.31224/2234, 2022  DOI: https://doi.org/10.31224/2234
  * [113] Lingyue Shen et al.  “An Energy Stable C^0 Finite Element Scheme for A Phase-Field Model of Vesicle Motion and Deformation”  In _SIAM Journal on Scientific Computing_ 44.1 SIAM, 2022, pp. B122–B145
  * [114] Zilong Song, Robert Eisenberg, Shixin Xu and Huaxiong Huang  “A Bubble Model for the Gating of K v Channels”  In _arXiv preprint arXiv:2204.13077_ , 2022
  * [115] Yiwei Wang and Chun Liu  “Some Recent Advances in Energetic Variational Approaches”  In _Entropy_ 24.5, 2022, pp. 721  DOI: 10.3390/e24050721

##  Appendix A Simulation results for case1

###  A.1 Effect of Oxygen

The effect of oxygen concentration in the reaction site is studied next. Fig.9
shows the dynamics (i.e., time dependence) of the concentration of ions in
different compartments at different oxygen concentrations. The dashed lines
are computed with default parameters shown in Tables 1-2. Panel (a) shows that
decreasing oxygen concentration at first decreases the reaction rate. The
decrease in reaction rate produces the accumulation of electrons (see panel
(b)) as they are supplied from the input source of constant electron current.
The proton concentration in the reaction site also increases: the reaction
rate decreases as the accumulated electron attracts more protons from the E242
site (see Fig. 10 (c)). Since the pump strength depends on the reaction rate,
the pump current $I_{Pump}$ also decreases (as shown in Fig. 10(a)). The
decrease in pump current induces the increase at E242 and the decrease at the
PLS (Proton Loading) site (see Fig. 10 (d)-(e)). Then the accumulated protons
and electrons modify the reaction rate $\mathcal{R}$ in Eq. (38j). The
reaction rate increases until $\mathcal{R}=-\frac{S_{v}I_{e}}{2F}$. That is in
fact the same equilibrium as determined by the default parameters.

At the same time, due to the accumulation of protons in E242, the chemical
potential $\mu_{E}$ is larger than then N side $\mu_{N}$. The result is an
activation of the rectifier making the current $I_{N2E}$ zero (see Fig. 10
(b)). This action depends on the rectifier, There is no rectifier between the
reaction site and N site, so the behavior is quite different.When the protons
are accumulated in B, the current $I_{N2B}$ is negative. Similarly, fewer
protons are pumped to the Protein Loading Site PLS site. The chemical
potential at that site $\mu_{X}$ is smaller than the P side $\mu_{P}$.
Negative current is also observed in Fig. 10 (f).

We confirm and extend the above observations by changing the oxygen
concentration from $10^{-6}\mu M$ to $10^{-2}\mu M$. As shown in Fig. 11 and
Fig. 12, the decrease of oxygen concentration (at the equilibrium) changes
just the concentrations of electron and proton at the BNC at the equilibrium
state. The reaction rate at equilibrium keeps constant due to constant
supplyment of electron flux $I_{e}$ and Eq.(37d). In this case, the proton
transportation follows the normal pattern as shown in Fig. 3 a.

Figure 9: Concentration with Different Oxygen concentration. (a) Reaction
rate; (b) Electron concentration $\rho_{e}$; (c) $[H]_{B}$; (d) $[H]_{E}$; (e)
$[H]_{X}$. The dash lines are results with default parameters.

Figure 10: Current with Different Oxygen concentration. (a) Pump $I_{pump}$;
(b) $I_{N2E}$; (c) $I_{N2B}$; (d) $I_{E2B}$; (e) $I_{E2X}$; (f) $I_{X2P}$; (g)
$I_{leak}$. The dash lines are results with default parameters.

Figure 11: Concentration at equilibrium states with different oxygen
concentration. (a) Electron concentration $\rho_{e}$; (b) $[H]_{E}$; (c)
$[H]_{B}$;(d) $[H]_{X}$.

Figure 12: Pump efficiency at equilibrium states with different oxygen
concentration $[H]_{N}$ and different threshold. (a) Reaction rate; (b)
$\mu_{X}-\mu_{E}$; (c) $I_{Pump}$;(d) $I_{E2X}/I_{e}$.

###  A.2 Effect of the Switch

In this section, the effect of the switch on the complex IV function is
studied when the oxygen level in the reaction site is very low $2.8E-5\mu M$
which is 1% of the default value. Especially, here we assume the switch is
perfect $SW_{0}=0$ and defective ones with $SW_{0}=1E-1,-2,-3fA$.

In Fig. 13, Panel (a) is consistent with the results of the last section: the
pump current decreases when the oxygen concentration is low. The pump current
recovers when electrons and protons accumulate. When the switch allows larger
counterflow, the pump current starts increases later but recovered faster to
the default value. The switch on the current $I_{N2E}$ is shown in Panel (b).
With small $SW_{0}=0,1E-2,-3$, the current is truncated if the counterflow
current magnitude is larger than $SW_{0}$. When the switch allows larger
counterflow, most of the accumulated protons at E242 flow back from E242 are N
side and less flow from E242 to reaction site B (Fig.13 (c)-(d)). And due the
decrease of chemical potential in E242, more protons are leaking back from
Protein Loading Site PLS site to the E242 (Fig.13 (e)) which induces more
protons flows from P side to the Protein Loading Site PLS site (Fig.13 (f)).
The corresponding dynamics of protons and electron are shown in Fig. 14.

Figure 13: Current with Oxygen concentration $2.8E-5\mu M$ and switch
threshold. (a) Pump $I_{pump}$; (b) $I_{N2E}$; (c) $I_{N2B}$; (d) $I_{E2B}$;
(e) $I_{E2X}$; (f) $I_{X2P}$; (g) $I_{leak}$. The dash lines are results with
default parameters.

Figure 14: Concentration with Different $SW_{0}$. (a) Electron concentration
$\rho_{e}$; (b) $[H]_{B}$; (c) $[H]_{E}$; (d) $[H]_{X}$. The dash lines are
results with default parameters.