Research

xxxxx, xxxxx, xxxx

A. Hillier

# On the Ambipolar Diffusion Formulation for Ion-neutral drifts in the non-
negligible drift velocity limit

A. S. Hillier1 1Department of Mathematics and Statistics, University of
Exeter, Exeter EX4 4QF U.K<EMAIL_ADDRESS>
###### Abstract

The ambipolar diffusion approximation is used to model partially ionised
plasma dynamics in a single fluid setting. To correctly apply the commonly
used version of ambipolar diffusion, a set of criteria should be satisfied
including the requirement that the difference in velocity between charges and
neutral species (known as drift velocity) is much smaller than the thermal
velocity, otherwise the drift velocity will drive a non-negligible level of
further collisions between the two species. In this paper we explore the
consequence of relaxing this assumption. We show that a new induction equation
can be formulated without this assumption. This formulation reduces to the
ambipolar induction equation in the case the drift velocity is small. In the
large drift velocity limit, the feedback of the drift velocity on the
collision frequency results in decreased diffusion of the magnetic field
compared with the standard ambipolar diffusion approximation for the same
parameters. This has a natural consequence of reducing the frictional heating
that can occur. Applying this to results from flux emergence simulations where
the expansion of the magnetic field leads to strong adiabatic cooling of the
partially ionised chromosphere resulted in a noticeable reduction in the
magnitude of the predicted drift velocities.

###### keywords:

xxxx, xxxx, xxxx

##  1 Introduction

In many astrophysical systems (e.g. stellar atmospheres, jets from young
stellar objects and protoplanetary disks), the conditions are such that there
is insufficient energy to fully ionise the plasma. Therefore, the majority of
the fluid is composed of neutral species. So even though magnetic fields are
recognised as being crucial in driving dynamics in these systems, under these
conditions, most of the species do not directly feel the Lorentz force of the
magnetic field. As a consequence, charged and neutral species feel different
forces. Consequently, this drives a drift in the velocity between these
different species [1, 2, 3, 4]. These physical processes are found to be
important in a wide range of astrophysical plasmas [5].

This drift between the charged and neutral species means that when the
particles collide, momentum is transferred between the plasma and neutral
fluids. Over time, this results in the momentum being redistributed such that
the two fluids move together. The collisions of the particles as they drift
past each other can result in frictional heating, meaning that coupling
process can be dynamically important for energy dissipation [6].

In situations where the dynamics of the neutrals and the ions happens at
sufficiently high frequency the neutrals and charges can act as separate
fluids, an example of this will be in the rapid transitions of physical
properties found in two-fluid shock fronts [7]. However, in many cases the
fluid is sufficiently coupled to be treated approximately as one fluid, but
with the drift between species taken into account via a modified Ohm’s law. In
this situation the induction equation becomes [8, 3]:

$\frac{\partial\mathbf{B}}{\partial t}=\nabla\times\left(\mathbf{v\times
B}+\eta_{\rm
AMB}\frac{(\mathbf{J}\times\mathbf{B})\times\mathbf{B}}{{B}^{2}}\right).$ (1)

There are a number of approximations that go into this diffusion
approximation, key ones are [8]:

  * •

The changes in the drift velocity occur at low frequency.

  * •

The temperature of the neutrals and the plasma is the same

  * •

The magnitude of the drift velocity is much smaller than the thermal velocity.

Of these approximations, it is worth highlighting the final in this list. This
is an often forgotten part of the approximations made when calculating the
collision frequency. Investigating the relaxation of this assumption is the
key goal of this paper.

In a low-$\beta$ regime (where plasma $\beta$ is the ratio of gas to magnetic
pressure) it is possible there will develop a strong Lorentz force that drives
a large enough drift velocity that the associated speed becomes greater than
the plasma thermal velocity. Large velocity drifts, evidenced by the magnitude
of the ambipolar diffusion term, have been observed in simulations of magnetic
flux emergence [9, 10]. These calculations have shown that the inclusion of
ambipolar diffusion drastically changed the rate at which the magnetic flux
emerged through the solar chromosphere due to the significant magnitude of the
diffusion. Similarly large ambipolar terms are found in the chromosphere
modelling simulations [11, 12, 13]. In these calculations large velocity
drifts were found to develop in the partially ionised chromosphere, including
in relation to magnetic flux emergence and spicule dynamics. With these large
velocity drifts significant heating was found. However, when ambipolar
diffusion gets large, and with it large velocity drifts are present in
simulations, there is the possibility that at least the third of the
assumptions behind the ambipolar diffusion approximation is invalidated.

To go beyond the first two of the listed approximations requires multi-fluid
modelling. For the final of these approximations, we can include the role of
super-thermal drift into the ambipolar diffusion approximation. In this paper,
we take advantage of the very accurate approximation that has been developed
for the collision frequency [14, 15] and change in collision cross-section
[16] to model the large drift velocity limit to present an extension to the
ambipolar diffusion term that allows it to hold for drift velocities greater
than the thermal velocity.

##  2 Determining the collisional velocity

Our first step is to look at the two-fluid momentum equations, so that we can
then simplify this into a single fluid model. These are

$\displaystyle\rho_{n}\frac{D\mathbf{v}_{\rm n}}{Dt}=-\nabla p_{\rm
n}-\rho_{n}\mathbf{g}-\alpha(\mathbf{v}_{\rm D},T_{\rm n},T_{\rm p})\rho_{\rm
n}\rho_{\rm p}\mathbf{v}_{\rm D},$ (2)
$\displaystyle\rho_{p}\frac{D\mathbf{v}_{\rm p}}{Dt}=-\nabla p_{\rm
p}+\frac{1}{c}\mathbf{J}\times\mathbf{B}-\rho_{p}\mathbf{g}+\alpha(\mathbf{v}_{\rm
D},T_{\rm n},T_{\rm p})\rho_{\rm n}\rho_{\rm p}\mathbf{v}_{\rm D},$ (3)

where $\mathbf{v}_{\rm D}=\mathbf{v}_{\rm n}-\mathbf{v}_{\rm p}$. By dividing
Equation 2 by $\rho_{n}$ and Equation 3 by $\rho_{p}$ then subtracting the
latter from the former approximately leads to this equation

$\frac{D\mathbf{v}_{\rm D}}{Dt}=\frac{1}{\rho_{p}}\nabla p_{\rm
p}-\frac{1}{c\rho_{p}}\mathbf{J}\times\mathbf{B}-\frac{1}{\rho_{n}}\nabla
p_{\rm n}-\alpha({v}_{\rm D},T_{\rm n},T_{\rm p})(\rho_{\rm n}+\rho_{\rm
p})\mathbf{v}_{\rm D}.$ (4)

The parameter $\alpha$ when multiplied by a density gives a collision
frequency, where these collision frequencies can either be calculated from
hard-sphere collisions between particles or through effective collisions
driven by charge exchange. For these two different effects, $\alpha$ takes
slightly different forms. When looking only at hard-sphere collisions,
$\alpha$ is given by [14, 17]

$\alpha\approx\frac{1}{m_{\rm n}+m_{\rm p}}\frac{8}{3}\sqrt{\frac{2k_{\rm
B}T_{\rm r}}{\pi m_{r}}}\left(1+\frac{9\pi s^{2}}{64}\right)^{1/2}\Sigma_{\rm
pn}({v}_{\rm D},T_{\rm n},T_{\rm p}),$ (5)

where

$\displaystyle T_{\rm r}=$ $\displaystyle\frac{m_{\rm p}T_{\rm n}+m_{\rm
n}T_{\rm p}}{m_{\rm n}+m_{\rm p}},$ (6) $\displaystyle m_{\rm r}=$
$\displaystyle\frac{m_{\rm n}m_{\rm p}}{m_{\rm n}+m_{\rm p}},$ (7)
$\displaystyle s=$ $\displaystyle v_{\rm D}\left(\frac{2k_{\rm B}T_{\rm
r}}{m_{r}}\right)^{-1/2},$ (8) $\displaystyle v_{\rm D}=$
$\displaystyle|\mathbf{v}_{\rm D}|=|\mathbf{v}_{\rm n}-\mathbf{v}_{\rm p}|.$
(9)

Here $\Sigma_{\rm pn}({v}_{\rm D},T_{\rm n},T_{\rm p})$ is the collision
velocity dependent cross-section (see, for example, [18] for details on the
collision velocity dependence) with $\Sigma_{\rm pn}=\Sigma_{\rm np}$. If one
is only interested in a hydrogen plasma, as we use in the investigations later
in this work, we can take $m_{\rm n}=m_{\rm p}=m_{\rm i}$, where $m_{\rm i}$
is the mass of a proton.

For charge exchange (in this case for hydrogen proton charge exchange) an
equivalent form is used [15]

$\alpha\approx\frac{1}{m_{\rm p}}\frac{8}{3}\sqrt{\frac{2k_{\rm B}T_{\rm
r}}{\pi m_{r}}}\left(1+\frac{9\pi s^{2}}{64}\right)^{1/2}\Sigma_{\rm
CX}({v}_{\rm D},T_{\rm n},T_{\rm p}).$ (10)

For the two different equations for $\alpha$, there are two key differences.
Firstly, this is that the cross section for hard sphere collisions and charge
exchange are different. Secondly, the factor ${1}/(m_{\rm n}+m_{\rm p})$ in
Equation 5 becomes ${1}/{m_{\rm p}}$ in Equation 10. Note that in the
literature there are other forms used for $\alpha$ when charge exchange is
considered [19, 20, 21, 22, 23]. As shown in Figure 1, the approximation used
in their works (solid green line), and the one adopted in ours (dashed blue
line) gives to very high level of accuracy the same dependence of
$\alpha/\Sigma$ for different drift velocities with at most $2$% difference
between the two models).

Figure 1: Plot of $\alpha/\Sigma$ (in units of cm/s) against ${v}_{\rm D}$
with $T=8000$ K for the approximation shown in equation 5 (dashed blue line)
and for the model applied in [19, 20, 21, 22, 23] for charge exchange (solid
green line).

##  3 Calculating $\mathbf{v}_{\rm D}$ in the strong coupling limit

In this section, for simplicity we assume that the cross-section is
independent of the drift velocity (many models take the collisional cross-
section to be a constant for example). We will look later at how the
requirement for the cross-section to independent of the drift velocity can be
relaxed. Here we will use the form of $\alpha$ suitable for hard sphere
collisions, but it is trivial to rework these derivations to take charge
exchange as the principle momentum transfer process into account.

The key assumption of ambipolar diffusion comes from the strong coupling limit
which implies that the time derivative of the drift velocity is negligible and
that the temperatures of the fluids must all be the same. By assuming, for
simplicity of the presentation of the mathematics, that we are dealing with a
hydrogen fluid (e.g. $m_{\rm p}=m_{\rm n}$), this leads to:

$\rho_{\rm T}\frac{1}{2m_{\rm p}}\frac{8}{3\sqrt{\pi}}\sqrt{\frac{4k_{\rm
B}T}{m_{p}}}\left(1+\frac{9\pi s^{2}}{64}\right)^{1/2}\Sigma_{\rm
pn}(T)\mathbf{v}_{\rm D}=-\frac{\mathbf{F}_{\rm
p}}{\rho_{p}}+\frac{\mathbf{F}_{\rm n}}{\rho_{n}},$ (11)

where we have taken (as is usual in the ambipolar diffusion formulation) that
$T_{\rm p}=T_{\rm n}=T_{\rm r}=T$, $\mathbf{F}_{\rm n}$ and $\mathbf{F}_{\rm
p}$ are denoting the force terms from the RHS of Equations 2 and 3
respectively, and taking that $\rho_{\rm T}\equiv\rho_{\rm p}+\rho_{\rm n}$.
This can be written as

$\left(1+\frac{9\pi s^{2}}{64}\right)^{1/2}s\hat{\mathbf{v}}_{\rm
D}=-X\hat{\mathbf{F}}_{\rho},$ (12)

where

$\mathbf{F}_{\rho}\equiv\frac{1}{\rho_{p}}\mathbf{F}_{\rm
p}-\frac{1}{\rho_{n}}\mathbf{F}_{\rm n}$ (13)

and with

$X=\frac{2m_{\rm p}}{\rho_{\rm T}}\frac{3\sqrt{\pi}}{8}\left({\frac{4k_{\rm
B}T}{m_{p}}}\right)^{-1}\Sigma_{\rm pn}(T)^{-1}|\mathbf{F}_{\rho}|,$ (14)

and the hat symbol denotes a unit vector implying

$\hat{\mathbf{v}}_{\rm D}=-\hat{\mathbf{F}}_{\rho}.$ (15)

Squaring both sides of Equation 12, we can write the following quadratic in
$s^{2}$

$\frac{9\pi}{64}s^{4}+s^{2}-X^{2}=0.$ (16)

This gives the following solution for $s^{2}$ (noting that $s^{2}\geq 0$)

$s^{2}=\frac{32}{9\pi}\left(\sqrt{1+\frac{9\pi}{16}X^{2}}-1\right).$ (17)

Therefore, $s$ is given by

$s=\sqrt{\frac{32}{9\pi}\left(\sqrt{1+\frac{9\pi}{16}X^{2}}-1\right)}.$ (18)

As such, our velocity drift is given by

$\mathbf{v}_{\rm
D}=-\sqrt{\frac{32}{9\pi}\left(\sqrt{1+\frac{9\pi}{16}X^{2}}-1\right)}\sqrt{\frac{4k_{\rm
B}T}{m_{\rm p}}}\hat{\mathbf{F}}_{\rho}.$ (19)

Figure 2: Plot of $s$ for a drift-velocity independent cross section against
$X$ (red line) as given by Equation 18. The two limits, small $X$ limit
(dashed line) with $s=X$ which is equivalent to the standard ambipolar
diffusion approximation and large $X$ limit (dash-triple dot line), are also
plotted.

Figure 2 shows the variation of $s$ against $X$ (red line). The dashed line
shows the curve expected for $s$ when $X^{2}$ is taken to be small, and the
triple-dot-dashed line shows the curve for $s$ when $X^{2}$ is assumed to be
large. The derivation of these two limits is explained in the rest of the
subsection. It is clear that as $X$ becomes greater than 1 the value of $s$
strongly diverges from that of the small-force limit.

###  3.1 The Ambipolar diffusion approximation in the Lorentz force dominant
limit

Generally for ambipolar diffusion in many astrophysical systems, the largest
term on the RHS of Equation 11 is the one given by the Lorentz force. In this
limit the drift velocity is given by:

$\mathbf{v}_{\rm D}=-v_{\rm
D}\frac{\mathbf{J}\times\mathbf{B}}{|\mathbf{J}\times\mathbf{B}|}.$ (20)

Using this drift velocity the induction equation can be rewritten in terms of
the neutral velocity

$\displaystyle\frac{\partial\mathbf{B}}{\partial t}=$
$\displaystyle\nabla\times(\mathbf{v}_{\rm
p}\times\mathbf{B})=\nabla\times(\mathbf{v}_{\rm
n}\times\mathbf{B}-\mathbf{v}_{\rm D}\times\mathbf{B})$ (21) $\displaystyle=$
$\displaystyle\nabla\times\left(\mathbf{v}_{\rm
n}\times\mathbf{B}\right)+\nabla\times\left(\sqrt{\frac{32}{9\pi}\left(\sqrt{1+\frac{9\pi}{16}X^{2}}-1\right)}\sqrt{\frac{4k_{\rm
B}T}{m_{\rm
p}}}\left(\frac{\mathbf{J}\times\mathbf{B}}{|\mathbf{J}\times\mathbf{B}|}\right)\times\mathbf{B}\right).$
(22)

An alternative way of formulating these equations can be to put the system
into the centre-of-mass reference frame [3] which gives

$\displaystyle\frac{\partial\mathbf{B}}{\partial t}=$
$\displaystyle\nabla\times(\mathbf{v}_{\rm
p}\times\mathbf{B})=\nabla\times(\mathbf{v}_{\rm CM}\times\mathbf{B}-\xi_{\rm
n}\mathbf{v}_{\rm D}\times\mathbf{B})$ (23) $\displaystyle=$
$\displaystyle\nabla\times\left(\mathbf{v}_{\rm
CM}\times\mathbf{B}\right)+\nabla\times\left(\xi_{\rm
n}\sqrt{\frac{32}{9\pi}\left(\sqrt{1+\frac{9\pi}{16}X^{2}}-1\right)}\sqrt{\frac{4k_{\rm
B}T}{m_{\rm
p}}}\left(\frac{\mathbf{J}\times\mathbf{B}}{|\mathbf{J}\times\mathbf{B}|}\right)\times\mathbf{B}\right).$
(24)

For the rest of the paper we will focus on the neutral frame-of-reference, but
as can be seen everything we present here can be easily adapted for systems
where the centre-of-mass frame of reference is more appropriate.

####  3.1.1 Small-force limit

To look into the behaviour of the drift velocity (and how it connects to
formulations of the induction equation, and the evolution of the magnetic
field this allows), it is informative to look at the formulation for $v_{\rm
D}$ in the limits of small and large Lorentz force.

The first limit we investigate is the limit of small forces, which can be
described as the situation where

$1\gg\frac{9\pi}{16}X^{2}.$ (25)

In this limit, we can perform a Taylor expansion, truncating at first order,
and on substitution into Equation 20 we find:

$\displaystyle\mathbf{v}_{\rm D}=$ $\displaystyle-X\left({\frac{4k_{\rm
B}T}{m_{\rm
p}}}\right)^{1/2}\frac{\mathbf{J}\times\mathbf{B}}{|\mathbf{J}\times\mathbf{B}|},$
(26) $\displaystyle=$ $\displaystyle-\frac{2m_{\rm p}}{\rho_{\rm
T}^{2}\xi_{\rm p}\Sigma_{\rm pn}(T)}\frac{3\sqrt{\pi}}{8}\left({\frac{4k_{\rm
B}T}{m_{\rm p}}}\right)^{-1/2}\frac{1}{c}\mathbf{J}\times\mathbf{B},$ (27)

which is equivalent to $s=X$. This is the formulation we would have expected
if we had made the assumption of small drift velocity in the original
derivation. In this case Equation 21 becomes:

$\displaystyle\frac{\partial\mathbf{B}}{\partial t}=$
$\displaystyle\nabla\times\left(\mathbf{v}_{\rm
n}\times\mathbf{B}+\frac{2m_{\rm p}}{\rho_{\rm T}^{2}\xi_{\rm p}\Sigma_{\rm
pn}(T)}\frac{3\sqrt{\pi}}{8}\left({\frac{4k_{\rm B}T}{m_{\rm
p}}}\right)^{-1/2}\frac{1}{c}(\mathbf{J}\times\mathbf{B})\times\mathbf{B}\right),$
(28) $\displaystyle=$ $\displaystyle\nabla\times\left(\mathbf{v}_{\rm
n}\times\mathbf{B}+\eta_{\rm
AMB}\frac{(\mathbf{J}\times\mathbf{B})\times\mathbf{B}}{B^{2}}\right),$ (29)

####  3.1.2 Large-force limit

An alternative limit to investigate, and one that takes our analysis into a
regime that is not accessible to standard ambipolar diffusion scenarios, is
the limit where the forces on the system are large, which implies we can have
large velocity drifts. We can simply state this as the regime where

$1\ll\sqrt{\frac{9\pi}{16}}X.$ (30)

In this case

$\displaystyle\mathbf{v}_{\rm D}=$ $\displaystyle-\sqrt{\frac{4kT}{\pi m_{\rm
p}}}\sqrt{\frac{8X}{3\sqrt{\pi}}}\frac{\mathbf{J}\times\mathbf{B}}{|\mathbf{J}\times\mathbf{B}|},$
(31) $\displaystyle=$ $\displaystyle-\sqrt{\frac{2m_{\rm p}}{\pi\rho_{\rm
T}^{2}\xi_{\rm p}\Sigma_{\rm
pn}(T)}}\frac{1}{\sqrt{c}}\frac{\mathbf{J}\times\mathbf{B}}{\sqrt{|\mathbf{J}\times\mathbf{B}|}},$
(32)

which again can be derived from Equation 11 in this case setting $X\gg 1$ in
Eqn 18 resulting in $s=\sqrt{8X/3\sqrt{\pi}}$. As we can clearly see, in this
limit the velocity drift scales as $\sqrt{|\mathbf{J}\times\mathbf{B}|}$
instead of $|{\mathbf{J}\times\mathbf{B}|}$ for the small-force limit. In this
limit, the induction equation becomes

$\displaystyle\frac{\partial\mathbf{B}}{\partial t}=$
$\displaystyle\nabla\times\left(\mathbf{v}_{\rm
n}\times\mathbf{B}+\sqrt{\frac{2m_{\rm p}}{\pi\rho_{\rm T}^{2}\xi_{\rm
p}\Sigma_{\rm
pn}(T)}}\frac{1}{\sqrt{c}}\frac{\mathbf{J}\times\mathbf{B}}{\sqrt{|\mathbf{J}\times\mathbf{B}|}}\times\mathbf{B}\right).$
(33)

###  3.2 Consequences for ambipolar diffusion heating

Based on the summation of the two-fluid energy internal energy equations [3],
the change in the total internal energy (as relevant for a single-fluid energy
equation) as a result of frictional heating ($H$) is given by:

$\displaystyle H=$ $\displaystyle\alpha(\mathbf{v}_{\rm D},T)\rho_{\rm
T}^{2}\xi_{\rm p}\xi_{\rm n}{v}_{\rm D}^{2},$ (34) $\displaystyle=$
$\displaystyle\frac{\rho_{\rm T}^{2}\xi_{\rm n}\xi_{\rm p}}{m_{\rm
p}}\frac{128}{27}\left(\frac{4k_{\rm B}T}{\pi m_{\rm
p}}\right)^{3/2}\Sigma_{\rm pn}(T)\frac{1}{\sqrt{2}}\left(1+\sqrt{1+\frac{9\pi
X^{2}}{16}}\right)^{1/2}\left(\sqrt{1+\frac{9\pi X^{2}}{16}}-1\right).$ (35)

where

$\frac{1}{\gamma-1}\frac{\partial p}{\partial t}=H,$ (36)

when there is only frictional heating in the system and $p=p_{\rm n}+p_{\rm
p}$. If we take

$Y(m_{\rm p},\rho_{\rm T},\xi_{\rm n},\xi_{\rm p},T)=\frac{\rho_{\rm
T}^{2}\xi_{\rm n}\xi_{\rm p}}{m_{\rm p}}\frac{128}{27}\left(\frac{4k_{\rm
B}T}{\pi m_{\rm p}}\right)^{3/2}\Sigma_{\rm pn}(T),$ (37)

then we can write this as

$\displaystyle H=$ $\displaystyle\frac{1}{\sqrt{2}}Y\left(1+\sqrt{1+\frac{9\pi
X^{2}}{16}}\right)^{1/2}\left(\sqrt{1+\frac{9\pi X^{2}}{16}}-1\right).$ (38)

Again we can look at the different limits, first taking the limit of small
$X^{2}$. This gives:

$\displaystyle H=$ $\displaystyle Y\frac{9\pi}{32}X^{2}$ (39) $\displaystyle=$
$\displaystyle\frac{m_{\rm p}\xi_{n}}{\rho_{\rm T}^{2}\xi_{\rm
p}}\frac{3\sqrt{\pi}}{4}\left(\frac{4k_{\rm B}T}{m_{\rm
p}}\right)^{-1/2}\Sigma_{pn}^{-1}\left|\frac{1}{c}\mathbf{J\times
B}\right|^{2}$ (40)

Alternatively, in the large force limit

$\displaystyle H=$
$\displaystyle\frac{Y}{\sqrt{2}}\left(\frac{9\pi}{16}\right)^{3/4}X^{3/2}$
(41) $\displaystyle=$ $\displaystyle 4\left(\frac{3}{2}\right)^{3/2}\xi_{\rm
n}\left(\frac{m_{\rm p}}{\rho_{\rm T}^{2}\xi_{\rm
p}}\right)^{1/2}\Sigma_{pn}^{-1/2}\left|\frac{1}{c}\mathbf{J\times
B}\right|^{3/2}.$ (42)

Therefore, when the Lorentz force in the system is small the heating scales
with $|\mathbf{J}\times\mathbf{B}|^{2}$. However, for large Lorentz force the
heating then scales as $|\mathbf{J}\times\mathbf{B}|^{3/2}$. Though not as
extreme a difference as the velocity difference scaling (and with it the
magnitude of ambipolar diffusion) with $|\mathbf{J}\times\mathbf{B}|$ changing
to $|\mathbf{J}\times\mathbf{B}|^{1/2}$. But this still illustrates an
important consequence, in low-$\beta$ systems the heating rates that are
calculated using the standard ambipolar diffusion approximation for single
fluid models (i.e. the small force limit of this paper) can easily be an
overestimate if the role of velocity drift in driving collisions is neglected.

##  4 Beyond singles species and constant cross-sections

The previous example looked at a fluid with only a single species (nominally
hydrogen but the analysis would apply to any species), and where the
collisional cross-section was independent of velocity drift. In this section
we present some extensions to this model.

###  4.1 Investigating the role of a non-constant collisional cross-section

Until now, we have worked on the assumption that the collisional cross-section
is independent of the drift velocity. However, this is not going to be the
case. For example, in the case of charge exchange, the cross section for
neutral-proton interactions can be approximated by the form [20]

$\Sigma_{\rm CX}=1.12\times 10^{-14}-7.15\times 10^{-16}\ln(v_{\rm
COLL}/100){\,\mathrm{c}m}^{2},$ (43)

with

$v_{\rm COLL}=\sqrt{\frac{16kT}{\pi m_{p}}+v_{\rm D}^{2}},$ (44)

note this form for the cross section is taken as a fit to data [16] with the
accuracy of the fit estimated as 10%. Looking at this form, it is clearly
awkward to manipulate algebraically, so if we return to the data for the
charge exchange cross section [16], by looking at their figure we can see that
then the magnitude of the cross section looks as if it varies as an
approximate power law over a large range of the data. By fitting the data of
[16] in the range of $v_{\rm COLL}=4.79\times 10^{5}$ to $3.68\times 10^{7}$
cm s-1 with a function of the form $\Sigma_{\rm CX}=Av_{\rm COLL}^{-1/n}$ with
$A$ a positive constant $n$ a positive integer, we attempt to approximate the
change in charge exchange cross section as a function of velocity. We found:

$\Sigma_{\rm CX}=10^{-12.85}v_{\rm COLL}^{-1/4}{\,\mathrm{cm}}^{2}=\Sigma_{\rm
CX,0}(T)\left(1+\frac{\pi}{4}s^{2}\right)^{-1/8},$ (45)

with $\Sigma_{\rm CX,0}(T)$ the cross section at zero drift velocity. Note
that this fit is within 5% of the atomic data. The accuracy can be seen in
Figure 3 (and has the benefit over that of [20] of being in the form of a
power law and as such allows easier algebraic manipulation). Due to the
restrictions in fitting, more accurate fits are possible, but the accuracy of
the fit used here is more than sufficient for the work. It is important to
note that more accurate models of the cross section exist, for example those
published in [18], where quantum mechanical effects are taken into account
resulting in non-monotonic changes in the cross-section at energies below 1
ev. However, we believe that the cross-section used here is sufficiently
accurate for our purposes.

Figure 3: Plot of $\Sigma_{\rm CX}$ against $v_{\rm COLL}$. The black line
gives the experimental data from [16] and the red line represents the fit used
in this study.

We can rewrite Equation 11 (assuming the only force is the Lorentz force) but
including a drift-velocity dependent cross-section giving

$\frac{\rho_{\rm T}}{m_{\rm p}}\frac{8}{3\sqrt{\pi}}\sqrt{\frac{4k_{\rm
B}T}{m_{p}}}\left(1+\frac{9\pi s^{2}}{64}\right)^{1/2}\Sigma_{\rm
CX,0}(T)\left(1+\frac{\pi}{4}s^{2}\right)^{-1/8}\mathbf{v}_{\rm
D}=-\frac{1}{\rho_{\rm p}}\mathbf{J\times B},$ (46)

or

$\left(1+\frac{9\pi
s^{2}}{64}\right)^{1/2}\left(1+\frac{\pi}{4}s^{2}\right)^{-1/8}s\hat{\mathbf{v}}_{\rm
D}=-X_{\rm CX}\frac{\mathbf{J\times B}}{|\mathbf{J\times B}|},$ (47)

with $X_{\rm CX}$ defined as

$X_{\rm CX}=\frac{m_{\rm p}}{\rho_{\rm T}^{2}\xi_{\rm
p}}\frac{3\sqrt{\pi}}{8}\left({\frac{4k_{\rm
B}T}{m_{p}}}\right)^{-1/2}\Sigma_{\rm CX,0}(T)^{-1}\left|\frac{\mathbf{J\times
B}}{c}\right|.$ (48)

By rearranging and putting to the 8th power we find

$\left(1+\frac{9\pi s^{2}}{64}\right)^{4}s^{8}=X_{\rm
CX}^{8}\left(1+\frac{\pi}{4}s^{2}\right).$ (49)

Rearranging Equation 49 for $s$ gives an 8th order polynomial in $s^{2}$. The
roots of this equation can be solved numerically, and then taking that we
require $s^{2}$ to be both purely real and positive to give physical solutions
we can find the appropriate solution and then determine $s$. We can also look
again at the limits, in this case when $s$ is either small or large. For
$s^{2}\ll 1$ we have

$s=X_{\rm CX}.$ (50)

Alternatively, for $s^{2}\gg 1$ we have

$s=X_{\rm
CX}^{4/7}\left({\frac{64}{9\pi}}\right)^{2/7}\left(\frac{\pi}{4}\right)^{1/14}.$
(51)

Therefore, unlike the constant cross section case, we have the dependence of
$v_{\rm D}$ on $|\mathbf{J\times B}|$ as scaling as $|\mathbf{J\times
B}|^{4/7}$ in the large drift velocity limit (equivalent to the large force
limit used earlier) instead of $|\mathbf{J\times B}|^{1/2}$ found in the case
of a constant cross section. Note this implies the scaling of the heating rate
with the magnitude of the Lorentz force of $|\mathbf{J\times B}|^{12/7}$.
Figure 4 shows $s$ against $X_{\rm CX}$, the red line is the solution for $s$,
and the dashed and triple dot dash lines are the small and large $s$ solutions
respectively.

Figure 4: Plot of $s$ for a drift-velocity dependent cross section against
$X_{\rm CX}$ (red line) calculated by taking the appropriate solution of
Equation 49. The two limits, small $s$ limit (dashed line) and large $s$ limit
(dash-triple dot line), are also plotted.

##  5 Summary and Conclusion

In this paper, we have provided a new analytic framework to relax a key
assumption behind the ambipolar diffusion formulation, something which is used
to model the influence of neutrals in the MHD induction equation. Namely, we
have removed the requirement that the difference in velocities between the
charged and neutral fluids has to be much less than the thermal velocity of
the plasma. The key consequence of removing this assumption is that it shows
that using the standard ambipolar diffusion approximation may overestimate the
magnitude of the diffusion of the magnetic field and the heating this
produces.

A key question remains, and this is: will there ever be a situation in the
solar atmosphere where the velocity drift gets large enough that the newly
proposed formulation is required? To date, the measured drift velocities,
though of order 1 km s-1 (for example [24]), have all been too small to
provide a meaningful departure from the ambipolar diffusion coefficient
predicted by the small drift velocity limit. However, simulations have shown
that there is the potential for huge velocity drifts. For example in the
simulations (including the standard ambipolar diffusion approximation) of
emerging flux by [25], they found drift velocities of $\sim 40$ km s-1 in
regions where the temperature had dropped to $\sim 2,000$ K. Putting these
parameters into Equation 8 gives a value of $s\approx 5$, which from Figure 2
we can see corresponds to a value of $X=5$ as in this limit $s=X$. Taking into
account that this large drift velocity would drive further particle collisions
acting as a limiter to the ambipolar diffusion we find a more accurate value
of $s\approx 2.5$ (equivalent to a drift velocity of $\sim 20$ km s-1). Though
this drift velocity is still very significant, it represents a non-trivial
reduction in magnitude.

It is reasonable to expect that a detailed investigation of the output of
single-fluid partially ionised plasma models of the solar atmosphere would
highlight examples where the changes in magnitude of the diffusion were even
larger. For example, in two-fluid simulations of magnetic reconnection [26,
27] it has been necessary to use this modified form of $\alpha$ as the effect
of drift-induced collisions became significant, so an important future work is
to assess during which physical processes the assumption of small drift
velocity breaks down. Beyond the solar atmosphere, there are many
astrophysical systems with low temperatures and densities. If this is combined
with being low plasma $\beta$, then the modified formulation presented here
may be applicable.

AH was supported in this research by his STFC Ernest Rutherford Fellowship
grant number ST/L00397X/2 and STFC research grants ST/R000891/1 and
ST/V000659/1. AH would like to acknowledge the discussions with members of
ISSI Team 457 “The Role of Partial Ionization in the Formation, Dynamics and
Stability of Solar Prominences”, which have helped improve the ideas in this
manuscript.

## References

  * [1] Khomenko E, Collados M, Díaz A, Vitas N. 2014 Fluid description of multi-component solar partially ionized plasma. Physics of Plasmas 21, 092901. (10.1063/1.4894106)
  * [2] Khomenko E. 2017 On the effects of ion-neutral interactions in solar plasmas. Plasma Physics and Controlled Fusion 59, 014038. (10.1088/0741-3335/59/1/014038)
  * [3] Khomenko E. 2020 pp. 69–116. In Multi-Fluid Extensions of MHD and Their Implications on Waves and Instabilities, pp. 69–116. Cham: Springer International Publishing. (10.1007/978-3-030-16343-3_3)
  * [4] Martínez-Sykora J, De Pontieu B, Hansteen V, Carlsson M. 2015 The role of partial ionization effects in the chromosphere. Philosophical Transactions of the Royal Society of London Series A 373, 20140268–20140268. (10.1098/rsta.2014.0268)
  * [5] Ballester JL, Alexeev I, Collados M, Downes T, Pfaff RF, Gilbert H, Khodachenko M, Khomenko E, Shaikhislamov IF, Soler R, Vázquez-Semadeni E, Zaqarashvili T. 2018 Partially Ionized Plasmas in Astrophysics. Space Sci. Rev. 214, 58. (10.1007/s11214-018-0485-6)
  * [6] Khomenko E, Collados M. 2012 Heating of the Magnetized Solar Chromosphere by Partial Ionization Effects. ApJ 747, 87. (10.1088/0004-637X/747/2/87)
  * [7] Hillier A, Takasao S, Nakamura N. 2016 The formation and evolution of reconnection-driven, slow-mode shocks in a partially ionised plasma. A&A 591, A112. (10.1051/0004-6361/201628215)
  * [8] Braginskii SI. 1965 Transport Processes in a Plasma. Reviews of Plasma Physics 1, 205.
  * [9] Leake JE, Arber TD. 2006 The emergence of magnetic flux through a partially ionised solar atmosphere. A&A 450, 805–818. (10.1051/0004-6361:20054099)
  * [10] Arber TD, Haynes M, Leake JE. 2007 Emergence of a Flux Tube through a Partially Ionized Solar Atmosphere. ApJ 666, 541–546. (10.1086/520046)
  * [11] Martínez-Sykora J, De Pontieu B, Carlsson M, Hansteen V. 2016 On the Misalignment between Chromospheric Features and the Magnetic Field on the Sun. ApJ 831, L1. (10.3847/2041-8205/831/1/L1)
  * [12] Martínez-Sykora J, De Pontieu B, Hansteen VH, Rouppe van der Voort L, Carlsson M, Pereira TMD. 2017a On the generation of solar spicules and Alfvénic waves. Science 356, 1269–1272. (10.1126/science.aah5412)
  * [13] Martínez-Sykora J, De Pontieu B, Carlsson M, Hansteen VH, Nóbrega-Siverio D, Gudiksen BV. 2017b Two-dimensional Radiative Magnetohydrodynamic Simulations of Partial Ionization in the Chromosphere. II. Dynamics and Energetics of the Low Solar Atmosphere. ApJ 847, 36. (10.3847/1538-4357/aa8866)
  * [14] Draine BT. 1986 Multicomponent, reacting MHD flows. MNRAS 220, 133–148. (10.1093/mnras/220.1.133)
  * [15] Zank GP, Adhikari L, Zhao LL, Mostafavi P, Zirnstein EJ, McComas DJ. 2018 The Pickup Ion-mediated Solar Wind. ApJ 869, 23. (10.3847/1538-4357/aaebfe)
  * [16] Barnett CF, Hunter HT, Fitzpatrick MI, Alvarez I, Cisneros C, Phaneuf RA. 1990 Atomic data for fusion. Volume 1: Collisions of H, H2, He and Li atoms and ions with atoms and molecules. Technical report.
  * [17] Oliver R, Soler R, Terradas J, Zaqarashvili TV. 2016 Dynamics of Coronal Rain and Descending Plasma Blobs in Solar Prominences. II. Partially Ionized Case. ApJ 818, 128. (10.3847/0004-637X/818/2/128)
  * [18] Vranjes J, Krstic PS. 2013 Collisions, magnetization, and transport coefficients in the lower solar atmosphere. A&A 554, A22. (10.1051/0004-6361/201220738)
  * [19] Pauls HL, Zank GP, Williams LL. 1995 Interaction of the solar wind with the local interstellar medium. J. Geophys. Res. 100, 21595–21604. (10.1029/95JA02023)
  * [20] Meier ET, Shumlak U. 2012 A general nonlinear fluid model for reacting plasma-neutral mixtures. Physics of Plasmas 19, 072508. (10.1063/1.4736975)
  * [21] Leake JE, Lukin VS, Linton MG, Meier ET. 2012 Multi-fluid Simulations of Chromospheric Magnetic Reconnection in a Weakly Ionized Reacting Plasma. ApJ 760, 109. (10.1088/0004-637X/760/2/109)
  * [22] Leake JE, Lukin VS, Linton MG. 2013 Magnetic reconnection in a weakly ionized plasma. Physics of Plasmas 20, 061202. (10.1063/1.4811140)
  * [23] Popescu Braileanu B, Lukin VS, Khomenko E, de Vicente Á. 2019 Two-fluid simulations of waves in the solar chromosphere. I. Numerical code verification. A&A 627, A25. (10.1051/0004-6361/201834154)
  * [24] Zapiór M, Heinzel P, Khomenko E. 2022 Doppler-velocity Drifts Detected in a Solar Prominence. ApJ 934, 16. (10.3847/1538-4357/ac778a)
  * [25] Nóbrega-Siverio D, Moreno-Insertis F, Martínez-Sykora J, Carlsson M, Szydlarski M. 2020 Nonequilibrium ionization and ambipolar diffusion in solar magnetic flux emergence processes. A&A 633, A66. (10.1051/0004-6361/201936944)
  * [26] Murtas G, Hillier A, Snow B. 2021 Coalescence instability in chromospheric partially ionized plasmas. Physics of Plasmas 28, 032901. (10.1063/5.0032236)
  * [27] Murtas G, Hillier A, Snow B. 2022 Collisional ionization and recombination effects on coalescence instability in chromospheric partially ionized plasmas. Physics of Plasmas 29, 062302. (10.1063/5.0087667)


Storage capacity of networks with discrete synapses and sparsely encoded
memories
Yu Feng1 and Nicolas Brunel1,2
1 Department of Physics, Duke University, Durham, NC 27710
2 Department of Neurobiology, Duke University, Durham, NC 27710

Abstract

Attractor neural networks (ANNs) are one of the leading theoretical frameworks
for the formation and retrieval of memories in networks of biological neurons.
In this framework, a pattern imposed by external inputs to the network is said
to be learned when this pattern becomes a fixed point attractor of the network
dynamics. The storage capacity is the maximum number of patterns that can be
learned by the network. In this paper, we study the storage capacity of fully-
connected and sparsely-connected networks with a binarized Hebbian rule, for
arbitrary coding levels. Our results show that a network with discrete
synapses has a similar storage capacity as the model with continuous synapses,
and that this capacity tends asymptotically towards the optimal capacity, in
the space of all possible binary connectivity matrices, in the sparse coding
limit. We also derive finite coding level corrections for the asymptotic
solution in the sparse coding limit. The result indicates the capacity of
networks with Hebbian learning rules converges to the optimal capacity
extremely slowly when the coding level becomes small. Our results also show
that in networks with sparse binary connectivity matrices, the information
capacity per synapse is larger than in the fully connected case, and thus such
networks store information more efficiently.

##  1 Introduction

It is widely believed that memories are stored in the brain through synaptic
modifications in an activity-dependent way. This idea has been implemented in
attractor neural network models, where the connectivity strength between
neurons is determined by Hebbian synaptic plasticity rules [1, 2]. In this
framework, a pattern is said to be learned if it becomes a fixed point
attractor of the network dynamics. An extensively studied question is, how
many patterns can be stored in such networks? Classical studies of memory
modeling synapses as continuous variables in networks of binary neurons have
shown that such networks can store a number of uncorrelated random patterns
$p$ that scales linearly that network size, $p_{max}=\alpha_{c}N$ where
$\alpha_{c}$ is of order 1 in the large $N$ limit [3, 2, 4]. However, there is
evidence suggesting that synapses in brain structures involved in memory, such
as the hippocampus and neocortex, are more digital than analog [5, 6, 8, 7].

A number of studies have addressed the question of the storage capacity of
networks with discrete synapses. Krauth and Mézard showed that networks can
potentially have a large capacity when all synapses are required to be binary
[9], using Gardner’s approach [4], with an upper bound for capacity
$\alpha_{cmax}=0.83$ discrete synapses, instead of $\alpha_{cmax}=2$ for
continuous synapses. Sompolinsky studied the storage capacity of a network
with a specific binarized Hebbian rule [10, 11], and showed its capacity is
remarkably close to the capacity of the Hopfield network [1], whose synapses
are continuous variables ($\alpha_{c}=0.10$ instead of 0.14). However, these
authors only studied the unbiased case, in which the coding level $f$ (i.e.
fraction of active neurons in a pattern) is 0.5, while neuronal activity in
areas involved in memory is typically very sparse (e.g. [12])- for instance,
the coding level in the human medial temporal lobe has been estimated to be
around 1% [13]. The upper bound for capacity in networks with arbitrary coding
levels and discrete synapses was computed by Gutfreund and Stein [14]. The
capacity in networks with Hebbian plasticity and binary synapses has only been
computed in the unbiased case, and the capacity for arbitrary coding levels
remains an open question. In this paper, we generalized Sompolinsky′s
calculation on the Hopfield model with binary synapses when coding level
$f=0.5$, to the model with fully connected or sparsely connected binary
synaptic connectivity with arbitrary coding levels. Our results show that the
network with binarized Hebbian rule has a similar capacity as the model with
continuous synapses for any coding level, and that this capacity tends
asymptotically towards the optimal capacity obtained by Gutfreund and Stein
[14], in the space of all possible binary connectivity matrices. Our results
also show that a network with sparse binary connectivity can have a larger
information capacity per synapse than a fully connected network, and thus can
allow a network to store information more efficiently.

##  2 Results

###  2.1 Storage capacity of fully-connected network with binary synapses

We consider a fully-connected neural network with $N$ binary (0,1) neurons.
The activity of neuron $i$ ($i=1,...,N$) is described by a binary variable,
$V_{i}=0,1$. Each neuron is connected to other neurons through the
connectivity matrix $W$. The activity of neuron $i$ at time $t$ is determined
by the asynchronous update rule (see Appendix VI for details):

$\displaystyle V_{i}\left(t+1\right)=\Theta[h_{i}\left(t\right)-\theta],$ (1)

where

$\displaystyle h_{i}\left(t\right)=\sum_{j\neq
i}^{N}W_{ij}V_{j}\left(t\right),$ (2)

is the local field, defined as the total input of neuron $i$, where $\theta$
is an activation threshold (constant independent of $V_{i}$), and $\Theta$ is
the Heaviside function. We also consider in Appendix a more general case of
dynamics with stochastic updates characterized by a temperature $T$, but we
focus in the main text in the zero temperature, deterministic limit.

The storage capacity of the network whose dynamics is defined by Eq.(1) and
Eq.(2) is determined by the connectivity matrix $W$. This connectivity matrix
$W$ depends on $p$ random uncorrelated patterns $\vec{\eta}^{\mu}$,
$\mu=1,...,p$, that are described by independent Bernoulli random variables:

$\displaystyle
P\left(\eta_{i}^{\mu}\right)=f\delta_{\left(\eta_{i}^{\mu},1\right)}+\left(1-f\right)\delta_{\left(\eta_{i}^{\mu},0\right)},$
(3)

where $\delta$ is the Kronecker delta function, and where $f$ is the coding
level (the fraction of active neurons). The storage capacity $\alpha$ is
defined as the maximal number of stored patterns $p$ divided by the network
size $N$, $\alpha=p/N$.

In this paper, we construct connectivity matrix $W$ from the patterns
$\vec{\eta}^{\mu}$ using a ‘clipped’ learning rule:

$\displaystyle
W_{ij}=\frac{\sqrt{p}}{N}F\left(\frac{1}{f(1-f)\sqrt{p}}\sum_{\mu=1}^{p}\left(\eta^{\mu}_{i}-f\right)\left(\eta^{\mu}_{j}-f\right)\right),$
(4)

where $F$ is given by:

$\displaystyle F\left(x\right)=\sqrt{\frac{\pi}{2}}\mbox{sign}(x),$ (5)

where the prefactor $\sqrt{\pi/2}$ is used for convenience. Thus, $W_{ij}$ is
only allowed to take two distinct values. With the nonlinear function $F(x)$
given by Eq. (5), $W_{ij}$ can be positive or negative. In neurobiological
networks, synaptic weights are sign-constrained, and their sign depends on
whether the presynaptic neuron is excitatory or inhibitory. The network with
the connectivity matrix given by Eqs. (4,5) leads to local fields of the form

$h_{i}=\sum_{j}W_{ij}V_{j}=\frac{\sqrt{2\pi
p}}{N}\sum_{j}\Theta(x_{ij})V_{j}-\sqrt{\frac{\pi
p}{2}}\frac{1}{N}\sum_{j}V_{j},$ (6)

where $\Theta(x)$ is the Heaviside function, and $x_{ij}$ is the argument of
$F$ in Eq. (4). Eq. (6) shows that the network is equivalent to a purely
excitatory network with binary weights (the first term in the r.h.s of Eq.
(6)) with an instantaneous linear inhibition (the second term in the r.h.s of
Eq. (6)).

Notice also that when $F(x)=x$, Eq.(4) yields the learning rule in the model
of Tsodyks and Feigel’man (TF) [2], where $W_{ij}$ is a continuous variable.
Therefore, we can interpret the learning rule of Eq.(4) as first learning
patterns $\vec{\eta^{\mu}}$ using the TF learning rule, and then clipping the
weight into discrete values at the end of the learning phase. In the
following, we call this model the Clipped Tsodyks-Feigel’man (CTF) model.

The nonlinearity of $F(x)$ in Eq.(5) makes the storage capacity more difficult
to calculate than the one of a network with a linear learning rule. In 1986,
Sompolinsky introduced a method to compute the storage capacity of Hopfield
networks with non-linear learning rules [10, 11]. In particular he showed that
in the large $N$ limit, these networks are equivalent to a linear learning
rule with an added random Gaussian noise,

$\displaystyle
W_{ij}=\frac{J}{Nf(1-f)}\sum_{\mu}^{p}\left(\eta_{i}^{\mu}-f\right)\left(\eta_{j}^{\mu}-f\right)+\delta_{ij},$
(7)

where $J$ is an embedding strength, and $\delta_{ij}$ is a random symmetric
Gaussian matrix. Both $J$ and the variance of the random Gaussian matrix
$\Delta_{0}^{2}={N\langle\delta_{ij}^{2}\rangle}/({J^{2}}{\alpha})$ can be
calculated as a function of $F(x)$ [11]. For $F(x)$ given by Eq.(5), the
embedding strength $J$ and $\Delta_{0}^{2}$ are given by (see details in
Appendix):

$J=1,\quad\Delta^{2}_{0}=\frac{\pi}{2}-1.$ (8)

####  2.1.1 Calculation of the storage capacity for arbitrary coding level $f$

To compute the storage capacity of a network with a learning rule given by
Eq.(7), we use standard methods and introduce the Hamiltonian

$H=\frac{1}{2}\sum_{i\neq j}W_{ij}V_{i}V_{j}+\theta\sum_{i}V_{i},$ (9)

where $W_{ij}$ is given Eq.(7). The typical free energy of the system can be
derived using the replica method [15, 3, 10]. The calculation allows us to
derive order parameters characterizing the system (such as the overlap of
network state with stored patterns), and its storage capacity. Using a replica
symmetric ansatz, the free energy of the system can be characterized by five
order parameters $m,Q,q,R,r$, where

$\begin{split}m&=\frac{1}{N}\sum_{i}\tilde{\eta}_{i}^{1}V_{i},\\\
Q&=\frac{1}{N}\sum_{i}V_{i},\\\ q&=\frac{1}{N}\sum_{i}V_{i}^{2},\end{split}$
(10)

and where $R$ and $r$ are conjugate variables of $Q$ and $q$, and are defined
by Eqs.(43) in Appendix II. The order parameter $m$ measures the retrieval
quality of a pattern stored in memory. Solutions with
$\tilde{m}\equiv\frac{m}{f(1-f)}\sim 1$ represent ‘retrieval states’ in which
the network goes to a fixed point close to one of the stored patterns.
Solutions with $\tilde{m}=0$ corresponds to no retrieval. The order parameter
$Q$ and $q$ represent the average neural activity and square of neural
activity of the network (see Appendix II for more details). The mean-field
equations of the system are obtained using a saddle point method. The full
equations are given by Eq.(48) in Appendix II for arbitrary coding levels and
temperature. In the zero-temperature limit, the equations simplify to the
following set of equations:

$\displaystyle\tilde{m}=\Phi(a_{1})-\Phi(a_{2}),$ (11)
$\displaystyle\tilde{r}=f\Phi(a_{1})+(1-f)\Phi(a_{2}),$ $\displaystyle
a_{1}=\frac{\tilde{\theta}-(1-f)\tilde{m}-Y}{\sqrt{\tilde{r}\alpha(1+\Delta_{0}^{2}(1-C)^{2})}},$
$\displaystyle
a_{2}=\frac{\tilde{\theta}+f\tilde{m}-Y}{\sqrt{\tilde{r}\alpha(1+\Delta_{0}^{2}(1-C)^{2})}},$
$\displaystyle Y=\frac{\alpha Cf}{2(1-C)}+\frac{1}{2}\alpha Cf\Delta_{0}^{2},$
$\displaystyle
C=\frac{f}{2\pi\alpha\tilde{r}}(fe^{-a_{1}^{2}/2}+(1-f)e^{-a_{2}^{2}/2})$

where
$\tilde{m}=m/f(1-f),\tilde{\theta}=\theta/f,\Delta_{0}^{2}=\Delta^{2}/\alpha$
and $\Phi\left(x\right)=\int_{x}^{\infty}Dz$.

Figure 1: Overlap as a function of storage load $\alpha$ for the Tsodyks-
Feigel’man model (TF) and the clipped Tsodyks-Feigel’man model (CTF), for
$f=0.02$. The solid lines represent the theoretical prediction and squares
represent the simulation results with a network of size $N=4000$ (mean and
standard deviation computed over five independent realizations). Dashed lines
mark the storage capacity for CTF and TF. For both models, the neuronal
activity threshold $\theta$ is chosen as the one that optimizes capacity. For
the parameters chosen here ($f=0.02$), $\tilde{\theta}\equiv\theta/f\approx
0.6$. Figure 2: Comparison between storage capacity of Hebbian rules and the
respective upper bounds. GUBC: Gardner Upper Bound for networks with
Continuous weights [4]. GUBB: Gutfreund and Stein Upper Bound for networks
with Binary weights [14]. See more details about capacity upper bounds in
Appendix IV. (A) Storage capacity of TF and CTF models as a function of coding
level $f$. When $f\to 0$, both capacities increase as $1/f\ln(1/f)$. Cyan and
magenta lines represent the storage capacity for a sparsely connected network
with continuous weights (STF) and a network with binary weights (SCTF). The
storage capacities of fully connected and sparsely connected networks converge
when $f\to 0$. Notice that for fully-connected network, storage capacity is
defined as $\alpha_{c}=p/N$ while for the sparsely-connected network, storage
capacity is defined as $\alpha_{s}=p/cN$, where $c\ll 1$.(B) Comparison
between storage capacity and upper bounds as a function of coding level $f$.
Even at $f=10^{-2}$, for both TF and CTF models, the capacity is only around a
third of the respective bounds, and thus the asymptotic solution Eq.(12) is
approached very slowly. (C) Comparison between the numerical solution and
asymptotic solutions. Solid lines are the numerical solutions of TF and CTF
models, the dotted lines with the same color are the corresponding asymptotic
solutions in the sparse coding limit (Eq.(13)), and dashed lines represent
asymptotic solutions with finite coding level corrections (Eq.(14,15)). (D)
Stored information per synapse as a function of coding level. When the coding
level $f$ goes to 0, the information stored in synapses increases but with an
extremely slow rate for both TF and CTF models. Dotted lines represents stored
information of asymptotic solutions in the sparse coding limit (i.e.,
$I(\alpha=(2f|\mbox{ln}f|)^{-1})$ for continuous weights case and
$I(\alpha=(\pi f|\mbox{ln}f|)^{-1})$ for binary weights case).

Once $f,\alpha,\theta$ are given, Eq.(11) can be solved numerically to obtain
the order parameters, including the rescaled overlap with the retrieved
pattern, $\tilde{m}$. Fig.1 shows $\tilde{m}$ as a function of $\alpha$ for
both TF and CTF models for $f=0.02$. The figure shows that analytical results
are in good agreement with simulations, using a network with $4,000$ neurons.
The maximal capacity of the network $\alpha_{c}$ is given by the largest value
of $\alpha$ for which there exist retrieval states (i.e. states with non-zero
$\tilde{m}$), optimized over the threshold $\theta$. Fig.1 shows that the
maximal capacity of the CTF model is lower than the one of the TF model, as
expected, but only by a factor of about 1.5. This maximal capacity
$\alpha_{c}$ is plotted as a function of $f$ in Fig.2A.

####  2.1.2 Sparse coding limit

In the biologically relevant sparse coding limit $f\to 0$, the mean-field
Eq.(49) take a rather simple form:

$\displaystyle\tilde{m}=\Phi\left(\frac{\tilde{\theta}-\tilde{m}(1-f)}{\sqrt{\tilde{r}\alpha(1+\Delta^{2}_{0})}}\right)-\Phi\left(\frac{\tilde{\theta}+\tilde{m}f}{\sqrt{\tilde{r}\alpha(1+\Delta^{2}_{0})}}\right),$
(12)
$\displaystyle\tilde{r}=f\Phi\left(\frac{\tilde{m}f}{\sqrt{\tilde{r}\alpha(1+\Delta^{2}_{0})}}\right)+(1-f)\Phi\left(\frac{\tilde{\theta}}{\sqrt{\tilde{r}\alpha(1+\Delta^{2}_{0})}}\right),$

where $\tilde{\theta}=\theta/f$ is a rescaled threshold, $\tilde{r}=r/f^{2}$,
and $\Phi\left(x\right)=\int_{x}^{\infty}Dz$ is the complementary cumulative
distribution function of the standard Gaussian distribution.

With additional analysis (see details in Appendix III), we find that the
maximum capacity in this limit is obtained when $\tilde{\theta}\sim 1$, and
the maximum capacity is:

$\begin{split}\alpha_{c}\simeq\frac{1}{\pi f\left|\log f\right|}.\end{split}$
(13)

which can be compared with the capacity of TF model obtained by [2],
$\alpha_{c}\simeq 1/(2f|\log f|)$. Thus, the two capacities differ by a factor
$\pi/2\sim 1.57$. We next address the question of how close this capacity is
to an upper bound in the space of all possible binary connectivity matrices.
This upper bound was computed by Gutfreund and Stein [14] for arbitrary coding
levels. The $f\rightarrow 0$ behavior could not be determined in a simple form
in that paper, however it was shown that the upper bound must be smaller or
equal than $1/(\pi f|\log f|)$. Our asymptotic result, Eq. (13), indicates
that the upper bound for binary connectivity matrices is indeed asymptotically
$1/(\pi f|\log f|)$, and thus the clipped TF model becomes asymptotically
optimal in the sparse coding limit. This is similar to what happens in
networks with continuous synapses, for which it was shown that the storage
capacity of the TF model tends to the upper bound of storage capacity obtained
by Gardner [4], in the space of all continuous synaptic matrices in the sparse
coding limit. Thus, in spite of their remarkable simplicity, both TF and CTF
models provide close to optimal learning rules for models with continuous and
discrete weights, respectively. For convenience, we summarize the storage
capacity of different models in Table 1 [9, 1, 10, 2, 4].

| capacity of Hebbian rule | Upper Bound
---|---|---
$f=0.5$, continuous $W_{ij}$ | $\sim 0.14$ | $2$
$f=0.5$, binary $W_{ij}$ | $\sim 0.1$ | $\sim 0.83$
$f\to 0$, continuous $W_{ij}$ | $(2f|\log f|)^{-1}$ | $(2f|\log f|)^{-1}$
$f\to 0$, binary $W_{ij}$ | $(\pi f|\log f|)^{-1}$ | $(\pi f|\log f|)^{-1}$
Table 1: Comparison between the storage capacity of Hebbian rules and upper
bounds computed using Gardner approach.

####  2.1.3 Leading correction to sparse coding limit

We can see that the capacities of networks with Hebbian rules in the unbiased
case ($f=0.5$) are much smaller than the corresponding upper bounds, while
they converge to the corresponding upper bounds in the sparse coding limit.
However, solving numerically mean-field Eqs.(49) for finite coding level $f$,
we find that the capacities of TF and CTF converge to the upper bounds
extremely slowly (see Fig. 2B). As shown in Fig.2B, the capacity of Hebbian
rules is only around $1/3$ compared to their corresponding upper bounds when
the coding level $f=10^{-2}$. This is mainly because the optimal threshold
$\tilde{\theta}$ approaches 1 extremely slowly when $f$ decreases (for
instance, for $f=0.02$, $\tilde{\theta}\sim 0.6$). With additional analysis
(see Appendix III), we derived the leading correction to the asymptotic
solution at the finite coding level :

$\begin{split}\alpha_{c}\simeq\frac{\tilde{\theta}^{2}_{opt}}{\pi f|\log
f|},\end{split}$ (14)

where the optimal threshold $\tilde{\theta}_{opt}$ is obtained by solving the
equation

$\begin{split}\frac{2\tilde{\theta}_{opt}^{2}\left|\log(1-\tilde{\theta}_{opt})\right|}{(1-\tilde{\theta}_{opt})^{2}}=\left|\log
f\right|.\end{split}$ (15)

Notice that when coding level $f\to 0$, $\tilde{\theta}_{opt}\to 1$ and
Eq.(14) recovers the asymptotic scaling of Eq.(13). The asymptotic solutions
Eq.(13) and Eq.(14) are compared in Fig.2C, and we can see that Eq.(14) agrees
with the numerical solutions very well when the coding level $f$ is small.
This result indicates that, in the biological sparse coding limit (i.e.,
coding level is small but finite), the capacities of models with Hebbian rules
are still notably smaller than the maximal capacities, in the space of all
possible connectivity matrices.

While the storage capacity in terms of numbers of patterns stored per synapse
diverges in the sparse coding limit, the information stored per pattern
decreases in that limit since it is proportional to the binary entropy of $f$.
As a result, the total information stored per synapse remains finite in the
sparse coding limit, both in the TF model [2] and in the corresponding Gardner
bound [4]. Fig.2D shows the information capacity in bits per synapse for
different models as a function of $f$,

$\begin{split}I=-\frac{\alpha}{\ln 2}\left(f\ln
f+(1-f)\ln(1-f)\right)\end{split}$ (16)

We find that when the coding level $f$ decreases, the information capacity of
TF and CTF increases quickly, while the corresponding upper bounds decrease
slowly. When $f$ goes to 0, the information capacity $I$ of TM and CTM further
increase and eventually converge to the optimal information capacity, but the
convergence rate is extremely low.

###  2.2 Storage capacity of sparsely-connected network with binary synapses

Cortical networks are characterized by low connection probabilities between
neurons (e.g. [16]). In the case we interpret the low synaptic efficacy state
to be zero, the network we have studied so far has a 50% connection
probability, much higher than observed connection probabilities in cortex,
which are of order 10% for excitatory neurons at short distances ($<100\mu$m).
This motivates the study of networks with sparser connectivity. Here, we study
two cases - one in which sparse connectivity is uncorrelated with learning,
and the other one where sparse connectivity is an outcome of learning with a
high synaptic threshold.

####  2.2.1 Sparse connectivity uncorrelated with learning

We first consider the case where learning occurs on top of a sparse random
Erdos-Renyi ‘structural’ connectivity matrix,

$\displaystyle
W_{ij}=\frac{c_{ij}\sqrt{p}}{Nc}F\left(\frac{1}{f(1-f)\sqrt{p}}\sum_{\mu=1}^{p}\left(\eta^{\mu}_{i}-f\right)\left(\eta^{\mu}_{j}-f\right)\right),$
(17)

where $F(x)$ is the same as the clipped function Eq.(5) for fully-connected
case, and $c_{ij}=1,0$ is a random binary matrix, with

$\displaystyle
P\left(c_{ij}\right)=c\delta_{\left(c_{ij},1\right)}+\left(1-c\right)\delta_{\left(c_{ij},0\right)},$
(18)

where $0<c\ll 1$ is the connection probability. The storage capacity of
learning rule Eq.(17) can be calculated similarly as the model in [17] in the
sparse connectivity limit $c\ll 1$, and the mean-field equations for finite
coding level $f$ are given as (see details in Appendix V):

$\begin{split}&\tilde{m}=\Phi\left(\frac{\tilde{\theta}-\tilde{m}(1-f)}{\sqrt{\alpha_{s}q(1+\Delta_{0}^{2})}}\right)-\Phi\left(\frac{\tilde{\theta}+\tilde{m}f}{\sqrt{\alpha_{s}q(1+\Delta_{0}^{2})}}\right),\\\
&q=f\Phi\left(\frac{\tilde{\theta}-\tilde{m}(1-f)}{\sqrt{\alpha_{s}q(1+\Delta_{0}^{2})}}\right)+(1-f)\Phi\left(\frac{\tilde{\theta}+\tilde{m}f}{\sqrt{\alpha_{s}q(1+\Delta_{0}^{2})}}\right),\end{split}$
(19)

where $\alpha_{s}=p/cN$ and where other order parameters are defined in
Eq.(10). The numerical solution of Eq.(19) are compared with fully-connected
case in Fig.2A. In the sparse coding case, the mean-field Eq.(19) coincide
with the mean-field Eq.(12), as expected [17]. We see that the capacity, in
terms of number of patterns stored divided by number of connections per
neuron, is larger in the sparsely connected case than in the fully connected
case, as expected from previous results in networks with continuous synapses
[18, 17].

####  2.2.2 Sparse connectivity induced by learning

In this section, we consider the case where sparse connectivity is obtained by
adding a threshold to the clipped function. Here we generalize the clipped
function Eq. (5) to:

$\begin{split}F_{T}(x)=\sqrt{2\pi}(\Theta(x-T)-M),\end{split}$ (20)

where M is given by

$\begin{split}M=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\Theta_{T}(x)e^{-x^{2}/2}dx,\end{split}$
(21)

where $\Theta_{T}(x)=\Theta(x-T)$, and $T$ is the threshold that can be used
to increase the sparseness of network connectivity, interpreting the low
synaptic state as a 0 state. With such a connectivity matrix, the connection
probability $R_{1}$ is given by:

$\begin{split}R_{1}=\frac{1}{2}\left(1-\mbox{erf}(\frac{T}{\sqrt{2}})\right)\end{split}$
(22)

In this case, the embedding strength $J$ and additional noise introduced by
clipped function Eq.(20)
$\Delta_{0}^{2}={N\langle\delta_{ij}^{2}\rangle}/{J^{2}}{\alpha}$ are

$\begin{split}J=e^{-T^{2}},\Delta_{0}^{2}=\frac{\pi}{2}e^{T^{2}}\left(1-\mbox{erf}\left(\frac{T}{\sqrt{2}}\right)^{2}\right)-1\end{split}$
(23)

Notice that the noise strength
$\Delta_{0}^{2}={N\langle\delta_{ij}^{2}\rangle}/{J^{2}}{\alpha}$ is an
increasing function of $T$. The storage capacity of the learning rule
Eq.(4,20) is determined by mean-field Eq.(12) for a given connection
probability $R_{1}$ and coding level $f$. As shown in Fig.3A, the storage
capacity $\alpha_{c}$ decreases when the threshold increases, and consequently
the connection probability $R_{1}$ decreases. Unsurprisingly, the storage
capacity decreases as the fraction of non-zero synapses decreases. However,
storing information with a smaller number of synapses also carries benefits in
terms of efficiency of information storage.

To quantify this efficency, we calculate the information stored in the network
per non-zero synapse, $I_{e}$:

$\begin{split}I_{e}=-\frac{\alpha}{R_{1}\ln 2}\left(f\ln
f+(1-f)\ln(1-f)\right)\end{split}$ (24)

The relation between $I_{e}$, $R_{1}$ and $f$ is shown in Fig.3B. We see that
the information capacity per synapse increases when the excitatory
connectivity becomes more sparse. This is mainly because we only keep
connections for which the Hebbian term (i.e. the argument of the clipping
function $F$ in Eq. (20) is large. In this way, the network can encode
information more efficiently when excitatory connections become sparse. Note
that maximizing storage capacity subject to a constraint of minimizing the
fraction of active synapses would lead to an optimal connection probability
$R_{1}^{*}$, whose precise value would depend on the cost of maintenance of an
active synapse. One can define a cost function $C=\alpha-\lambda R_{1}$, where
second term represents the cost of maintenance of active synapses. For a given
$\lambda$, we can obtain $R_{1}^{*}$ that maximize $C$. As shown in Fig.4B,
the more costly synaptic maintenance is, the sparser the resulting
connectivity. And the optimal connection probability $R_{1}^{*}\sim 0.1$ when
$\lambda\sim 10$.

Figure 3: Storage capacity and information capacity for a network with sparse
excitatory binary connections. (A) Storage capacity as a function of
connection probability $R_{1}$ and coding level $f$. The storage capacity
decreases when $R_{1}$ decreases. (B) Information capacity per active synapse,
as a function of $f$ and $R_{1}$. The information per synapse $I_{e}$
increases when $f$ and $R_{1}$ decrease. This indicates that for learning rule
Eqs.(4,20), both sparse coding and sparse connectivity can improve the coding
efficiency of the network. This result also indicates that the network can
have an optimal $R_{1}$ to balance storage capacity and coding efficiency.
Figure 4: Connection probability that optimizes capacity subject to a synapse
maintenance cost $\lambda$. (A) Cost function $C$ as a function of $R_{1}$ and
$\lambda$. (B) Optimal connection probability $R_{1}^{*}$ for different
$\lambda$. We can see that the more costly synaptic maintenance is, the
sparser the resulting connectivity. In both (A) and (B), coding level is set
to be $f=0.01$.

##  3 Discussion

We have calculated the storage capacity of an attractor neural network endowed
with binary synaptic weights at arbitrary coding levels. Our results show that
a network with a binarized Hebbian learning rule has a capacity that is close
to the capacity of a network with continuous weights at any coding level,
since the decrease in capacity is only about 1.5 compared to continuous
weights. Our results generalize the results obtained by Sompolinsky for a
coding level $f=0.5$ [11], to arbitrary coding levels. Furthermore, our
analysis shows that the storage capacity of CTF tends in the sparse coding
limit to the upper bound of storage capacity, in the space of all possible
binary connectivity matrices. We also provide a finite coding level correction
for this asymptotic solution, and the results indicate the capacities of TF
and CTF converge extremely slowly to the optimal capacity when the coding
level decreases, since the corrections are of order $1/\sqrt{\log(1/f)}$. In
particular, for $f=0.01$ [13], the capacity of the clipped model is only about
a third of the upper bound. Our results also show that sparse connectivity
matrices can allow these networks to have a larger information capacity per
synapse and thus encode information more efficiently.

The binary connectivity matrices used in this paper were constructed using a
clipped function whose argument is an analog variable containing information
about all stored patterns. This assumes that the synapse can store continuous
information during the learning phase, before binarizing this information. An
alternative scenario is that the synapse is required to be discrete during all
learning phases. Tsodyks, Amit and Fusi studied models with discrete synapses
under an online learning setting in which synapses only have information about
the currently shown pattern to make a transition between states. They showed
that this leads to a drastic decrease in storage capacity when the coding
level is $f=0.5$ [19, 20, 21], since in that case the total number of stored
patterns can scale at most as $\sqrt{N}$, implying a vanishing amount of
information stored per synapse in the large $N$ limit. Later work found that a
storage capacity of order 1 bit/synapse can be recovered in the sparse coding
limit ($f\sim\log(N)/N$), even when synapses are required to be discrete
during all phases of learning [21, 22].

Another scenario studied by multiple authors consists in synapses with binary
weights with multiple hidden states (describing e.g. different configurations
of protein interaction networks on the post-synaptic side) [25, 26, 24, 23].
With appropriate structure of hidden states, such synapses can greatly extend
the time for which synaptic connectivity can remain correlated with a pattern
shown at a particular time. This scenario has been primarily studied using a
signal-to-noise analysis quantifying the degree of correlation of the synaptic
matrix with patterns presented to the network. To our knowledge, this scenario
has never been implemented in attractor network models, and thus the storage
capacity in these multi-state models is still an open question. More
experimental data will be necessary to understand which class of models best
captures synaptic plasticity in neurobiological synapses.

## Appendix

### I. Calculation of $J$ and $\Delta_{0}$

To compute $J$ and $\Delta_{0}$ in Eq. (8) we use the same strategy as
Sompolinsky [10, 11]. We first calculate the average overlap between a given
pattern and the local field when the network is retrieving that pattern. Let
us denote $\tilde{\eta}_{i}^{\mu}=\eta_{i}^{\mu}-f$. For the learning rule
given by Eq. (7), we have:

$\begin{split}\langle\tilde{\eta}^{1}_{i}\sum_{j}\tilde{\eta}^{1}_{j}W_{ij}\rangle=Jf\left(1-f\right)\end{split}$
(25)

Similarly, the average overlap for clipped learning rule Eq. (4) is:

$\begin{split}\left\langle\tilde{\eta}^{1}_{i}\sum_{j}\tilde{\eta}^{1}_{j}W_{ij}\right\rangle&=\frac{\sqrt{p}}{N}\left\langle\sum_{j}\tilde{\eta}^{1}_{i}\tilde{\eta}^{1}_{j}F\left(\frac{\sum_{\mu}^{p}\tilde{\eta}^{\mu}_{i}\tilde{\eta}^{\mu}_{j}}{\sqrt{p}{f\left(1-f\right)}}\right)\right\rangle\\\
&=f\left(1-f\right)\left\langle xF\left(x\right)\right\rangle\\\ .\end{split}$
(26)

where

$\begin{split}x=\frac{\sum_{\mu}^{p}\tilde{\eta}^{\mu}_{i}\tilde{\eta}^{\mu}_{j}}{\sqrt{p}{f\left(1-f\right)}}.\end{split}$
(27)

In the large $p$ limit, $x$ becomes a random variable drawn from a standard
Gaussian distribution, $x\sim N\left(0,1\right)$ according to the Central
Limit Theorem. Thus, from Eqs. (25,26), we obtain the embedding strength $J$
as:

$J=\langle xF\left(x\right)\rangle,$ (28)

In order to obtain the variance $\delta_{ij}$ of Eq. (7), we calculate the
variance of synaptic weights for both linear and clipped learning rules. In
the linear case, we have:

$\begin{split}N^{2}\langle
W_{ij}^{2}\rangle&=N^{2}\frac{J^{2}}{f^{2}\left(1-f\right)^{2}N^{2}}\left\langle\left(\sum_{\mu}^{p}\tilde{\eta}_{i}^{\mu}\tilde{\eta}_{j}^{\mu}\right)^{2}\right\rangle\\\
&=pJ^{2}\left\langle\left(\sum_{\mu}^{p}\frac{\tilde{\eta}_{i}^{\mu}\tilde{\eta}_{j}^{\mu}}{f\left(1-f\right)\sqrt{p}}\right)^{2}\right\rangle\\\
&=N\alpha J^{2}.\end{split}$ (29)

The $<\cdot>$ in Eq. (29) goes to 1 when $p\to\infty$. For the clipped
$F\left(x\right)$, the variance of the weights is:

$\begin{split}N^{2}\langle
W_{ij}^{2}\rangle&=N^{2}\frac{1}{N^{2}}p\left\langle
F^{2}\left(\sum_{\mu}^{p}\frac{\tilde{\eta}_{i}^{\mu}\tilde{\eta}_{j}^{\mu}}{f\left(1-f\right)\sqrt{p}}\right)\right\rangle\\\
&=p\langle F^{2}\left(x\right)\rangle\\\ &=N\alpha\tilde{J}^{2},\end{split}$
(30)

where we denote $\langle F^{2}\left(x\right)\rangle$ as $\tilde{J}^{2}$. From
Eqs. (29,30), we can see the additional noise introduced by a nonlinear
$F\left(x\right)$ is:

$\langle\delta_{ij}^{2}\rangle=\frac{\alpha}{N}\left(\tilde{J}^{2}-J^{2}\right).$
(31)

Let $\Delta^{2}_{0}$ denote
${N\langle\delta_{ij}^{2}\rangle}/{J^{2}}{\alpha}$, we have:

$\Delta^{2}_{0}=\left(\frac{\tilde{J^{2}}}{J^{2}}-1\right).$ (32)

For $F\left(x\right)$ given by Eq. (5), we obtain the embedding strength and
noise parameter as

$J=1,\quad\Delta^{2}_{0}=\pi/2-1.$ (33)

### II. Storage capacity of the fully-connected CTF model

The Hamiltonian for the learning rule (7) is

$H=\frac{1}{2}\sum_{i\neq j}W_{ij}V_{i}V_{j}+\theta\sum_{i}V_{i},$ (34)

where $W_{ij}$ is determined from Eqs.(7,33). We can calculate the
corresponding free energy by using the replica method (see e.g. [3, 2]). To
compute the average logarithm of the partition function over the distribution
of all random binary patterns $\langle\log Z\rangle$ directly, we can use the
relation:

$\langle\log Z\rangle=\lim\limits_{n\to 0}\frac{\langle Z^{n}\rangle-1}{n}.$
(35)

For the Hamiltonian in equation (34), we have:

$\begin{split}\langle\langle
Z^{n}\rangle\rangle\propto&\Big{\langle}\Big{\langle}\mbox{Tr}_{V^{a}}\exp\left(\frac{\beta
J}{f(1-f)2N}\sum_{\mu
a}(\sum_{i}\tilde{\eta}_{i}^{\mu}V_{i}^{a})^{2}+\frac{\beta}{2}\sum_{ija}\delta_{ij}V_{i}^{a}V_{j}^{a}-\beta\theta\sum_{ia}V_{i}^{a}\right)\Big{\rangle}\Big{\rangle},\end{split}$
(36)

where $a=1,\ldots,n$ is the replica index and the double brackets mean the
average over both $V_{i}^{a}$ and $\eta_{i}^{\mu}$. We are interested in the
overlaps between the network state and the patterns stored in memory. We
assume that the network has a macroscopic overlap with a single stored pattern
($\mu=1$ the one currently being retrieved by the network) and define the
following order parameters:

$\begin{split}m_{a}&\equiv\frac{1}{N}\sum_{i}\tilde{\eta}_{i}^{1}V_{i}^{a},\\\
m_{\mu
a}&\equiv\frac{1}{\sqrt{N}}\sum_{i}\tilde{\eta}_{i}^{\mu}V_{i}^{a}.\end{split}$
(37)

The partition function can be rewritten in terms of these order parameters as

$\begin{split}\left\langle\left\langle
Z^{n}\right\rangle\right\rangle\propto&\Big{\langle}\Big{\langle}Tr_{V^{a}}\int\left(\prod_{\mu
a}dm_{a}^{\mu}\frac{\beta N}{2\pi}\right)\\\ \times&\exp\Big{(}-\frac{\beta
JN}{2f(1-f)}\sum_{a}m_{a}^{2}+\frac{\beta
J}{f(1-f)}\sum_{a}m_{a}\sum_{i}\tilde{\eta}_{i}^{1}V_{i}^{a}\\\ -&\frac{\beta
J}{2f(1-f)}\sum_{\mu a}\left(m_{a}^{\mu}\right)^{2}+\frac{\beta
J}{\sqrt{N}f(1-f)}\sum_{\mu
a}m_{a}^{\mu}\sum_{i}\tilde{\eta}_{i}^{\mu}V_{i}^{a}\\\
+&\frac{\beta}{2}\sum_{ija}V_{i}^{a}V_{j}^{a}\delta_{ij}-\beta\theta\sum_{ia}V_{i}^{a}\Big{)}\Big{\rangle}\Big{\rangle}.\end{split}$
(38)

The terms including $\eta_{i}^{\mu}$ and $\delta_{ij}$ in equation (38) can be
averaged,

$\begin{split}\Big{\langle}\exp\left(\frac{\beta J}{\sqrt{N}f(1-f)}\sum_{\mu
a}m_{a}^{\mu}\sum_{i}\tilde{\eta}_{i}^{\mu}V_{i}^{a}\right)\Big{\rangle}\
\propto\ \exp\left(\frac{\beta^{2}J^{2}}{2Nf(1-f)}\sum_{i\mu
ab}V_{i}^{a}V_{j}^{a}m_{a}^{\mu}m_{b}^{\mu}\right),\end{split}$ (39)
$\begin{split}\Big{\langle}\exp\left(\frac{\beta}{2}\sum_{ija}V_{i}^{a}V_{j}^{a}\delta_{ij}\right)\Big{\rangle}\
\propto\
\exp\left(N\beta^{2}J^{2}\Delta^{2}\sum_{ab}\left(\frac{1}{N}\sum_{i}V_{i}^{a}V_{i}^{b}\right)\left(\frac{1}{N}\sum_{i}V_{i}^{a}V_{i}^{b}\right)\right),\end{split}$
(40)

where $\Delta^{2}={N\langle\delta^{2}_{ij}\rangle}/{J^{2}}$. We then introduce
the order parameters:

$\begin{split}Q_{a}&=\frac{1}{N}\sum_{i}V_{i}^{a},\\\
q_{ab}&=\frac{1}{N}\sum_{i}V_{i}^{a}V_{i}^{b},\end{split}$ (41)

and use the integral representation of $\delta$ function :

$\begin{split}\delta\left(NQ_{a}-\sum_{i}V_{i}^{a}\right)&=\int\frac{dR_{a}}{2\pi}e^{-R_{a}\left(NQ_{a}-\sum_{i}V_{i}^{a}\right)},\\\
\delta\left(Nq_{ab}-\sum_{i}V_{i}^{a}V_{i}^{b}\right)&=\int\frac{dr_{ab}}{2\pi}e^{-r_{ab}\left(Nq_{ab}-\sum_{i}V_{i}^{a}V_{i}^{b}\right)}.\end{split}$
(42)

Combining Eqs. (38,39,40,41,42), the partition function can be written as:

$\begin{split}\langle\langle
Z^{n}\rangle\rangle\propto\int\left(\prod_{a}dQ_{a}dR_{a}\right)\left(\prod_{a\langle
b}q_{ab}r_{ab}\right)e^{-N\beta g\left(\beta,m,Q,q,R,r\right)},\end{split}$
(43)

where

$\begin{split}g=&\frac{J}{2f(1-f)}\sum_{a}m_{a}^{2}+\theta\sum_{a}Q_{a}-\frac{1}{\beta}\sum_{a}R_{a}Q_{a}-\beta
J^{2}\Delta^{2}\sum_{ab}q_{ab}^{2}-\frac{1}{\beta}\sum_{a\langle
b}r_{ab}q_{ab}\\\ &-\frac{1}{\beta}\log Tr_{V}^{a}\exp\left(\frac{\beta
J}{f(1-f)}\sum_{a}\tilde{\eta}^{1}m_{a}V^{a}-\sum_{a}R_{a}V^{a}-\sum_{a\langle
b}r_{ab}V^{a}V^{b}\right)\\\
&-\frac{\alpha}{\beta}\log\int\left(\prod_{a}dm_{a}\right)\exp\left(-\frac{\beta
J}{2f(1-f)}\sum_{a}m_{a}^{2}+\frac{\beta^{2}J^{2}}{2f(1-f)}\sum_{a\neq
b}m_{a}m_{b}q_{ab}\right).\end{split}$ (44)

The free energy per neuron is given as:

$\begin{split}G/N=\lim\limits_{n\to 0}\frac{1}{n}\mbox{min}\
g\left(\beta,m,Q,q,R,r\right).\end{split}$ (45)

In the large $N$ limit, min $g\left(\beta,m,Q,q,R,r\right)$ is dominated by
its value at saddle points. Next we will give the saddle point equations using
replica symmetric $ansatz$. We assume that saddle-point values of the order
parameters are not dependent on their replica index:

$\begin{split}&m_{a}=m,\\\ &Q_{a}=Q,R_{a}=R,\\\ &q_{ab}=q,r_{ab}=r\
\left(a\neq b\right).\end{split}$ (46)

Now the free energy per neuron is simplified to:

$\begin{split}G/N&=\frac{m^{2}}{2f(1-f)}+\frac{\alpha}{2\beta}\\{log\left(1-\beta
J\left(Q-q\right)\right)-\frac{\beta Jq}{1-\beta J\left(Q-q\right)}\\}\\\
&-\frac{rq}{2\beta}+\frac{RQ}{\beta}+\theta Q-\frac{1}{4}\beta
J^{2}\Delta^{2}\left(q^{2}-Q^{2}\right)\\\ &-\frac{1}{\beta}\int
Dz\log\left(1+\exp\left(\frac{\beta
J}{f(1-f)}m\tilde{\eta}+R-\frac{r}{2}+\sqrt{r}z\right)\right).\end{split}$
(47)

The saddle point equations are obtained by setting the derivatives of $G/N$ to
0:

$\displaystyle\frac{m}{f(1-f)}=\int
Dz\Big{\langle}\Big{\langle}\tilde{\eta}K\left(\beta\left(\frac{Jm\tilde{\eta}}{f(1-f)}+\frac{1}{\beta}\left(R-\frac{\alpha\beta^{2}\tilde{r}+\beta^{2}J^{2}\Delta^{2}q}{2}\right)+\sqrt{\alpha\tilde{r}+J^{2}\Delta^{2}q}z\right)\right)\Big{\rangle}\Big{\rangle},$
(48) $\displaystyle
R-\frac{\alpha\beta^{2}\tilde{r}+\beta^{2}J^{2}\Delta^{2}q}{2}=\frac{\alpha}{2}\frac{\beta(Q-q)J}{1-\beta
J\left(Q-q\right)}-\beta\theta+\frac{1}{2}\beta^{2}\left(Q-q\right)J^{2}\Delta^{2},$
$\displaystyle\tilde{r}=\frac{J^{2}q}{\left(1-J\beta\left(Q-q\right)\right)^{2}},$
$\displaystyle Q=\int
Dz\Big{\langle}\Big{\langle}K\left(\beta\left(\frac{Jm\tilde{\eta}}{f(1-f)}+\frac{1}{\beta}\left(R-\frac{\alpha\beta^{2}\tilde{r}+\beta^{2}q}{2}\right)+\sqrt{\alpha\tilde{r}+J^{2}\Delta^{2}q}z\right)\right)\Big{\rangle}\Big{\rangle},$
$\displaystyle q=\int
Dz\Big{\langle}\Big{\langle}K^{2}\left(\beta\left(\frac{Jm\tilde{\eta}}{f(1-f)}+\frac{1}{\beta}\left(R-\frac{\alpha\beta^{2}\tilde{r}+\beta^{2}q}{2}\right)+\sqrt{\alpha\tilde{r}+J^{2}\Delta^{2}q}z\right)\right)\Big{\rangle}\Big{\rangle},$

where
$\tilde{r}=\frac{1}{\beta^{2}\alpha}\left(r-\beta^{2}J^{2}\Delta^{2}q\right)$,
$K\left(x\right)=\left(1+\exp\left(-x\right)\right)^{-1}$, and
$Dz=dz\frac{\exp\left(-x^{2}/2\right)}{\sqrt{2\pi}}$. In the zero temperature
limit $\beta\to\infty$, these saddle points equations can be simplified to:

$\displaystyle\tilde{m}=\Phi(a_{1})-\Phi(a_{2}),$ (49)
$\displaystyle\tilde{r}=f\Phi(a_{1})+(1-f)\Phi(a_{2}),$ $\displaystyle
a_{1}=\frac{\tilde{\theta}-(1-f)\tilde{m}-Y}{\sqrt{\tilde{r}\alpha(1+\Delta_{0}^{2}(1-C)^{2})}},$
$\displaystyle
a_{2}=\frac{\tilde{\theta}+f\tilde{m}-Y}{\sqrt{\tilde{r}\alpha(1+\Delta_{0}^{2}(1-C)^{2})}},$
$\displaystyle Y=\frac{\alpha Cf}{2(1-C)}+\frac{1}{2}\alpha Cf\Delta_{0}^{2},$
$\displaystyle
C=\frac{f}{2\pi\alpha\tilde{r}}(fe^{-a_{1}^{2}/2}+(1-f)e^{-a_{2}^{2}/2})$

where
$\tilde{m}=m/f(1-f),\tilde{\theta}=\theta/f,\Delta_{0}^{2}=\Delta^{2}/\alpha$
and $\Phi\left(x\right)=\int_{x}^{\infty}Dz$.

### III. Sparse coding limit

In the sparse coding limit $f\to 0$, $C$ goes to zero and $Y$ goes to zero.
Eqs. (49) become:

$\displaystyle\tilde{m}=\Phi\left(\frac{\tilde{\theta}-\tilde{m}(1-f)}{\sqrt{\tilde{r}\alpha\left(1+\Delta^{2}_{0}\right)}}\right)-\Phi\left(\frac{\tilde{\theta}+f\tilde{m}}{\sqrt{\tilde{r}\alpha\left(1+\Delta^{2}_{0}\right)}}\right),$
(50)
$\displaystyle\tilde{r}=f\Phi\left(\frac{\tilde{\theta}-f\tilde{m}}{\sqrt{\tilde{r}\alpha\left(1+\Delta^{2}_{0}\right)}}\right)+\Phi\left(\frac{\tilde{\theta}}{\sqrt{\tilde{r}\alpha\left(1+\Delta^{2}_{0}\right)}}\right),$

In the small $f$ limit, $\tilde{m}\sim 1$ requires:

$\begin{split}&\Phi\left(\frac{\tilde{\theta}-\tilde{m}(1-f)}{\sqrt{\alpha
f\pi/2}}\right)\sim
1,\quad\Phi\left(\frac{\tilde{\theta}+\tilde{m}f}{\sqrt{\alpha
f\pi/2}}\right)\ll 1\end{split}$ (51)

Furthermore, $\alpha$ is maximized when $\tilde{r}$ is minimized, which
requires the stronger condition

$\begin{split}&\Phi\left(\frac{\tilde{\theta}+\tilde{m}f}{\sqrt{\alpha
f\pi/2}}\right)\ll f,\end{split}$ (52)

which leads to $\tilde{r}\sim f$. Using
$\lim\limits_{x\to+\infty}\Phi\left(x\right)\simeq\frac{1}{x\sqrt{2\pi}}\exp\left(\frac{-x^{2}}{2}\right)$,
in the small $f$ limit, Eq.(52) gives:

$\begin{split}\sqrt{\frac{\alpha
f}{2\tilde{\theta}^{2}}}\exp\left(-\frac{\tilde{\theta}^{2}}{\pi
f\alpha}\right)\ll f.\end{split}$ (53)

Rewriting $\alpha={k}/{f\log\left(f^{-1}\right)}$, we find that Eq. (53) is
satisfied provided $k<{\theta^{2}}/{\pi}$. Thus, the maximum storage capacity
$\alpha$ increase with $\tilde{\theta}^{2}$ as:

$\begin{split}\alpha_{c}\simeq\frac{\tilde{\theta}^{2}}{\pi f|\log
f|}.\end{split}$ (54)

In the sparse coding limit, the optimal threshold is obtained at
$\tilde{\theta}=1$ (the maximum value of $\tilde{\theta}$), and thus

$\begin{split}\alpha_{c}=\frac{1}{\pi f|\log f|}.\end{split}$ (55)

This storage capacity coincide with the optimal capacity obtained by Gutfreund
for the Ising interaction case (see Appendix IV for details).

We next ask the question of how close the threshold can be to 1, when the
coding level $f$ is small but finite. The threshold needs to be sufficiently
far from one, so that the argument of the function $\Phi$ in the first
condition in Eq. (51) is large and negative. We find that for thresholds that
are close to 1, the maximal capacity is

$\begin{split}\alpha=\frac{(1-\tilde{\theta})^{2}}{\pi
f\left|\log(1-\tilde{\theta})\right|)}.\end{split}$ (56)

Eqs.(54,56) give the optimal threshold $\theta_{opt}$ (i.e., the optimal value
of $\tilde{\theta}$) as a solution to the equation

$\begin{split}\frac{2\theta_{opt}^{2}\left|\log(1-\theta_{opt})\right|}{(1-\theta_{opt})^{2}}=\left|\log
f\right|,\end{split}$ (57)

The maximum storage capacity is

$\begin{split}\alpha_{c}\simeq\frac{\theta^{2}_{opt}}{\pi f|\log
f|}.\end{split}$ (58)

where $\theta_{opt}$ is given as a function of $f$ by Eq. (57). When $f\to 0$,
$\theta_{opt}\to 1$ and Eq.(58) becomes the Gutfreund bound Eq.(55). However,
this convergence is extremely slow, as shown in Fig. 2.

### IV. Bounds for capacity

Gardner Upper Bound for networks with Continuous weights (GUBC). This bound
was calculated by Gardner in 1987 for networks with continuous weights. The
upper bound is obtained when the volume of the space of solutions for the
weights $\\{J_{ij}\\}$ vanishes (see details in [4]). By solving the Eq.(37)
and Eq.(38) in [4], one can obtain the GUBC for arbitrary coding levels. This
result is shown by the red curve in Fig.2.

In the sparse coding limit, the asymptotic solution is given in Eq.(40) in
[4], where:

$\begin{split}\alpha_{cmax}=\frac{1}{2f|\mbox{ln}f|}.\end{split}$ (59)

Gutfreund and Stein Upper Bound for networks with Binary weights (GUBB). This
bound was calculated by Gutfreund and Stein in 1990 [14]. They extended
Gardner’s formalism to the case of networks with binary weights. Using a
replica symmetric ansatz, the solution space of binary weights vanishes when
capacity reaches:

$\begin{split}\alpha_{cmax}=\frac{2}{\pi}GUBC\end{split}$ (60)

In the sparse coding limit, Eq.(59) and Eq.(60) give:

$\begin{split}\alpha_{cmax}=\frac{1}{\pi f|\mbox{ln}f|}.\end{split}$ (61)

However, it can be shown that replica symmetry is broken, and Eq.(60) is an
overestimate. A better estimate of the upper bound for networks with binary
weights is given by zero entropy condition (see [14] for details), obtained by
solving Eqs.(20-24) and Eqs.(29-34) in [14]. This zero entropy line is shown
by the blue curve in Fig.2.

By numerically solving Eqs.(20-24, 29-34) in [14] and Eqs.(37-38) in [4], one
can see that the zero entropy line is getting close to the Gardner line
(Eq.(60)) when the coding level decrease. This numerical result indicates that
Eq.(61) is also a good estimate for the zero entropy line in the sparse coding
limit.

### V. Storage capacity of the sparsely-connected CTF model

Using similar calculations as in Appendix I, one can obtain that the nonlinear
learning rule Eq.(17) can be transformed to a linear learning rule:

$\displaystyle
W_{ij}=\frac{c_{ij}}{cNf(1-f)}\sum_{\mu}^{p}\left(\eta_{i}^{\mu}-f\right)\left(\eta_{j}^{\mu}-f\right)+\delta_{ij},$
(62)

where the variance of
$\delta_{ij}=\frac{\alpha\Delta_{0}^{2}}{N}=\frac{\alpha}{N}(1+\frac{\pi}{2})$,
and where the connection probability $c\ll 1$. In this case, the local field
$h_{i}$ can be written as:

$\displaystyle h_{i}$ $\displaystyle=\sum_{i\neq j}^{N}W_{ij}V_{j}$ (63)
$\displaystyle=\frac{1}{cf(1-f)N}\tilde{\eta}_{i}^{1}\sum_{j}^{N}c_{ij}\tilde{\eta}_{j}^{1}V_{j}$
$\displaystyle+\frac{1}{cf(1-f)N}\sum_{\mu>1}\sum_{j}^{N}c_{ij}\tilde{\eta}_{i}^{\mu}\tilde{\eta}_{j}^{\mu}V_{j}$
$\displaystyle+\sum_{j}^{N}\delta_{ij}V_{j},$

where $\tilde{\eta}_{i}^{1}$ denote the pattern that is currently being
retrieved by the network. When $N$ and $p$ are large, the second and third
term in the r.h.s. of Eq. (63) $Y_{1}$ and follow a Gaussian distribution with
zero mean. Introducing order parameters $\tilde{m}=m/f(1-f)$ and $q$ defined
in Eq.(19), Eq.(63) can be simplified to:

$\displaystyle h_{i}$ $\displaystyle=\tilde{\eta}_{i}^{1}\tilde{m}+Y,$ (64)

where the variance of $Y$ is:

$\displaystyle var(Y)$ $\displaystyle=\alpha_{s}(\Delta_{0}^{2}+1)q$ (65)

and $\alpha_{s}=p/cN$.

$\displaystyle\tilde{m}$
$\displaystyle=\frac{1}{Nf(1-f)}\sum_{j}^{N}\tilde{\eta}_{j}^{1}\Theta(\tilde{\eta}_{j}^{1}\tilde{m}+Y-\theta),$
(66)

$\displaystyle q$
$\displaystyle=\frac{1}{N}\sum_{j}^{N}\Theta^{2}(\tilde{\eta}_{j}^{1}\tilde{m}+Y-\theta).$
(67)

Averaging Eq.(66,67) over $\eta$ and the Gaussian noise $Y$, the mean-field
equations for order parameters $\tilde{m}$ and $q$ are

$\begin{split}&\tilde{m}=\Phi\left(\frac{\tilde{\theta}-\tilde{m}(1-f)}{\sqrt{\alpha_{s}q(1+\Delta_{0}^{2})}}\right)-\Phi\left(\frac{\tilde{\theta}+\tilde{m}f}{\sqrt{\alpha_{s}q(1+\Delta_{0}^{2})}}\right),\\\
&q=f\Phi\left(\frac{\tilde{\theta}-\tilde{m}(1-f)}{\sqrt{\alpha_{s}q(1+\Delta_{0}^{2})}}\right)+(1-f)\Phi\left(\frac{\tilde{\theta}+\tilde{m}f}{\sqrt{\alpha_{s}q(1+\Delta_{0}^{2})}}\right),\end{split}$
(68)

Note that these equations coincide with the fully connected case in the sparse
coding limit, as in the case of continuous weights [2, 17].

### VI. Numerical simulations

The simulations in Fig.1 used a network with 4,000 neurons and coding level
$f=0.02$. The overlaps $\tilde{m}$ are averaged over five independent
realizations. Simulations consist in a learning phase in which the
connectivity matrix $W_{ij}$ is built, and a retrieval phase in which network
dynamics run until it reaches a fixed point. For each input pattern $\mu$, we
choose as initial conditions $\\{V_{i}=\eta_{i}^{\mu}\\}$. Overlaps $m$ are
obtained by averaging over all $m^{\mu}$.

## References

  * [1] Hopfield, John J. ”Neural networks and physical systems with emergent collective computational abilities.” Proceedings of the national academy of sciences 79.8 (1982): 2554-2558.
  * [2] Tsodyks, Mikhail V., and Mikhail V. Feigel’man. ”The enhanced storage capacity in neural networks with low activity level.” EPL (Europhysics Letters ) 6.2 (1988): 101.
  * [3] Amit, Daniel J., Hanoch Gutfreund, and Haim Sompolinsky. ”Statistical mechanics of neural networks near saturation.” Annals of physics 173.1 (1987): 30-67.
  * [4] Gardner, Elizabeth. ”The space of interactions in neural network models.” Journal of physics A: Mathematical and general 21.1 (1988): 257.
  * [5] Petersen, Carl CH, et al. ”All-or-none potentiation at CA3-CA1 synapses.” Proceedings of the National Academy of Sciences 95.8 (1998): 4732-4737.
  * [6] O’Connor, Daniel H., Gayle M. Wittenberg, and Samuel S-H. Wang. ”Graded bidirectional synaptic plasticity is composed of switch-like unitary events.” Proceedings of the National Academy of Sciences 102.27 (2005): 9679-9684.
  * [7] Hruska, Martin, et al. ”Synaptic nanomodules underlie the organization and plasticity of spine synapses.” Nature neuroscience 21.5 (2018): 671-682.
  * [8] Dorkenwald, Sven, et al. ”Binary and analog variation of synapses between cortical pyramidal neurons.” BioRxiv (2019).
  * [9] Krauth, Werner, and Marc Mézard. ”Storage capacity of memory networks with binary couplings.” Journal de Physique 50.20 (1989): 3057-3066.
  * [10] Sompolinsky, Haim. ”The theory of neural networks: The Hebb rule and beyond.” Heidelberg colloquium on glassy dynamics. Springer, Berlin, Heidelberg, 1987.
  * [11] Sompolinsky, Haim. ”Neural networks with nonlinear synapses and a static noise.” Physical Review A 34.3 (1986): 2571.
  * [12] Lee, Jae Sung, et al. ”The statistical structure of the hippocampal code for space as a function of time, context, and value.” Cell 183.3 (2020): 620-635.
  * [13] Waydo, Stephen, et al. ”Sparse representation in the human medial temporal lobe.” Journal of Neuroscience 26.40 (2006): 10232-10234.
  * [14] Gutfreund, Hanoch, and Yaakov Stein. ”Capacity of neural networks with discrete synaptic couplings.” Journal of Physics A: Mathematical and General 23.12 (1990): 2613.
  * [15] Mézard, M., Parisi, G., and Virasoro, M. A. (1987). Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications (Vol. 9). World Scientific Publishing Company.
  * [16] Markram, Henry, et al. ”Physiology and anatomy of synaptic connections between thick tufted pyramidal neurones in the developing rat neocortex.” The Journal of physiology 500.2 (1997): 409-440.
  * [17] Tsodyks, M. V. ”Associative memory in asymmetric diluted network with low level of activity.” EPL (Europhysics Letters) 7.3 (1988): 203.
  * [18] Derrida, Bernard, Elizabeth Gardner, and Anne Zippelius. ”An exactly solvable asymmetric neural network model.” EPL (Europhysics Letters) 4.2 (1987): 167.
  * [19] Tsodyks, M. V. ”Associative memory in neural networks with binary synapses” Modern Physics Letters B 4 (11), 713-71613..
  * [20] Amit, Daniel J., and Stefano Fusi. ”Constraints on learning in dynamic synapses.” Network: Computation in Neural Systems 3.4 (1992): 443-464.
  * [21] Amit, Daniel J., and Stefano Fusi. ”Learning in neural networks with material synapses.” Neural Computation 6.5 (1994): 957-982.
  * [22] Dubreuil, Alexis M., Yali Amit, and Nicolas Brunel. ”Memory capacity of networks with stochastic binary synapses.” PLoS computational biology 10.8 (2014): e1003727.
  * [23] Benna, Marcus K., and Stefano Fusi. ”Computational principles of synaptic memory consolidation.” Nature neuroscience 19.12 (2016): 1697-1706.
  * [24] Lahiri, Subhaneil, and Surya Ganguli. ”A memory frontier for complex synapses.” Advances in neural information processing systems 26 (2013): 1034-1042.
  * [25] Fusi, Stefano, Patrick J. Drew, and Larry F. Abbott. ”Cascade models of synaptically stored memories.” Neuron 45.4 (2005): 599-611.
  * [26] Fusi S, Abbott LF. Limits on the memory storage capacity of bounded synapses. Nat Neurosci. 2007 Apr;10(4):485-93