text:
# The Dynamical Systems Approach to theEquations of a Linearly ViscousCompressible Barotropic Fluid

## 1. Introduction

The long-time behaviour of solutions to the evolutionary equations arising in the mathematical fluid mechanics has been the subject of many theoretical studies. This type of problems is apparently related to the phenomena of turbulence, and there is still a significant gap between many formal “scenarios” and mathematically rigorous results.

The dynamical systems in question are usually related to a system of partial differential equations and, consequently, they are defined in an infinite dimensional phase space. On the other hand, the presence of dissipative terms in the equations due to viscosity results in the existence of global attractors-compact invariant sets attracting uniformly in time all trajectories emanating from a given bounded set. Such a theory is well developed for the incompressible linearly viscous fluids, and the reader may consult the monographs of BABIN and VISHIK , TEMAM  or CONSTANTIN et al.  for the recent state of art.

On the other hand, much less seems to be known for the compressible fluids. While there is an existence theory of the weak solutions for the incompressible Navier-Stokes equations due to LERAY , its “compressible” counterpart appeared only recently in the work of LIONS . Even in the incompressible case, there is a qualitative difference between the two-dimensional case solved by LADYZHENSKAYA , and the three-dimensional case representing one of the most challenging unsolved problems of modern mathematics. It is worth-noting that a similar gap divides the one and more-dimensional problems for the compressible fluids.

The time evolution of the fluid density $`\varrho =\varrho (t,x)`$ and the velocity $`\stackrel{}{u}=\stackrel{}{u}(t,x)`$ is governed by the Navier-Stokes system of equations:

$$_t\varrho +\mathrm{div}(\varrho \stackrel{}{u})=0,$$
(1.1)
$$_t(\varrho \stackrel{}{u})+\mathrm{div}(\varrho \stackrel{}{u}\stackrel{}{u})+p=\mathrm{div}S+\varrho \stackrel{}{f},$$
(1.2)
where $`p`$ is the pressure, $`S`$ the viscous stress tensor, and $`\stackrel{}{f}`$ a given external force.

In what follows, we consider linearly viscous (Newtonian) fluids where the viscous stress is related to the velocity by the constitutive law

$$S=\mu \left(\stackrel{}{u}+\stackrel{}{u}^T\right)+\lambda \mathrm{div}\stackrel{}{u}I,$$
(1.3)
where the viscosity coefficients satisfy

$$\mu >0,\lambda +\mu 0.$$
(1.4)
Generally speaking, the pressure $`p`$ depends on the density and the internal energy (temperature) of the fluid. If it is the case, the system (1.1), (1.2) is not closed and should be complemented by the energy equation. Unfortunately, however, the available global existence results for this full system allow for only for small initial data (cf. MATSUMURA and NISHIDA ).

There are physically relevant situations when one can assume the flow is barotropic, i.e., the pressure depends solely on the density. This is the case when either the temperature (the isothermal case) or the entropy (the isentropic case) are supposed to be constant. The typical constitutive relation between the pressure and the density then reads

$$p=p(\varrho )=a\varrho ^\gamma ,a>0,$$
(1.5)
where $`\gamma =1`$ in the isothermal case, and $`\gamma >1`$ represents the adiabatic constant in the isentropic regime. More general and even non-monotone pressure-density constitutive laws arise in nuclear plasma physics (see ). In the barotropic regime, the equations (1.1), (1.2) form a closed system and complemented by suitable initial and boundary conditions represent a (at least formally) well-posed problem.

If the problem is posed on a spatial domain $`\mathrm{\Omega }R^N`$, one usually assumes that the fluid adheres completely to the boundary which is mathematically expressed by the no-slip boundary conditions for the velocity:

$$\stackrel{}{u}|_\mathrm{\Omega }=0.$$
(1.6)
Note that for viscous fluids such a condition is in a very good agreement with physical experiment.

In accordance with the deterministic principle, the state of the system at any time $`t>t_0`$ should be given by the initial conditions

$$\varrho (t_0)=\varrho _I,(\varrho \stackrel{}{u})(t_0)=\stackrel{}{q}_I.$$
(1.7)
The function $`\varrho _I`$ is non-negative and the momentum $`\stackrel{}{q}_I`$ satisfies the compatibility condition

$$\stackrel{}{q}_I=0\text{a.a. on the set}\{\varrho _I=0\}.$$
(1.8)
The reason why we impose the initial conditions for the momentum $`\varrho \stackrel{}{u}`$ rather than for the velocity $`\stackrel{}{u}`$ is that the former quantity is weakly continuous with respect to time while the instantaneous values of the velocity are determined only almost anywhere with respect to time. Clearly such a problem does not arise provided the initial density $`\varrho _I`$ is strictly positive. However, it is an interesting open problem whether or not this property is preserved at any positive time for any distributional solution of the problem satisfying the natural energy estimates (cf. HOFF and SMOLLER ).

## 2. Finite energy weak solutions and well-posedness

In order to study the long-time behaviour, one should first make sure that the class of objects one deals with is not void. More precisely, one should be able to prove the existence of global-in-time solution for any initial data $`\varrho _I`$, $`\stackrel{}{q}_I`$ satisfying some physically relevant hypothesis.

Multiplying the equations of motion by $`\stackrel{}{u}`$ and integrating by parts one deduces the energy inequality

$$\frac{\mathrm{d}}{\mathrm{d}t}E[\varrho ,\stackrel{}{u}]+_\mathrm{\Omega }\mu |\stackrel{}{u}|^2+(\lambda +\mu )|\mathrm{div}\stackrel{}{u}|^2\mathrm{d}x_\mathrm{\Omega }\varrho \stackrel{}{f}\stackrel{}{u}dx,$$
(2.1)
where the total energy $`E`$ is given by the formula

$$E[\varrho ,\stackrel{}{u}]=_\mathrm{\Omega }\varrho |\stackrel{}{u}|^2+P(\varrho )\mathrm{d}x$$

with

$$P^{}(\varrho )\varrho P(\varrho )=p(\varrho ).$$

Note that in the most common isentropic case, the function $`P`$ can be taken in the form

$$P(\varrho )=\frac{a}{\gamma 1}\varrho ^\gamma .$$

The energy inequality is the main (and almost the only one) source of a priori estimates. Accordingly, “reasonble” solutions of the problem (1.1), (1.2) defined on a bounded time interval $`IR`$ should belong to the class

$$\varrho 0,\varrho L^{\mathrm{}}(I;L^\gamma (\mathrm{\Omega })),\stackrel{}{u}L^2(I;W_0^{1,2}(\mathrm{\Omega },R^N)).$$
(2.2)
The energy inequality (2.1) represents an additional constraint imposed on any solution $`\varrho `$, $`\stackrel{}{u}`$ of the problem. Similarly as in the theory of the variational (weak) solutions developed for the incompressible case by Leray, it is not clear if it is satisfied for any weak solution of the problem.

Following DiPERNA and LIONS  we shall say that $`\varrho `$, $`\stackrel{}{u}`$ is a renormalized solution of the continuity equation (1.1) if the identity

$$_tb(\varrho )+\mathrm{div}(b(\varrho )\stackrel{}{u})+\left(b^{}(\varrho )\varrho b(\varrho )\right)\mathrm{div}\stackrel{}{u}=0$$
(2.3)
holds in the sense of distributions for any function $`bC^1(R)`$ such that

$$b^{}(\varrho )=0\text{for all}\varrho \mathrm{const}(b).$$
(2.4)
Similarly as for the energy inequality, it is not known if any weak solution $`\varrho `$, $`\stackrel{}{u}`$ of (1.1) satisfies (2.3).

Definition 2.1 We shall say that $`\varrho `$, $`\stackrel{}{u}`$ is a finite energy weak solution of the problem (1.1 \- 1.6) on a set $`I\times \mathrm{\Omega }`$ if the following conditions are satisfied:

* The functions $`\varrho `$, $`\stackrel{}{u}`$ belong to the function spaces determined in (2.2);
* the energy inequality (2.1) is satified in $`𝒟^{}(I)`$;
* the continuity equation (1.1) as well as its renormalized version (2.3) hold in $`𝒟^{}(I\times R^N)`$ provided $`\varrho `$, $`\stackrel{}{u}`$ were extended to be zero outside $`\mathrm{\Omega }`$;
* the momentum equation (1.2) is satisfied in $`𝒟^{}(I\times \mathrm{\Omega })`$.

The most general available existence result reads as follows:

Theorem 2.1 Let $`\mathrm{\Omega }R^N`$, $`N=2,3`$ be a bounded Lipschitz domain. Let $`I=(0,T)`$, and let the initial data $`\varrho _I`$, $`\stackrel{}{q}_I`$ satisfy (1.8) together with

$$\varrho _IL^\gamma (\mathrm{\Omega }),\frac{|\stackrel{}{q}_I|^2}{\varrho _I}L^1(\mathrm{\Omega }).$$

Let $`\stackrel{}{f}`$ be a bounded measurable function of $`tI`$, $`x\mathrm{\Omega }`$. Finally, let the pressure $`pC[0,\mathrm{})C^1(0,\mathrm{})`$ be given by a constitutive law

$$p=p(\varrho ),\frac{1}{a}\varrho ^{\gamma 1}bp^{}(\varrho )a\varrho ^{\gamma 1}+b,\text{for all}\varrho >0,$$

where $`a>0`$, $`b0`$, and

$$\gamma >\frac{N}{2}.$$

Then the problem (1.1 \- 1.6) admits a finite energy weak solution $`\varrho `$, $`\stackrel{}{u}`$ on $`I\times \mathrm{\Omega }`$ satisfying the initial conditions (1.7).

In his pioneering work, LIONS  proved Theorem 2.1 for $`\mathrm{\Omega }`$ regular, $`p`$ monotone, and $`\gamma \frac{3}{2}`$ if $`N=2`$, and $`\gamma \frac{9}{5}`$ for $`N=3`$. The hypotheses concerning $`\gamma `$ were relaxed in , , the case of a general bounded domain $`\mathrm{\Omega }`$ treated in , and the hypothesis of monotonicity of the pressure removed in , .

## 3. Ultimate boundedness

The first issue to be discussed when describing the long-time asymptotics of a given dynamical system is ultimate boundedness or dissipativity. This means there exists an absorbing set bounded in a suitable topology. Here “suitable topology” is of course that one induced by the total energy $`E`$. One of possible results in this direction is contained in the following theorem.

Theorem 3.1 Let $`\mathrm{\Omega }R^N`$, $`N=2,3`$ be a bounded Lipschitz domain. Let $`𝐟`$ be a bounded measurable function such that

$$\mathrm{ess}\underset{tR,x\mathrm{\Omega }}{sup}|\stackrel{}{f}(t,x)|F.$$

Let the pressure $`p`$ be given by the isentropic constitutive law

$$p=a\varrho ^\gamma \text{with}\gamma >1\text{if}N=2,\gamma >\frac{5}{3}\text{for}N=3.$$

Finally, set

$$_\mathrm{\Omega }\varrho dx=m>0.$$

Then there exists a constant $`E_{\mathrm{}}`$ depending solely on $`m`$ and $`F`$ having the following property:

Given $`E_I`$, there exists a time $`T=T(E_I)`$ such that

$$E[\varrho ,\stackrel{}{u}](t)E_{\mathrm{}}\text{for a.a.}t>T$$

provided

$$\mathrm{ess}\underset{t0+}{lim\; sup}E[\varrho ,\stackrel{}{u}](t)E_I$$

and $`\varrho `$, $`\stackrel{}{u}`$ is a finite energy weak solution of the problem (1.1 \- 1.6) on $`(0,\mathrm{})\times \mathrm{\Omega }`$.

The proof of Theorem 3.1 can be found in  and . The reader will have noticed the “critical” exponent $`\gamma >\frac{5}{3}`$ which is larger than in the existence Theorem 2.1 and, as a matter of fact, does not include any physically relevant case. Indeed the value $`\frac{5}{3}`$ happens to be the largest adiabatic constant corresponding to a monoatomic gas.

Under the hypotheses of Theorem 3.1, the dissipative mechanism induced by viscocity is strong enough to prevent any “resonance” phenomena, i.e., the existence of unbounded solutions driven by a bounded external force. Another interesting feature is the existence of periodic solutions provided $`\stackrel{}{f}`$ is periodic in time.

Theorem 3.2 In addition to the hypotheses of Theorem 3.1, assume that $`\stackrel{}{f}`$ is time periodic, i.e.,

$$\stackrel{}{f}(t+\omega ,x)=\stackrel{}{f}(t,x)\text{for a.a.}tR,x\mathrm{\Omega }$$

with a certain period $`\omega >0`$.

Then there exists at least one finite energy weak solution of the problem (1.1 \- 1.6) on $`R\times \mathrm{\Omega }`$ which is $`\omega `$periodic in time, i.e.,

$$\varrho (t+\omega )=\varrho (t),(\varrho \stackrel{}{u})(t+\omega )=(\varrho \stackrel{}{u})(t)\text{for all}tR.$$

See  for the proof.

## 4. Asymptotic compactness

The property of asymptotic compactness of a given dynamical system plays a crucial role in the proof of existence of a global (compact) attractor. Here, the main problem is the density component which is bounded only in $`L^\gamma (\mathrm{\Omega })`$, and, consequently, compact only with respect to the weak topology on this space. Moreover, given the hyperbolic character of the continuity equation, one cannot hope any possible oscillations of the density to be killed at a finite time. However, the amplitude of possible oscillations is decreasing in time uniformly on trajectories emanating from a given bounded set.

Theorem 4.1 Let $`\mathrm{\Omega }R^N`$, $`N=2,3`$ be a bounded Lipschitz domain. Let the pressure $`p`$ be given by the isentropic constitutive relation

$$p=a\varrho ^\gamma ,a>0,\gamma >1\text{if}N=2,\gamma >\frac{5}{3}\text{for}N=3.$$

Let $`\stackrel{}{f}_n`$ be a sequence of functions such that

$$\mathrm{ess}\underset{tR,x\mathrm{\Omega }}{sup}|\stackrel{}{f}_n(t,x)|F\text{independently of}n=1,2,\mathrm{}.$$

Finally, let $`\varrho _n`$, $`\stackrel{}{u}_n`$ be a sequence of finite energy weak solutions to the problem (1.1 \- 1.6) on $`(0,\mathrm{})\times \mathrm{\Omega }`$ such that

$$_\mathrm{\Omega }\varrho _ndx=m,$$

$$\mathrm{ess}\underset{t0+}{lim\; sup}E[\varrho _n,\stackrel{}{u}_n](t)E_I$$

independently of $`n=1,2,\mathrm{}.`$

Then any sequence of times $`t_n\mathrm{}`$ contains a subsequence such that

$$\varrho _n(t_n+t)\varrho (t)\text{strongly in}L^1(\mathrm{\Omega })\text{and weakly in}L^\gamma (\mathrm{\Omega }),$$

$$(\varrho _n\stackrel{}{u}_n)(t_n+t)(\varrho \stackrel{}{u})(t)\text{weakly in}L^1(\mathrm{\Omega },R^N)$$

for any $`tR`$.

Moreover,

$$_J_\mathrm{\Omega }|\varrho _n(t_n+t,x)\varrho (t,x)|^\gamma dx𝑑t0,$$

$$_J_\mathrm{\Omega }|(\varrho _n\stackrel{}{u}_n)(t_n+t,x)(\varrho \stackrel{}{u})(t,x)|dx𝑑t0$$

for any bounded interval $`JR`$. The limit functions $`\varrho `$, $`\stackrel{}{u}`$ represent a globally defined (for $`tR`$) finite energy weak solution of the problem (1.1 \- 1.6) driven by a force

$$\stackrel{}{f}=\underset{n\mathrm{}}{lim}\stackrel{}{f}_n(t_n+)\text{in the weak star topology of the space}L^{\mathrm{}}(R\times \mathrm{\Omega },R^N),$$

and such that

$$\mathrm{ess}\underset{tR}{sup}E[\varrho ,\stackrel{}{u}](t)<\mathrm{}.$$

The proof of Theorem 4.1 is given in  (see also ).

## 5. The long-time behaviour, attractors

The results presented in Sections 3., 4. allow us to develop a theory of attractors analogous to that one for the incompressible flows (see e.g. TEMAM ). Assume, for the sake of simplicity, that the driving force $`\stackrel{}{f}`$ is independent of $`t`$. We introduce

$$𝒜=\{[\varrho _I,\stackrel{}{q}_I]|\varrho _I=\varrho (0),\stackrel{}{q}_I=(\varrho \stackrel{}{u})(0)\text{where}\varrho ,\stackrel{}{u}\text{is a finite energy weak solution}$$

$$\text{of the problem (}\text{1.1}\text{ - }\text{1.6}\text{) on}R\times \mathrm{\Omega }\text{with}E[\varrho ,\stackrel{}{u}]L^{\mathrm{}}(R)\}.$$
(5.1)
In other words, the set $`𝒜`$ is formed by all globally defined (for $`tR`$) trajectories whose energy is uniformly bounded.

The next result shows that $`𝒜`$ is a global attractor in the sense of FOIAS and TEMAM .

Theorem 5.1 Let $`\mathrm{\Omega }R^N`$, $`N=2,3`$ be a bounded Lipschitz domain. Let $`\stackrel{}{f}=\stackrel{}{f}(x)`$ be a bounded measurable function independent of time, and let the pressure $`p`$ be given by the isentropic constitutive law

$$p=a\varrho ^\gamma ,a>0,\gamma >1\text{for}N=2,\gamma >\frac{5}{3}\text{if}N=3.$$

Then the set $`𝒜`$ defined by (5.1) is compact in the space

$$L^\alpha (\mathrm{\Omega })\times L_{weak}^{\frac{2\gamma }{\gamma +1}}(\mathrm{\Omega }),$$

and

$$\underset{[\varrho ,\stackrel{}{u}](E_I)}{sup}\left[\underset{[\varrho _I,\stackrel{}{q}_I]𝒜}{inf}\left(\varrho (t)\varrho _I_{L^\alpha (\mathrm{\Omega })}+\left|_\mathrm{\Omega }((\varrho \stackrel{}{u})(t)\stackrel{}{q}_I)\varphi dx\right|\right)\right]0$$

$$\text{for}t\mathrm{}$$

for any $`1\alpha <\gamma `$ and any $`\varphi L^{\frac{2\gamma }{\gamma 1}}(\mathrm{\Omega },R^N)`$, where the symbol $`(E_I)`$ stands for the set of all finite energy weak solutions of the problem (1.1 \- 1.6) on $`(0,\mathrm{})\times \mathrm{\Omega }`$ such that

$$\mathrm{ess}\underset{t0}{lim\; sup}E[\varrho ,\stackrel{}{u}](t)E_I.$$

The proof can be found in .

To conclude, we give another result concerning the long-time behaviour of solutions on condition that the driving force is a gradient of a scalar potential.

Theorem 5.2 Let $`\mathrm{\Omega }R^N`$, $`N=2,3`$ be a bounded Lipschitz domain. Let the pressure $`p`$ be given by the isentropic constitutive relation

$$p=a\varrho ^\gamma ,a>0,\gamma >\frac{N}{2}.$$

Let the driving force $`\stackrel{}{f}`$ be of the form

$$\stackrel{}{f}=F,$$

where $`F=F(x)`$ be a scalar potential which is globally Lipschitz on $`\overline{\mathrm{\Omega }}`$ and such that the upper level sets

$$[F>k]=\{x\mathrm{\Omega }|F(x)>k\}$$

are connected for any $`kR`$.

Then any finite energy weak solution $`\varrho `$, $`\stackrel{}{u}`$ of the problem (1.1 \- 1.6) satisfies

$$\varrho (t)\varrho _s\text{in}L^\gamma (\mathrm{\Omega }),(\varrho \stackrel{}{u})(t)0\text{in}L^1(\mathrm{\Omega })\text{as}t\mathrm{},$$

where $`\varrho _s`$ is a solution of the static problem

$$a\varrho _s^\gamma =\varrho _s\stackrel{}{f}\text{on}\mathrm{\Omega }.$$

Acknowledgement  The work was supported by Grant No. 201/02/0854 of GA ČR.


format: text

text:
# Analysis of the intraspinal calcium dynamics and its implications on the plasticity of spiking neurons

## I Introduction

Calcium ions are ubiquitous mediators of metabolic change in many cellular systems. In the cortex, it is believed to be the primary chemical signal for the induction of bidirectional synaptic plasticity , which is widely accepted to be the basis for many forms of learning, memory and development.

Experimentally, bi-directional synaptic plasticity can be induced by various different induction protocols , including the recently characterized spike-time-dependent plasticity (STDP) . In this induction method, if a presynaptic action potential occurs within a window of a few tens milliseconds before a postsynaptic action potential ($`\mathrm{\Delta }t>0`$), long-term potentiation (LTP) is elicited; if the order is switched ($`\mathrm{\Delta }t<0`$), long-term depression (LTD) happens. Interestingly, STDP and many other induction protocols share a common property, which is their dependence on the amount of integrated activation of the N-methyl-D-aspartate receptors (NMDAR) during conditioning. NMDAR are natural coincidence detectors; their activation relies both on the efficiency of neurotransmitter release from the presynaptic cell, and on the level of depolarization of the postsynaptic cell. Most remarkably, they are permeable to calcium ions, which suggests a functional link between calcium influx through NMDAR and the induction of LTP and LTD.

Recently, several dynamical models have been proposed to account for STDP , the relationship between STDP and rate-based models has been examined using averaging methods , and there have been attempts to account for receptive field formation on the basis of STDP . We have recently proposed a unified theory of synaptic plasticity  based on the following assumptions: (1) the level of intracellular calcium concentration determines the sign and magnitude of synaptic plasticity : when calcium falls below a threshold $`\theta _d`$, no plasticity occurs, when it falls between $`\theta _d`$ and $`\theta _p`$, $`\theta _p>\theta _d`$, LTD is induced, and for calcium above $`\theta _p`$, LTP is induced; (2) the relevant sources of calcium are the NMDA channels; and (3) dendritic back-propagating action potentials (BPAPs) contributing to STDP have a slow ”after-depolarizing” tail component. We have shown that these simple assumptions are sufficient to account for the various experimental plasticity-induction protocols. In addition, this model has also produced previously uncharacterized predictions, such as: (a) the shape of the STDP learning curve should be frequency-dependent; and (b) there should be a novel form of spike time-dependent LTD for $`\mathrm{\Delta }t`$ larger than the LTP-inducing intervals. Recent experimental results are consistent with these predictions .

In this paper, we present a systematic calculation of the calcium transients as a function of pre- and post-spike timings, as well as of the neuronal firing rates. In section II we derive the solution for the mean calcium transients, in the simplest case where the BPAP consists of a single exponential, and show that its dependence on $`\mathrm{\Delta }t`$ is rather intuitive. A more realistic approach is taken in section III, where the BPAP is composed of a sum of two exponentials with different time constants. This more complex assumption is in closer agreement with our previous work  and is consistent with experimental observations . The analysis shows that the slow component contributes for different calcium levels in the baseline condition (pre only) and the post-pre condition ($`\mathrm{\Delta }t<0`$), while the fast component sharpens the difference of the calcium levels in the transition between post/pre ($`\mathrm{\Delta }t<0`$) and pre/post ($`\mathrm{\Delta }t>0`$) conditions. The sharpness of this transition is shown to be further enhanced by the appropriate choice of the calcium decay time-constant. For completeness, in section IV, we study the effects of calcium-accumulation due to repetitive firing which affects the rate-dependence of some plasticity-induction protocols including the STDP learning window. Finally, in section V we take into account the stochastic properties of synaptic transmission and estimate the trial-by-trial variability of the calcium transients. This variability can be significant if the number of postsynaptic NMDAR in each spine is small. Our analysis demonstrates that this variability depends $`\mathrm{\Delta }t`$; increasing with $`\mathrm{\Delta }t`$ for $`\mathrm{\Delta }t>0`$. These results could significantly alter the form of the STDP curve under different experimental conditions.

## II Analysis of calcium influx through NMDAR

The postsynaptic intraspinal calcium ions are removed from these compartments both through diffusion and through binding to calcium buffers. Let these phenomena be characterized by a single time constant $`\tau `$, and assume that the calcium dynamics follow a first-order linear differential equation,

$$\frac{d[Ca]}{dt}=I(t)\frac{1}{\tau }[Ca].$$
(1)
where $`[Ca]`$ is the intracellular calcium concentration, $`I(t)`$ is the inward calcium current, and $`\tau `$ is the time constant of the passive decay process. This simple dynamics is to the first approximation consistent with experiments  . Given a initial condition $`[Ca](0)`$, equation 1 has a solution of the form:

$`[Ca](t)=e^{t/\tau }\left[{\displaystyle _0^t}e^{t/\tau }I(t^{})𝑑t^{}+[Ca](0)\right].`$ (2)
Notice that if the calcium current is a superposition of separate current sources, the calcium concentration can also be decomposed into such a sum.

In this paper we consider the NMDA channels as the sole source of calcium currents. These channels are voltage dependent , and have a long open time that can usually be described as a sum of two exponentials , thus its currents can be generally described as:

$$I(t,\{t^{pre}\},\{t^{post}\})=\overline{g}(V)f,$$
(3)
where $`\{t^{pre}\}`$ and $`\{t^{post}\}`$ are the sets of times of pre- and postsynaptic action-potentials (APs) that occured prior to time $`t`$ and $`\overline{g}`$ is the mean total conductance of the population of NMDAR in a synapse. The $``$-function represents the voltage dependence of calcium influx through NMDAR, and $`f`$ is the fraction of NMDAR in open state. It is easy to see that $`V=V(t,\{t^{pre}\},\{t^{post}\})`$ and $`f=f(t,\{t^{pre}\})`$, where the time dependence of $`f`$ is a double exponential.

In this section, however, we will focus on the simple scenario of a single pair of pre and postsynaptic APs, where the dynamics of the NMDAR and the BPAP are described as single exponentials with time constants $`\tau _\text{N}`$ and $`\tau _\text{B}`$, respectively. We will leave the more realistic two-component-BPAP case to the next section. Unless otherwise mentioned, we replace the stochastic variable $`f`$ by its mean without change of notation. Thus, for a single presynaptic spike at time $`t^{pre}`$, we can write,

$$f(t,t^{pre})=\mu \mathrm{\Theta }(tt^{pre})e^{(tt^{pre})/\tau _\text{N}},$$
(4)
where $`\mathrm{\Theta }`$ is the Heaviside function and $`\mu `$ is the probability that a channel will open given a presynaptic AP. Let us assume that the postsynaptic depolarization is dominated by the BPAP. Therefore, $`V`$ in equation 3 can be written as:

$`V(t,t^{post})`$ $`=V_R+\text{BPAP(t)}`$ (5)
$`=V_R+V_B\mathrm{\Theta }(tt^{post})e^{(tt^{post})/\tau _\text{B}},`$ (6)
where $`V_R`$ is resting membrane potential and $`V_B`$ is the magnitude of the BPAP due to a postsynaptic AP occuring at time $`t^{post}`$. In general, $``$ is a non-linear and non-monotonic function of $`V`$. For tractability, we will make a major simplification by linearizing the $``$-function in the interval $`[V_R,V_B]`$,

$$(V)=a+bV.$$
(7)
The original $``$ function used and its linear fit are displayed in figure 10. In appendix B, we assess the validity of this approximation. By using equations 4 through 7, equation 3 becomes:

$`I(t,t^{pre},`$ $`t^{post})`$
$`=`$ $`\overline{g}\left(a+bV_R+bV_B\mathrm{\Theta }(tt^{post})e^{(tt_{post})/\tau _\text{B}}\right)`$
$`\times \mu \mathrm{\Theta }(tt^{pre})e^{(tt^{pre})/\tau _\text{N}}`$
$`=`$ $`\overline{g}\mu \left(a+bV_R\right)\mathrm{\Theta }(tt^{pre})e^{(tt^{pre})/\tau _\text{N}}`$
$`+\overline{g}\mu bV_B\mathrm{\Theta }(tt^{post})\mathrm{\Theta }(tt^{pre})`$
$`\times e^{(tt^{pre})/\tau _\text{N}}e^{(tt^{post})/\tau _\text{B}}`$
$`=`$ $`I^{pre}(t)+I^{+/}(t).`$ (8)
The current separates into two components, one that depends only on the timing of the presynaptic AP ($`I^{pre}`$) and a second, associative term, that depends on the relative timing between the pre- and postsynaptic APs ($`I^{+/}`$). The associative term will give rise to the differences between pre/post and post/pre conditions, which is essential to the STDP learning.

Let $`\mathrm{\Delta }t=t^{post}t^{pre}`$. It is straightforward to show that:

$`I^{+/}=\{\begin{array}{cc}I_{peak}^+\mathrm{\Theta }(tt^{post})e^{(tt^{post})/\tau _1},\hfill & \text{if }\mathrm{\Delta }t>0\hfill \\ & \\ I_{peak}^{}\mathrm{\Theta }(tt^{pre})e^{(tt^{pre})/\tau _1},\hfill & \text{otherwise},\hfill \end{array}`$ (12)
where $`\tau _1^1=\tau _\text{B}^1+\tau _\text{N}^1`$, and $`I_{peak}^+=I_{peak}^+(\mathrm{\Delta }t)`$ and $`I^{}=I^{}(\mathrm{\Delta }t)`$ denote the peak magnitudes of the pre/post component of the calcium current for $`\mathrm{\Delta }t>0`$ and $`\mathrm{\Delta }t<0`$ respectively:

$`I_{peak}^+(\mathrm{\Delta }t)`$ $`=\overline{g}\mu bV_Be^{\mathrm{\Delta }t/\tau _\text{N}}`$ (13a)
$`I_{peak}^{}(\mathrm{\Delta }t)`$ $`=\overline{g}\mu bV_Be^{\mathrm{\Delta }t/\tau _\text{B}}.`$ (13b)
The dependence of the peak current on $`\mathrm{\Delta }t`$ is quite simple (FIG. 1). For $`\mathrm{\Delta }t>0`$ it decays exponentially with the time constant of the NMDAR open probability $`\tau _\text{N}`$; whereas for $`\mathrm{\Delta }t<0`$ it decays with $`\tau _\text{B}`$, the width of the BPAP. It is straightforward to see that, if $`\mathrm{\Delta }t<0`$ yields calcium levels corresponding to LTD, and if to each calcium concentration is associated a degree of plasticity, there should be a $`\mathrm{\Delta }t>0`$ region that also corresponds to LTD. The difference is that, since $`\tau _\text{B}<<\tau _\text{N}`$, $`\theta _d`$ in the post/pre region is reached for values of $`\mathrm{\Delta }t`$ much closer to zero than it is for the pre/post region. Such difference will be even greater in the case where the BPAP’s dynamic is described by two time constants instead of one (see section III).

Using equations 8 and 12, we can rewrite equation 2 in terms of each of the components of the calcium current,

$$\begin{array}{c}[Ca](t)=[Ca]^{pre}(t)+\mathrm{\Theta }(\mathrm{\Delta }t)[Ca]^+(t)\hfill \\ \hfill +\mathrm{\Theta }(\mathrm{\Delta }t)[Ca]^{}(t)+[Ca](0)e^{t/\tau }.\end{array}$$
(14)
with

$`[Ca]^{pre}(t)=e^{t/\tau }{\displaystyle I^{pre}(t^{})e^{t^{}/\tau }𝑑t^{}}`$ (15a)
$`[Ca]^+(t)=e^{t/\tau }{\displaystyle I^+(t^{})e^{t^{}/\tau }𝑑t^{}}`$ (15b)
$`[Ca]^{}(t)=e^{t/\tau }{\displaystyle I^{}(t^{})e^{t^{}/\tau }𝑑t^{}}.`$ (15c)
Integrating explicitly, we have,

$`[Ca]^{pre}(t)=\overline{g}\mu (a+bV_R)\tau _2\mathrm{\Theta }(tt^{pre})`$ (15a)
$`\times \left(e^{(tt^{pre})/\tau _\text{N}}e^{(tt^{pre})/\tau }\right)`$
$`[Ca]^+(t)=I_{peak}^+\tau _3\mathrm{\Theta }(tt^{post})`$ (15b)
$`\times \left(e^{(tt^{post})/\tau _1}e^{(tt^{post})/\tau }\right)`$
$`[Ca]^{}(t)=I_{peak}^{}\tau _3\mathrm{\Theta }(tt^{pre})`$ (15c)
$`\times \left(e^{(tt^{pre})/\tau _1}e^{(tt^{pre})/\tau }\right),`$
where $`\tau _2^1=\tau ^1\tau _\text{N}^1`$ and $`\tau _3^1=\tau ^1\tau _\text{N}^1\tau _\text{B}^1`$.

Examples of calcium transients are displayed in FIG. 2, where $`[Ca]^{pre/post}=[Ca]^{pre}+[Ca]^+`$ and $`[Ca]^{post/pre}=[Ca]^{pre}+[Ca]^{}`$ . Naturally, the $`[Ca]^{pre}`$ term is identical for the pre-post (FIG. 2a, $`\mathrm{\Delta }t`$ = 10 ms) and the post-pre (FIG. 2b, $`\mathrm{\Delta }t`$ = -10 ms) conditions (dashed line). The difference between these conditions lies in the associative term (dash-dot line). $`[Ca]^+`$ raises at the point $`t=t^{post}`$ and $`[Ca]^{}`$ raises at $`t=t^{pre}`$, their relative amplitudes being determined by $`I_{peak}^{+/}`$ (FIG. 1).

## III A back-propagating action potential with two components

The model described in the previous section assumes a BPAP with a single exponential falling phase. However, there are experimental indications that an additional slow after-depolarizing component exists in some cell bodies  and dendrites . The dynamics of the BPAP defines to what extent the postsynaptic spike-timing interacts with the presynaptic one; in particular, a slower decaying process ensures that such interaction spans for a period of tens of milliseconds, as required by the STDP learning rule. In this section, we assume a BPAP composed by a sum of two exponentials, with, respectively, a slow $`\tau _\text{B}^s`$ and a fast $`\tau _\text{B}^f`$ time constants. Equation 6 should therefore be re-written as:

$`V(t,t^{post})=`$ $`V_R+V_B\mathrm{\Theta }(tt^{post})`$ (6’)
$`\times \left[v^fe^{(tt^{post})/\tau _\text{B}^f}+v^se^{(tt^{post})/\tau _\text{B}^s}\right]`$
where $`v^s`$, $`v^f`$, $`\tau _\text{B}^s`$ and $`\tau _\text{B}^f`$ are the relative amplitudes and time constants of the slow and the fast components of the BPAP respectively, $`(v^s,v^f)>0,v^s+v^f=1`$. Analogously to the previous section, the calcium currents from equation 12 can be separated into two components, $`I^{pre}`$ and $`I^{+/}`$. $`I^{pre}`$ is identical to what we have derived before, whereas $`I^{+/}`$ reads:

$`I^{+/}(t)=\{\begin{array}{cc}\mathrm{\Theta }(tt^{post})I_{peak}^+[v^fe^{(tt^{post})/\tau _1^f}\hfill & \\ +v^se^{(tt^{post})/\tau _1^s}],\text{if }\mathrm{\Delta }t>0\hfill & \\ & \\ \mathrm{\Theta }(tt^{pre})[I_{peak}^fv^fe^{(tt^{pre})/\tau _1^f}\hfill & \\ +I_{peak}^sv^se^{(tt^{pre})/\tau _1^s}],\text{otherwise }\hfill & \end{array}`$ (20)
where $`I_{peak}^+`$ is as defined in equation 13a, $`I_{peak}^f=\overline{g}\mu bV_Be^{\mathrm{\Delta }t/\tau _\text{B}^f}`$, $`I_{peak}^s=\overline{g}\mu bV_Be^{\mathrm{\Delta }t/\tau _\text{B}^s}.`$ The calcium current $`I^{pre/post}`$ peaks at $`t=t^{post}`$ for $`\mathrm{\Delta }t>0`$ with amplitude $`I_{peak}^+`$, and at $`t=t^{pre}`$ for $`\mathrm{\Delta }t<0`$ with amplitude $`I_{peak}^{}=v^sI_{peak}^s+v^fI_{peak}^f`$. Thus, $`I_{peak}^{}`$ decays more abruptly than its single-component-BPAP counterpart (FIG. 3a).

We can again integrate each component of the calcium currents separately. We find for this case:

$`[Ca]^{pre}(t)=\overline{g}\mu (a+bV_R)\tau _2\mathrm{\Theta }(tt^{pre})`$ (15a’)
$`\times \left(e^{(tt^{pre})/\tau _\text{N}}e^{(tt^{pre})/\tau }\right)`$
$`[Ca]^+(t)=I_{peak}^+\mathrm{\Theta }(tt^{post})`$ (15b’)
$`\times \{v^f\tau _3^f(e^{(tt^{post})/\tau _1^f}e^{(tt^{post})/\tau })`$
$`+v^s\tau _3^s(e^{(tt^{post})/\tau _1^s}e^{(tt^{post})/\tau })\}`$
$`[Ca]^{}(t)=\mathrm{\Theta }(tt^{pre})`$ (15c’)
$`\times \{I_{peak}^fv^f\tau _3^f(e^{(tt^{pre})/\tau _1^f}e^{(tt^{post})/\tau })`$
$`+I^sv^s\tau _3^s(e^{(tt^{pre})/\tau _1^s}e^{(tt^{pre})/\tau })\}`$
where $`(\tau _1^f)^1=(\tau _\text{B}^f)^1+\tau _\text{N}^1`$, $`(\tau _3^f)^1=\tau ^1\tau _\text{N}^1(\tau _\text{B}^f)^1`$, and $`\tau _1^s`$ and $`\tau _3^s`$ are defined similarly to $`\tau _\text{B}^s`$, but with $`\tau _\text{B}^s`$ replacing $`\tau _\text{B}^f`$.

The effects of the BPAP slow component on the calcium transients is shown in FIG. 3b and 4a. Since both $`[Ca]^+`$ and $`[Ca]^{}`$ have the same linear dependence on the magnitude $`v_s`$ of this slow component, it primarily contributes to the separation between the pre-only condition (diamond), and the associative conditions (circle and square).

A longer decaying tail for the postsynaptic variable provides a longer range of time-interval in which interaction with the presynaptic variables is possible. Thus, the dynamics of the presynaptic variable, such as the calcium time constant $`\tau `$, will also influence the relative magnitudes of the calcium transients for the pre-post and post-pre conditions. Indeed, the combination of a faster calcium kinetics with the different BPAP time courses enhances the difference between these magnitudes (FIG. 4b). Therefore, a two-component-BPAP model produces a sharper transition of the learning curve between the post-pre LTD and the pre-post LTP. Initial estimates of the calcium time constant in the spines were of the order of 70-125ms . However, these values are obtained through calcium imaging, whose results may depend on the kinetics of the calcium indicators used, limiting the accuracy of the estimates. Recent experimental techniques, in which the effect of the calcium buffers are minimized, have lead to time constants of as low as 20 ms .

## IV Rate-dependence of the calcium transients

In the previous sections we have analyzed the calcium transients evoked by one pair of pre and postsynaptic spikes, these have significant implications for STDP at low frequency. However, STDP at low frequency is not the sole method for inducing synaptic plasticity, and in some cases even STDP requires that the pairs of pre and postspikes are delivered above a certain frequency . Thus it is important to analyze the contribution of spike delivery frequency $`\nu `$, to the calcium concentration, in addition to the relative timing effect. We analyze the case in which both pre and postsynaptic neurons fire regularly at the same frequency $`\nu `$, the onset of the postsynaptic spike-train being shifted from the presynaptic spike train by $`\mathrm{\Delta }t<1/\nu `$. Notice, however, that in the most general case the number of presynaptic and postsynaptic spikes as well as their timing differences will vary.

Recall that the NMDA current in equation 3 is the product between the fraction of opened channels $`f(t)`$ and the linearized Mg-block function $`(V(t))=a+bV(t)`$, where $`V(t)=V_R+\text{BPAP}(t)`$ omits depolarization due to the excitatory postsynaptic potentials (EPSPs). For a sequence $`\{t_1^{pre},t_2^{pre},\mathrm{},t_N^{pre}\}`$ of presynaptic spike-times and $`\{t_1^{post},t_2^{post},\mathrm{},t_M^{post}\}`$ of postsynaptic spike-times, we can write:

$$f(t)=\underset{n=1}{\overset{N}{}}A_n\mathrm{exp}\left(\frac{tt_n^{pre}}{\tau _\text{N}}\right)\mathrm{\Theta }\left(tt_n^{pre}\right),$$
(21)
$`\text{BPAP}(t)=V_B`$ $`{\displaystyle \underset{m=1}{\overset{M1}{}}}[\mathrm{exp}({\displaystyle \frac{tt_m^{post}}{\tau _\text{B}}})\mathrm{\Theta }(tt_m^{post})`$
$`.\times \mathrm{\Theta }(t+t_{m+1}^{post})]`$ (22)
$`+\mathrm{exp}\left({\displaystyle \frac{tt_M^{post}}{\tau _\text{B}}}\right)\mathrm{\Theta }\left(tt_M^{post}\right),`$
where $`N`$ and $`M`$ are the total number of pre and post-spikes before time $`t`$, and $`A_n`$ is the increase of the fraction of opened channels upon each pre-spike $`n`$. Let $`f(t)`$ be at level $`B_n`$ immediately before the pre-spike $`n`$ occurs. If the NMDA opening is proportional to the amount of non-opened channels, then $`A_n=\mu (1B_n)`$, where $`\mu `$ is the fraction of previously closed channels that open due to the presynaptic spike, $`\mu <1`$. Writing $`B_n`$ explicitly for each $`n`$, it is easy to see that it satisfies the following expression:

$`B_n=`$ $`\mu {\displaystyle \underset{i=1}{\overset{n1}{}}}e^{i/\nu \tau _\text{N}}`$
$`=`$ $`\mu e^{1/\nu \tau _\text{N}}{\displaystyle \frac{1(1\mu )^ne^{n/\nu \tau _\text{N}}}{1(1\mu )e^{1/\nu \tau _\text{N}}}}`$ (23)
and

$$B_{\mathrm{}}=\underset{n\mathrm{}}{lim}B_n=\frac{\mu e^{1/\nu \tau _\text{N}}}{1(1\mu )e^{1/\nu \tau _\text{N}}}$$
(24)
With this, equation 21 can be written as:

$`f(t)=`$ $`{\displaystyle \underset{n=1}{\overset{N1}{}}}\left(B_n+A_n\right)\mathrm{exp}\left({\displaystyle \frac{tt_n^{pre}}{\tau _\text{N}}}\right)`$
$`\times \mathrm{\Theta }(tt_n^{pre})\mathrm{\Theta }(t+t_{n+1}^{pre})`$
$`+`$ $`\left(A_N+B_N\right)\mathrm{exp}\left({\displaystyle \frac{tt_N^{pre}}{\tau _\text{N}}}\right)\mathrm{\Theta }(tt_N^{pre})`$
$`=`$ $`{\displaystyle \underset{n}{\overset{N1}{}}}\left[B_n(1\mu )+\mu \right]\mathrm{exp}\left({\displaystyle \frac{tt_n^{pre}}{\tau _\text{N}}}\right)`$
$`\times \mathrm{\Theta }(tt_n^{pre})\mathrm{\Theta }(t+t_{n+1}^{pre})`$ (25)
$`+`$ $`\left[B_N(1\mu )+\mu \right]\mathrm{exp}\left({\displaystyle \frac{tt_N^{pre}}{\tau _\text{N}}}\right)\mathrm{\Theta }(tt_N^{pre}).`$
We can define, analogously to equation 8 and the expressions 15a through 15c from section II, the following variables: $`I^{pre}=\overline{g}(a+bV_R)f(t)`$, $`I^{+/}=\overline{g}b\text{BPAP}(t)f(t)`$ and:

$`[Ca]^{pre}(t)`$ $`=e^{t/\tau }{\displaystyle I^{pre}(t^{})e^{t^{}/\tau }𝑑t^{}},`$ (26a)
$`[Ca]^{+/}(t)`$ $`=e^{t/\tau }{\displaystyle I^{pre/post}(t^{})e^{t^{}/\tau }𝑑t^{}}.`$ (26b)
Using equation 25, equation 26a reads:

$$\begin{array}{c}[Ca]^{pre}=\overline{g}e^{t/\tau }\left(a+bV_R\right)\tau _2\hfill \\ \hfill \times \{(e^{1/\nu \tau _2}1)\underset{n=1}{\overset{N1}{}}[B_n(1\mu )+\mu ]e^{t_n^{pre}/\tau }\\ \hfill +\left(e^{t/\tau _2}e^{t_N^{pre}/\tau _2}\right)\left[B_N(1\mu )+\mu \right]\\ \hfill .\times e^{t_N^{pre}/\tau _\text{N}}\}.\end{array}$$
(27)
We use the approximation $`B_n=B_{\mathrm{}}`$, and write the time indexes with respect to the timing of the first spike:

$`t_n^{pre}=t_1^{pre}+{\displaystyle \frac{n1}{\nu }},`$ (28a)
$`t_N^{pre}=t_1^{pre}+{\displaystyle \frac{\nu t1}{\nu }}\text{ and}`$ (28b)
$`N=(tt_1^{pre})\nu +1,`$ (28c)
where $``$, is the floor operation. Evaluating the sum, we have:

$`[Ca]^{pre}=`$ $`\overline{g}e^{t/\tau }(a+bV_R)\tau _2\left[B_{\mathrm{}}(1\mu )+\mu \right]`$
$`\times \{(e^{1/(\nu \tau _2)}1)e^{t_1^{pre}/\tau }\left[{\displaystyle \frac{1\mathrm{exp}\left(\frac{N1}{\nu \tau }\right)}{1\mathrm{exp}\left(\frac{1}{\nu \tau }\right)}}\right]`$
$`.+\mathrm{exp}({\displaystyle \frac{t_N^{pre}}{\tau _\text{N}}}+{\displaystyle \frac{t}{\tau _2}})\mathrm{exp}\left({\displaystyle \frac{t_N^{pre}}{\tau }}\right)\},`$ (29)
where $`N`$ and $`t_N^{pre}`$ are defined as in equations 28b and 28c.

To calculate the contribution to the calcium transients due to the interaction of the pre and postsynaptic spikes, we substitute the expression 22 and 25 into equation 26b. Let $`t_m^{post}t_n^{pre}+\mathrm{\Delta }t`$ and $`t_M^{post}t_N^{pre}+\mathrm{\Delta }t`$. The indexes $`m`$ and $`M`$ will always refer to the post-spikes while $`n`$ and $`N`$ will always refer to the pre-spikes; for simplicity, we drop the superscripts so that $`t_n^{pre}=t_n`$ and $`t_m^{post}=t_m`$. Because the depolarization due to each post-spike does not build up, the product between the sums given by expressions 25 and 22 will yield only four terms, due to the interaction of the pairs $`\{(n,m1);(n,m);(N,M1);(N,M)\}`$. These terms are:

$`{\displaystyle _0^t}f(t^{})\text{BPAP}(t^{})dt^{}=V_B[B_{\mathrm{}}(1\mu )+\mu ]\{.`$
$`{\displaystyle \underset{n=2}{\overset{N1}{}}}{\displaystyle _{t_n}^{t_m}}\mathrm{exp}\left({\displaystyle \frac{t^{}t_n}{\tau _\text{N}}}{\displaystyle \frac{t^{}t_{m1}}{\tau _\text{B}}}+{\displaystyle \frac{t^{}}{\tau }}\right)𝑑t^{}`$ (30a)
$`+{\displaystyle \underset{n=1}{\overset{N1}{}}}{\displaystyle _{t_m}^{t_{n+1}}}\mathrm{exp}\left({\displaystyle \frac{t^{}t_n}{\tau _\text{N}}}{\displaystyle \frac{t^{}t_m}{\tau _\text{B}}}+{\displaystyle \frac{t^{}}{\tau }}\right)𝑑t^{}`$ (30b)
$`+{\displaystyle _{t_N}^t_{}}\mathrm{exp}\left({\displaystyle \frac{t^{}t_N}{\tau _\text{N}}}{\displaystyle \frac{tt_{M1}}{\tau _\text{B}}}+{\displaystyle \frac{t^{}}{\tau }}\right)𝑑t^{}`$ (30c)
$`+\delta _{M,N}{\displaystyle _{t_M}^t}\mathrm{exp}({\displaystyle \frac{t^{}t_N}{\tau _\text{N}}}{\displaystyle \frac{t^{}t_M}{\tau _\text{B}}}+{\displaystyle \frac{t^{}}{\tau }})dt^{}.\},`$ (30d)
where $`t_{}=\mathrm{min}(t,t_M)`$ and we have used that $`B_n=B_N=B_{\mathrm{}}`$. Recalling expressions 28a through 28c, one can easily perform the integrations 30a through 30d. The remaining sums are finite power series, therefore they converge. Thus, the calcium concentration due to the interaction between pre and postsynaptic spikes will be:

$$\begin{array}{c}[Ca]^{+/}=\overline{g}e^{t/\tau }bV_B\left[B_{\mathrm{}}(1\mu )+\mu \right]\hfill \\ \hfill \left\{\text{30a}\text{ + }\text{30b}\text{ + }\text{30c}\text{ + }\text{30d}\right\}\end{array}$$
(31)
with each of the terms being:

$$\begin{array}{c}\text{30a}=\tau _3\left[\mathrm{exp}\left(\frac{\mathrm{\Delta }t}{\tau _2}\frac{1}{\nu \tau _\text{B}}\right)\mathrm{exp}\left(\frac{\mathrm{\Delta }t1/\nu }{\tau _\text{B}}\right)\right]\hfill \\ \hfill \times e^{t_1^{pre}/\tau }\left[\frac{1\mathrm{exp}\left(\frac{N1}{\nu \tau }\right)}{1\mathrm{exp}\left(\frac{1}{\nu \tau }\right)}\right]\end{array}$$
(32)
$$\begin{array}{c}\text{30b}=\tau _3\left[\mathrm{exp}\left(\frac{\mathrm{\Delta }t}{\tau _\text{B}}+\frac{1}{\nu \tau _3}\right)\mathrm{exp}\left(\frac{\mathrm{\Delta }t}{\tau _2}\right)\right]\hfill \\ \hfill \times e^{t_1^{pre}/\tau }\left[\frac{1\mathrm{exp}\left(\frac{N1}{\nu \tau }\right)}{1\mathrm{exp}\left(\frac{1}{\nu \tau }\right)}\right]\end{array}$$
(33)
$`\text{30c}=\{\begin{array}{cc}\tau _3\left[\mathrm{exp}\left(\frac{\mathrm{\Delta }t}{\tau _2}\frac{1}{\nu \tau _\text{B}}\right)\mathrm{exp}\left(\frac{\mathrm{\Delta }t1/\nu }{\tau _\text{B}}\right)\right]\hfill & \\ \times \mathrm{exp}\left(\frac{t_N^{pre}}{\tau }\right),\text{ if }\left(t_N^{pre}+\mathrm{\Delta }t\right)<t\hfill & \\ & \\ \tau _3\left[\mathrm{exp}\left(\frac{t_N^{pre}}{\tau _1}+\frac{t}{\tau _3}\right)\mathrm{exp}\left(\frac{t_N^{pre}}{\tau }\right)\right]\hfill & \\ \times \mathrm{exp}\left(\frac{\mathrm{\Delta }t1/\nu }{\tau _\text{B}}\right),\text{if }\left(t_N^{pre}+\mathrm{\Delta }t\right)>t\text{.}\hfill & \end{array}`$
$`\text{30d}=\{\begin{array}{cc}\tau _3\left[\mathrm{exp}\left(\frac{t_M^{post}}{\tau _1}+\frac{t}{\tau _3}\right)\mathrm{exp}\left(\frac{t_M^{post}}{\tau }\right)\right]\hfill & \\ \times \mathrm{exp}(\frac{\mathrm{\Delta }t}{\tau _\text{N}}),\text{ if }\left(t_N^{pre}+\mathrm{\Delta }t\right)<t\hfill & \\ & \\ 0,\text{if }\left(t_N^{pre}+\mathrm{\Delta }t\right)>t\text{.}\hfill & \end{array}`$
Plots of these expressions are shown in FIG. 5. As expected, the overall calcium is frequency and timing-dependent. There is temporal integration of calcium levels, so that higher frequency stimulations build up the intracellular calcium concentration (FIG. 5a). In addition, the $`\mathrm{\Delta }t`$-dependence shown in previous sections (FIG. 5b) for single pairs of pre and post spikes is retained. The amount of temporal integration will naturally depend on the specific values of $`\tau `$, $`\tau _\text{N}`$ and $`\tau _\text{B}`$, and so will the difference between the pre-post and the post-pre conditions, as analyzed before. For example, slower dynamics would result in more time integration at moderate frequencies. However, given these constants, it is possible to set a threshold of calcium concentration between LTD and LTP (or between no-plasticity and LTD) that would correspond to low- and high-frequency STDPs, respectively.

We now extract an expression for the dependence of calcium transients on $`\mathrm{\Delta }t`$ and on $`\nu `$ in the limit $`t\mathrm{}`$. Let $`t=t_N^{pre}+\delta `$, where the new variable $`\delta `$ tracks the time since the last presynaptic spike, and define $`\gamma _1=\overline{g}(a+bV_R)\tau _2[B_{\mathrm{}}(1\mu )+\mu ]`$ and $`\gamma _2=\overline{g}bV_B\tau _3[B_{\mathrm{}}(1\mu )+\mu ]`$. For $`t\mathrm{}`$ and $`\mathrm{\Delta }t<\delta `$, we have:

$`[Ca](\delta )=`$ $`\gamma _1e^{\delta /\tau }\left\{e^{\delta /\tau _2}+\left[{\displaystyle \frac{e^{1/\nu \tau _2}e^{1/\nu \tau }}{e^{1/\nu \tau }1}}\right]\right\}`$
$`+\gamma _2e^{\delta /\tau }e^{\mathrm{\Delta }t/\tau _\text{B}}\{.\left[{\displaystyle \frac{e^{1/\nu \tau _3}e^{\mathrm{\Delta }t/\tau _3}}{e^{1/\nu \tau }1}}\right]`$
$`+\mathrm{exp}\left({\displaystyle \frac{1}{\nu \tau _\text{B}}}+{\displaystyle \frac{1}{\nu \tau }}\right)\left[{\displaystyle \frac{e^{\mathrm{\Delta }t/\tau _3}1}{e^{1/\nu \tau }1}}\right]`$
$`+(e^{\delta /\tau _3}e^{\mathrm{\Delta }t/\tau _3}).\}`$ (34)
And for $`\mathrm{\Delta }t>\delta `$:

$`[Ca](\delta )=`$ $`\gamma _1e^{\delta /\tau }\left\{e^{\delta /\tau _2}+\left[{\displaystyle \frac{e^{1/\nu \tau _2}e^{1/\nu \tau }}{e^{1/\nu \tau }1}}\right]\right\}`$
$`+\gamma _2e^{\delta /\tau }e^{\mathrm{\Delta }t/\tau _\text{B}}\{.[e^{\mathrm{\Delta }t/\tau _3}1]\left[{\displaystyle \frac{e^{1/\nu \tau _\text{B}}}{e^{1/\nu \tau }1}}\right]`$
$`+\left[{\displaystyle \frac{e^{1/\nu \tau _3}e^{\mathrm{\Delta }t/\tau _3}}{e^{1/\nu \tau }1}}\right]`$
$`+(e^{\delta /\tau _3}1)e^{1/\nu \tau _\text{B}}.\}`$ (35)
These expressions show how the peaks of the calcium transients depend on the timing and on the frequency of the pre and postsynaptic spikes. As the frequency increases (FIG. 6a) the $`\mathrm{\Delta }t`$-dependent curves move up due to temporal integration, and the difference between pre-post and post-pre (FIG. 6b) becomes smaller.

## V Stochastic properties of calcium influx

In previous sections we have calculated the mean calcium transients under various conditions and assumptions. However, synaptic transmission is a random process and the trial-by-trial calcium transients deviate from the mean. One source of variability is the stochasticity of the presynaptic neurotransmitter release. A failure of release at a given synapse would result in eliminating the postsynaptic calcium transient at that synapse. Although this source of variability can be significant, its consequences are simple to predict, because the absence of release will result in no synaptic plasticity in the specific spine. Another source of variablity is the stochasticity of the postsynaptic opening of NMDA channels. Because the number of NMDAR at each synapse can be small ($`10`$) , these fluctuations could be considerable.

In this section, we calculate the variability due to $`Z`$ NMDAR in the postsynaptic terminal. We assume the following simple Markov model: the state vector $`𝐏_n(t)`$ denotes the probability of being in each one of the $`n`$ possible states and the transition matrix $`R_n`$ represents transition probabilities from state to state. For a Markov processes, the evolution equation has the form:

$$\frac{d𝐏_n}{dt}=R_n𝐏_n,$$
(36)
with the solution:

$$𝐏_n(t)=𝐏_n(0)\mathrm{exp}\left(_0^tR_n𝑑t^{}\right),$$
(37)
where $`𝐏_n(0)`$ is the initial condition.

In the simplest model the process has only two states, unbound $`u`$ and open-bound $`o`$:

$$u\begin{array}{c}k_1G\\ \\ \\ k_1\end{array}o$$
(38)
where $`G`$ is the glutamate concentration, $`k_1`$ and $`k_1`$ are the forward and backward time constants respectively.

This model is clearly simplified; typically, a five-state model is used to account for NMDAR kinetics. However, it is sufficient for describing the single exponential decay kinetics we have assumed for the NMDAR since, during decay, $`G=0`$ (see equation 41, below). The transition matrix for this model is:

$$R_2=\left(\begin{array}{cc}k_1G& k_1\\ k_1G& k_1.\end{array}\right)$$
(39)
where $`G`$ is the neurotransmitter concentration. We assume that $`G`$ has a constant, non-zero value for only a limited amount of time. We are primarily interested in the falling phase ($`G=0`$) of the NMDAR current.

We denote by $`p^u`$ and $`p^o`$ the probability of being in the unbound and open states respectively. Since $`p^u+p^o=1`$ the differential equation for $`p^o`$ reads:

$$\dot{p^o}=k_1Gp^o(k_1G+k_1).$$
(40)
For a constant value of $`G`$, the solution of this equation is:

$$p^o(t)=\frac{k_1G}{k_1G+k_1}+Ce^{(k_1G+k_1)t}.$$
(41)
where the constant $`C`$ is determined by the initial conditions.

During the falling phase ($`G=0`$) the time constant has the value $`\tau _\text{N}=1/k_1`$. The solution takes the form $`p^o(t)=\mu e^{t/\tau _\text{N}}`$, where $`\mu `$ can now be defined more precisely as the open probability at the start of the falling phase. The total calcium current through $`Z`$ identical NMDAR is $`I(t)=_{i=1}^Zs_ig(V(t))`$, where $`g`$ is the single channel NMDAR conductance, and $`s_i=1`$ when receptor in open state 0 when in closed and unbound states. To simplify the notation we can rewrite equation 7 as $`=V_0+V_1\mathrm{\Theta }(tt^{post})e^{(tt^{post})/\tau _\text{B}}`$, where $`V_0=a+bV_R`$ and $`V_1=bV_B`$.

Calcium concentration is assumed to be governed by equation 1, therefore the average calcium concentration will be:

$$[Ca](t)=ge^{t/\tau }_{\mathrm{}}^t𝑑t^{}e^{t^{}/\tau }\underset{i=1}{\overset{Z}{}}s_i(t^{})(t^{}),$$
(42)
where $``$ denotes averaging over the probability distribution of $`s_i(t)`$.

However, $`s_i(t)=Zs_i(t)=Zp^o(t)`$, thus the calcium concentration given $`Z`$ NMDAR is:

$$[Ca](t)_Z=Zge^{t/\tau }_{\mathrm{}}^t𝑑t^{}e^{t^{}/\tau }p^o(t^{})(t^{}).$$
(43)
Expression 43 is similar to equation 2 since $`I(t)=Zg(t)p^o(t)`$. We have calculated this integral in section II for the different cases (pre, pre-post, post-pre). These assumptions are equivalent to the simple Markov model assumed here. This result is shown in equations 42 and 43. We note that $`[Ca](t)_Z=Z[Ca](t)_1`$.

The variance for $`Z`$ NMDAR, has the form $`\sigma _Z^2(t)=[Ca]^2(t)_Z[Ca](t)_Z^2`$. The second moment of the calcium concentration for $`Z`$ NMDAR is:

$`[Ca]^2(t)_Z=`$ $`g^2e^{2t/\tau }{\displaystyle _{\mathrm{}}^t}𝑑t^{}𝑑t^{\prime \prime }e^{t^{}/\tau }e^{t^{\prime \prime }/\tau }(t^{})(t^{\prime \prime })`$
$`\times {\displaystyle \underset{i=1}{\overset{Z}{}}}s_i(t^{}){\displaystyle \underset{j=1}{\overset{Z}{}}}s_j(t^{\prime \prime })`$
$`=`$ $`g^2e^{2t/\tau }{\displaystyle _{\mathrm{}}^t}𝑑t^{}𝑑t^{\prime \prime }e^{t^{}/\tau }e^{t^{\prime \prime }/\tau }(t^{})(t^{\prime \prime })`$
$`\times {\displaystyle \underset{i,j=1}{\overset{Z}{}}}p_{ij}^{oo}(t^{},t^{\prime \prime }),`$ (44)
where $`p_{ij}^{oo}(t,t^{})`$ is the joint probability of having channel $`i`$ open at time $`t`$ and channel $`j`$ open at time $`t^{}`$. We assume independent identical channels, therefore for $`ij`$ we have $`p_{ij}^{oo}(t,t^{})=p^o(t)p^o(t^{})`$, and when $`i=j`$, $`p_{ii}^{oo}(t,t^{})`$ is defined as $`p^{oo}(t,t^{})`$. We therefore have that:

$`[Ca]^2(t)_Z=`$ $`g^2e^{2t/\tau }{\displaystyle _{\mathrm{}}^t}𝑑t^{}𝑑t^{\prime \prime }e^{t^{}/\tau }e^{t^{\prime \prime }/\tau }(t^{})(t^{\prime \prime })`$
$`\times \left[Z(Z1)p^o(t^{})p^o(t^{\prime \prime })+Zp^{oo}(t^{},t^{\prime \prime })\right]`$
$`=`$ $`g^2Z(Z1)[Ca](t)^2+g^2Ze^{2t/\tau }`$
$`\times {\displaystyle _{\mathrm{}}^t}dt^{}dt^{\prime \prime }e^{t^{}/\tau }e^{t^{\prime \prime }/\tau }(t^{})(t^{\prime \prime })p^{oo}(t^{},t^{\prime \prime })`$ (45)
Thus the variance of calcium concentration for $`Z`$ NMDAR is:

$`\sigma _Z^2(t)=`$ $`g^2Ze^{2t/\tau }{\displaystyle _{\mathrm{}}^t}𝑑t^{}𝑑t^{\prime \prime }e^{t^{}/\tau }e^{t^{\prime \prime }/\tau }(t^{})(t^{\prime \prime })`$
$`\times \left[p^{oo}(t^{},t^{\prime \prime })p^o(t^{})p^o(t^{\prime \prime })\right]`$
$`=`$ $`g^2Ze^{2t/\tau }{\displaystyle _{\mathrm{}}^t}𝑑t^{}𝑑t^{\prime \prime }(t^{})(t^{\prime \prime })e^{t^{}/\tau }e^{t^{\prime \prime }/\tau }p^{oo}(t^{},t^{\prime \prime })`$
$`Z[Ca](t)^2.`$ (46)
We can use the Bayes rule to rewrite $`p^{oo}(t,t^{})=p^{oo}(t|t^{})p^o(t^{})`$. Thus there are two types of variables which we need to calculate from the Markov process, $`p^o(t)`$ and $`p^{oo}(t|t^{})`$. The solution for $`p^o(t)`$ is given by equation 41. We define the time $`t=0`$ as the beginning of the falling phase ($`G=0`$), or $`t^{pre}=0`$. Further, we assume an instantaneous rise time. So for $`t<0`$, $`p^o(t)=0`$ and at $`t=0`$ $`p^o(0)=\mu `$. Thus for $`t>0`$ $`p^o(t)=\mu e^{t/\tau _\text{N}}`$.

Similarly, $`p^{oo}(t|t^{})`$ is the probability of being in open state at time $`t`$ given that the channel was in open state at time $`t^{}`$. We calculate this for the case $`t,t^{}0`$. If $`t^{}>t`$ this can be solved in a similar manner to $`p^o`$, since in this regime it is a stationary Markov process, it must be a function of $`(tt^{})`$, and at $`t=t`$, $`p^{oo}(t=t^{}|t^{})=1`$. During the falling phase there is one directional movement from $`o`$ to $`u`$, since $`G=0`$. Therefore, if we know that at time $`t^{}`$ the channel is in an open state it must have been in an open state, at any time $`t<t^{}`$. Therefore:

$$p^{oo}(t|t^{})=\{\begin{array}{cc}e^{(tt^{})/\tau _\text{N}}:\hfill & t>t^{}\hfill \\ 1:\hfill & tt^{}\hfill \end{array}$$
(47)
Note that this holds for $`(t,t^{})>0`$, and that by definition at times smaller than zero the receptor is unbound. Hence:

$$p^{oo}(t,t^{})=\{\begin{array}{cc}\hfill \mu e^{t/\tau _\text{N}}:& t>t^{}\hfill \\ \hfill \mu e^{t^{}/\tau _\text{N}}:& t<t^{}\hfill \end{array}$$
(48)
This expression can now be substituted into equation 46, and by using the linear approximation of $``$ we obtain an expression for $`\sigma _Z^2(t)`$ that can be computed analytically. To aid in calculation we denote:

$$\begin{array}{c}\sigma _Z^2(t)=Z\mu g^2e^{2t/\tau }(V_0^2\alpha _1+2V_0V_1\alpha _2+V_1^2\alpha _3)\hfill \\ \hfill Z[Ca](t)^2\end{array}$$
(49)
where:

$$\alpha _1=_0^t𝑑t^{}_0^t𝑑t^{\prime \prime }p^{oo}(t^{},t^{\prime \prime })e^{t^{}/\tau }e^{t^{\prime \prime }/\tau }\mathrm{\Theta }(t)$$
(50a)
$$\alpha _2=_0^t𝑑t^{}_0^t𝑑t^{\prime \prime }e^{(t^{}t^{post})/\tau _\text{B}}p^{oo}(t^{},t^{\prime \prime })$$
(50b)
$`\alpha _3=`$ $`{\displaystyle _0^t}𝑑t^{}{\displaystyle _0^t}𝑑t^{\prime \prime }e^{(t^{}t^{post})/\tau _\text{N}}e^{(t^{\prime \prime }t^{post})/\tau _\text{B}}p^{oo}(t^{},t^{\prime \prime })`$
$`\times e^{t^{}/\tau }e^{t^{\prime \prime }/\tau }\mathrm{\Theta }(t^{}t^{post})\mathrm{\Theta }(t^{\prime \prime }t^{post})`$ (50c)
In appendix A we analytically calculate the form of the terms $`\alpha _1`$, $`\alpha _2`$ and $`\alpha _3`$. Although the analytical form is complex, the consequences are simple and significant. As expected, $`CV_Z=CV_1/\sqrt{Z}`$. Thus the variability is significant only at relatively low $`Z`$. In FIG. 7 we show $`[Ca](t)_Z`$ and $`\sigma _Z`$, for $`\mathrm{\Delta }t=`$ -10, 10 and 60 ms. In all cases both the mean and $`\sigma _Z`$ change over time. The peak of $`[Ca](t)_Z`$ for $`\mathrm{\Delta }t=10`$ ms and $`\mathrm{\Delta }t=60`$ ms are similar. However, the magnitude of the variability, $`\sigma _Z`$, is significantly larger at $`\mathrm{\Delta }t=60`$ ms than at $`\mathrm{\Delta }t=10`$ ms. To quantify the dependence of the relative variability as a function of $`\mathrm{\Delta }t`$, we measure the coefficient of variation, $`CV_Z=\sigma _Z/[Ca](t)_Z`$, where both $`[Ca](t)_Z`$ and $`\sigma _Z(t)`$ are measured at the time $`t`$ where $`[Ca](t)_Z`$ is maximal.

When $`\mathrm{\Delta }t<0`$ the $`CV`$ is low and decreases as $`\mathrm{\Delta }t0`$ from below. For $`\mathrm{\Delta }t>0`$ the $`CV`$ increases as $`\mathrm{\Delta }t`$ increases (FIG 8). It is interesting to compare the $`CV`$ for values of $`\mathrm{\Delta }t<0`$ and $`\mathrm{\Delta }t>0`$ which have similar peak calcium levels. For example for $`\mathrm{\Delta }t=10`$ ms, $`CV_{10}=0.34`$ and for $`\mathrm{\Delta }t=60`$ ms, which has a similar peak calcium level, $`CV_{10}=0.51`$. Therefore $`CV_{10}(60)/CV_{10}(10)1.5`$. Since $`CV_Z=CA_1/\sqrt{Z}`$, for this set of parameters, but for every $`Z`$ $`CV_Z(60)/CV_Z(10)1.5`$.

The variability of calcium transients also depends on $`\mu `$ (FIG. 9), decreasing as $`\mu `$ increases, for all $`\mathrm{\Delta }t`$. The variable $`\mu `$ is the probability of glutamate binding to postsynaptic receptors given a presynaptic spike, and can be taken to be the presynaptic probability of release. Thus, as the presynaptic probability of release increases, $`CV`$ decreases.

## VI Discussion

Calcium transients due to NMDA currents are believed to play a major role in many forms of synaptic plasticity. We have recently shown that a model based on these transients can account for various different plasticity-induction protocols . In this paper, we compute analytically the calcium dynamics evoked by pairs of pre- and postsynaptic spikes under different conditions, so that the variables that control calcium transients and impact synaptic plasticity can be investigated.

We showed that the peak of the calcium transients depends on the $`\mathrm{\Delta }t`$: for $`\mathrm{\Delta }t>0`$, the peak concentration decays with the time constant $`\tau _\text{N}`$ of the NMDAR, and for $`\mathrm{\Delta }t<0`$, it decays with the time constant $`\tau _\text{B}`$ of the BPAP. Therefore, if spike time dependent plasticity, indeed depends on the calcium transients, $`\tau _\text{N}`$ and $`\tau _\text{B}`$ will determine the width of the pre-post LTP and the post-pre LTD windows, respectively. These results were confirmed by simulations, in which the approximations we used in the calculations were relaxed. The $`\mathrm{\Delta }t`$ dependence of the peak calcium transient is similar in the more complex and realistic case where the BPAP dynamics is composed of a fast and a slow component: the difference lies in the $`\mathrm{\Delta }t<0`$ case, which now depends on the combination of the two time constants. The two-component BPAP allows a sharper transition between the post-before-pre and the pre-before-post STDP window, which is further enhanced by a fast calcium dynamics. We also showed that, besides the timing of the pre and post-spikes, the peak values of the calcium transients also depend on the frequency $`\nu `$ of the pre and postsynaptic conditioning, increasing for greater values of $`\nu `$. However, at higher frequencies, the dependence on $`\mathrm{\Delta }t`$ decreases. The amount of calcium build-up at a given frequency depends on the system parameters; slower time constants results in more temporal integration. Temporal integration of calcium transients, as described here, can explain why, in some cases, the induction of STDP is frequency-dependent.

One of the predictions of the unified calcium model is that there exists a pre-before-post LTD for values of $`\mathrm{\Delta }t`$ greater than the ones that elicit LTP . This has been shown experimentally by some investigators ; others have only placed a small number of data points in this region . In all published STDP experiments, there is a large variability in the magnitude and the sign of synaptic plasticity across the different values of $`\mathrm{\Delta }t`$. It would therefore be difficult to assess the existence of this type of LTD without a large amount of data. The variability encountered in the experiments indicates that it should be important to examine the stochastic properties of the calcium transients, in addition to their mean.

We calculated the variance of the calcium transients and showed that for a small number $`Z`$ of NMDAR, it can be quite significant. There are indications that the number of NMDAR in the dendrite is indeed small ($`10`$) . Further, we showed that the $`CV`$ increases monotonically with $`\mathrm{\Delta }t`$, for $`\mathrm{\Delta }t>0`$. These fluctuations can have significant implications on the outcome of downstream processes dependent on calcium, such as phosphorylation or other second-messenger cascades upon which plasticity relies.

Previously, several authors have simulated and analyzed dynamical models of synaptic plasticity that can lead to STDP . These models typically bypass detailed descriptions of the physiology and biochemistry underlying synaptic plasticity. Recently, biophysical models have been receiving increasing focus . These more mechanistic approaches have allowed one to show that various different induction protocols can be accounted for by a single set of assumptions regarding the dependence of plasticity on calcium concentrations . In these models, plasticity is a highly non-linear function of calcium. Therefore, fluctuations of calcium transients will not only add spread to the resulting LTD and LTP, but will also shift their mean (the mean of the solution is not equivalent to the solution given the mean calcium level). For example, in the unified calcium model , the sign and magnitude of synaptic plasticity is determined by a saturated U-shaped function of calcium concentration. When calcium falls below a threshold $`\theta _d`$, no plasticity occurs, when it falls between $`\theta _d`$ and $`\theta _p`$, $`\theta _p>\theta _d`$, LTD is induced, and for calcium above $`\theta _p`$, LTP is induced. Thus, in the LTD region, if the spread of the calcium is of the order of $`\theta _p\theta _d`$, occasional LTP can still happen. We showed an explicit example for the cases of $`\mathrm{\Delta }t=10`$ ms and $`\mathrm{\Delta }t=60`$ ms, where the mean peak calcium concentrations are similar, but the variability for the latter is larger. If this mean calcium is in the LTD range, it is likely that the observed phenomenon is that the LTD is stronger and more robust for $`\mathrm{\Delta }t=10`$ ms than for $`\mathrm{\Delta }t=60`$ ms.

This work is based on several simplifying assumptions. This allows us a qualitative understanding of the aspects of calcium dynamics that are crucial for the induction of synaptic plasticity. Incorporation of more realistic features, such as those assumed in the unified calcium model, would probably yield quantitatively more precise results. For example, a double-component NMDAR decay would allow for more calcium build-up without sacrificing the temporal asymmetry between pre-post and post-pre conditions. However, we believe that, despite the constraints of analytical tractability, this study can serve as a basis for future models that aim to formalize the biophysical basis of synaptic plasticity.

## VII Acknowledgments

HZS would like to thank the members of the IBNS at Brown University as well as M. Tsodyks, D. Golomb and Y. Kubota for useful discussions and comments, and the ITP at Santa Barbara where some of this work was carried out. This work was partly supported by the Burroughs-Wellcome Fund and the Charles A. Dana Foundation.

## Appendix A Variance calculation

Using the definition in equation 50a through V we obtain:

$`\alpha _1=`$ $`[{\displaystyle _0^t}dt^{}{\displaystyle _t^{}^t}dt^{\prime \prime }e^{t^{\prime \prime }/\tau _\text{N}}e^{t^{}/\tau }e^{t^{\prime \prime }/\tau }`$
$`+{\displaystyle _0^t}dt^{}{\displaystyle _0^t^{}}dt^{\prime \prime }e^{t^{}/\tau _\text{N}}e^{t^{}/\tau }e^{t^{\prime \prime }/\tau }]\mathrm{\Theta }(t)`$
$`=`$ $`[{\displaystyle _0^t}dt^{}e^{t^{}/\tau }\tau _2(e^{t/\tau _2}e^{t^{}/\tau _2})`$
$`+{\displaystyle _0^t}dt^{}e^{t^{}/\tau _2}\tau (e^{t^{}/\tau }1)dt^{}]\mathrm{\Theta }(t)`$
$`=`$ $`[e^{t/\tau _2}\tau _2\tau ((e^{t/\tau }1)\tau _2{\displaystyle _0^t}dt^{}e^{t^{}/\tau _4}`$
$`+\tau {\displaystyle _0^t}dt^{}e^{t^{}/\tau _4}\tau {\displaystyle _0^t}dt^{}e^{t^{}/\tau _2}]\mathrm{\Theta }(t)`$
$`=`$ $`[e^{t/\tau _4}(\tau _2\tau \tau _2\tau _4+\tau \tau _4)`$
$`2e^{t/\tau _2}\tau _2\tau +(\tau _2\tau _4\tau \tau _4+\tau _2\tau )]\mathrm{\Theta }(t),`$ (51)
where $`\tau _4^1=2\tau ^1\tau _\text{N}^1`$.

$`\alpha _2=`$ $`e^{t^{post}/\tau _\text{B}}[\tau _5(\tau \tau _2)(e^{t/\tau _5}e^{\stackrel{~}{t}/\tau _5})`$
$`\tau \tau _3(e^{t/\tau _3}e^{\stackrel{~}{t}/\tau _3})+\tau _2\tau _6e^{t/\tau _2}(e^{t/\tau _6}e^{\stackrel{~}{t}/\tau _6})]`$
$`\times \mathrm{\Theta }(t\stackrel{~}{t})`$ (52)
where $`\tau _3`$ is as defined before, $`\tau _5^1=2\tau ^1\tau _\text{N}^1\tau _\text{B}^1`$, $`\tau _6^1=\tau ^1\tau _\text{B}^1`$ and

$`\stackrel{~}{t}=\{\begin{array}{cc}\hfill t^{post}:& t^{post}>0\hfill \\ \hfill 0:& t^{post}0\hfill \end{array}`$ (55)
Finally,

$`\alpha _3=2e^{2t^{post}/\tau _\text{B}}\tau _6\tau _3`$ $`\left(e^{t/\tau _6}e^{t/\tau _3}e^{\stackrel{~}{t}/\tau _6}+e^{\stackrel{~}{t}/\tau _6}\right)`$
$`\times \mathrm{\Theta }(t\stackrel{~}{t})`$ (56)
## Appendix B Numerical estimation of the error resulting from the linear approximation of $``$

In the simulations, we used the following form for the voltage-dependence of the NMDAR.:

$$(V)=\frac{VV_{rev}}{1+\frac{e^{0.062V}}{3.57}}$$
(57)
where $`V_{rev}=130`$ mV, fit performed in the interval $`[70,10]`$ mV.

Such functional dependence captures the qualitative non-linear dependence of the NMDAR on the voltage, although the reversal potential for calcium channels might be inaccurate at the values close to it. In our analysis, we approximated $``$ with a linear dependence on voltage. The fit used is shown in FIG. 10a. We see that over the relevant range (\[-70 0\]) this is a relatively good approximation.

In panels b-d of FIG. 10 we compare the numerically extracted calcium transients (dashed line), using the non-linear $``$ and the analytical results obtained for the linear approximation (solid line). At rest (FIG. 10b), the linear approximation significantly underestimates the magnitudes of the transients, because this is where the linear fit departs the most from the true curve. However, in the working regime (\[-10, 10\] mV), the approximation can be considered adequate.

Let equation 57 be rewritten as $`(v)=a+bV`$ as in equation 7. The parameters of linearization are $`a=0.1031/\overline{g}`$ and $`b=0.0015/\overline{g}`$, where $`\overline{g}=10^3`$ is the average NMDA conductance.


text:
# Infrared cutoff dependence of the critical flavor number in three-dimensional QED

## Abstract

We solve, analytically and numerically, a gap equation in parity invariant QED3 in the presence of an infrared cutoff $`\mu `$ and derive an expression for the critical fermion number $`N_c`$ as a function of $`\mu `$. We argue that this dependence of $`N_c`$ on the infrared scale might solve the discrepancy between continuum Schwinger-Dyson equations studies and lattice simulations of QED3.

Parity invariant quantum electrodynamics in 2+1 dimensions with $`N`$ flavors of four-component massless fermions (QED3) QED3  has been attracting a lot of interest for almost two decades. While it was often regarded as a nice polygon for studying nonperturbative phenomena in gauge field theories, such as dynamical mass generation, recently the model has found applications in condensed matter physics, in particular in high-$`T_c`$ superconductivity supercond .

QED3 is ultraviolet finite and has a built-in intrinsic mass scale given by the dimensionful gauge coupling $`e`$, or $`\alpha =e^2N/8`$, which plays a role similar to the $`\mathrm{\Lambda }`$ scale parameter in QCD. In the leading order in $`1/N`$ expansion, it was found that at large momenta ($`p\alpha `$) the effective coupling between fermions and gauge bosons vanishes (asymptotic freedom) whereas it has a finite value or infrared stable fixed point at $`p\alpha `$ apbokawij86 . This behavior was shown to be robust against the introduction of higher order $`1/N`$ corrections nash ; guhare01 . Studies of the gap equation for a fermion dynamical mass in the leading order in the $`1/N`$ expansion have shown that massless QED3 exhibits chiral symmetry breaking whenever the number of fermion species $`N`$ is less than some critical value $`N_c`$, which is estimated to be in the region $`3<N_c<5`$ ($`N_c=32/\pi ^23.2`$ in the simplest ladder approximation apnawij88 ). A renormalization group analysis gives $`3<N_c<4`$ terao . Below such a critical $`N_c`$, the $`U(2N)`$ flavor symmetry is broken down to $`U(N)\times U(N)`$, the fermions acquire a dynamical mass, and $`2N^2`$ Goldstone bosons appear; for $`N>N_c`$ the particle spectrum of the model consists of interacting massless fermions and a photon. In addition, it was argued that the dynamical symmetry breaking phase transition at $`N=N_c`$ is a conformal phase transition miya97 ; gumish98 . The last one is characterized by a scaling function for the dynamical fermion mass with an essential singularity apnawij88 ; miya97 .

On the other hand, there is an argument due to Appelquist et al. Appelquist  that $`N_c3/2`$. The argument is based on the inequality $`f_{\mathrm{IR}}f_{\mathrm{UV}}`$ where $`f`$ is the thermodynamic free energy which is estimated in both infrared and ultraviolet regimes by counting massless degrees of freedom.

The above-described version of QED3, the so-called parity invariant noncompact version, was studied on a lattice lattice1 . Recent lattice simulations of QED3 with $`N2`$ have found no decisive signal for chiral symmetry breaking lattice2 . In particular, for $`N=2`$, Ref. lattice2  reports an upper bound for the chiral condensate $`\overline{\psi }\psi `$ to be of order $`5\times 10^5`$ in units of $`e^4`$. We recall that for quenched, $`N=0`$, QED3 both numerical quenched\_num  and analytical dim\_reg\_qed  studies have shown that chiral symmetry is always broken and numerically $`\overline{\psi }\psi 5\times 10^3`$. The above result of Ref. lattice2  seemingly contradicts studies based on Schwinger-Dyson (SD) equations which advocate dynamical symmetry breaking for $`N=2`$ and seems to favor the estimate of Appelquist et al. (we mention also that work in progress is being done on lattice simulations of QED3 with a single, $`N=1`$, fermion flavor lattice3 ).

One of the major problems in studying dynamical symmetry breaking using lattice simulations is to obtain control of finite size effects, which play a nontrivial role due to the presence of a massless photon lattice2 . In terms of the intrinsic lattice spacing $`a`$, three dimensionless length scales appear in lattice simulations: the size of the lattice $`L1050`$, the dimensionless lattice coupling constant $`\beta 1/(e^2a)`$, and the bare fermion mass $`m_0a`$. In order to establish dynamical symmetry breaking the role of all these scales in the problem should be well under control. At present, lattice discretization appears to be well understood; however, the finite size effects appear nontrivial and are most likely the source of the discrepancy between lattice and continuum studies of QED3 lattice2  .

In this paper, we present a possible explanation for the discrepancy between recent lattice and continuum studies of QED3 by studying analytically and numerically the gap equation in the presence of the infrared (IR) cutoff $`\mu `$, where $`\mu `$ is inversely related to the size $`L`$ of a lattice. Since the characteristic ratio $`e^2/\mu `$ for continuum QED3 turns out to be proportional to the ratio $`L/\beta `$ in lattice simulations lattice2 , a comparison between two approaches can be made. We will show that the presence of an IR cutoff reduces the value of the critical number $`N_c`$ and derive the relationship between $`N_c`$ and $`\mu `$. Recently, a similar gap equation, but with a massive photon, was studied numerically in Refs. cheng ; pefr03  in connection with applications of QED3 to high-$`T_c`$ superconductors (for an earlier numerical study in QED3, see Ref. kona92 ).

The gap equation of massless QED3 with IR cutoff $`\mu `$ reads (compare with Eq. (3.3) of Ref.apbokawij86 )

$`\mathrm{\Sigma }(p)={\displaystyle \frac{\lambda \alpha }{2p}}{\displaystyle _\mu ^{\mathrm{}}}𝑑k{\displaystyle \frac{k\mathrm{\Sigma }(k)}{k^2+\mathrm{\Sigma }^2(k)}}\mathrm{ln}{\displaystyle \frac{k+p+\alpha }{|kp|+\alpha }},`$ (1)
where $`\lambda =8/N\pi ^2`$ and the Landau gauge is used. Equation (1) is the simplest approximation to the SD equation for the fermion self-energy which neglects corrections to the fermion wave-function renormalization and vertex corrections. Extensive studies showed that the nature of chiral symmetry breaking which emerges from Eq. (1) with $`\mu =0`$ remains qualitatively unchanged after including higher order corrections. In what follows we solve Eq. (1) numerically but in order to be able to treat it analytically one needs to make further approximations. Because the integrand is damped at large momenta of integration, the main contribution comes from the region with momenta $`k\alpha `$; thus, expanding the logarithm we come to the simplified gap equation

$`\mathrm{\Sigma }(p)={\displaystyle \frac{\lambda }{p}}{\displaystyle _\mu ^\alpha }𝑑k{\displaystyle \frac{k\mathrm{\Sigma }(k)}{k^2+\mathrm{\Sigma }^2(k)}}\mathrm{min}(k,p).`$ (2)
The scale $`\mu `$ can be identified with the inverse size of the lattice in the temporal direction as $`\mu \pi /La`$ footnote  since both scales $`\mu `$ and $`\pi /La`$ represent a minimal fermion momentum in continuum and lattice theories, respectively \[$`\mu `$ is related to the photon mass $`m_a`$ ($`\mu =m_a^2/\alpha `$) in the model of Ref. pefr03 \]. Although this argument for the identification of $`\mu `$ and $`\pi /La`$ (up to a factor of order 1) is certainly heuristic, it is not unreasonable.

On the other hand, the ultraviolet (UV) cutoff $`\alpha `$ in Eq. (2) is in fact a physical scale which separates a nonperturbative dynamics at $`k\alpha `$ from a perturbative one, with momenta much larger than $`\alpha `$, where the dimensionless running coupling $`\overline{\alpha }(k)=\alpha /(k+\alpha )`$ apbokawij86  is weak. Neglecting those perturbative contributions in Eq. (1) results in an effective UV cutoff $`\alpha `$ in Eq. (2).

We shall solve the above equation (2) using the well-established linearized approximation when the momentum-dependent dynamical mass in the denominator $`k^2+\mathrm{\Sigma }^2(k)`$ is replaced by the constant dynamical mass at the lower limit: $`\mathrm{\Sigma }^2(k)\mathrm{\Sigma }_0^2`$, $`\mathrm{\Sigma }_0=\mathrm{\Sigma }(\mu )`$. This approximation is known to work well, especially near the phase transition point where $`\mathrm{\Sigma }_0\alpha `$. As we show, the linearized equation (2) with an IR cutoff can be dealt with analytically and an exact expression for $`N_c`$ as a function of $`\alpha /\mu `$ can be derived, following the method developed in Refs. MGS ; RNC ; guhataya02 .

In the linearized approximation, Eq. (2) can be reduced to a hypergeometric differential equation

$`u(1u)\mathrm{\Sigma }^{\prime \prime }(u)+{\displaystyle \frac{3}{2}}(1u)\mathrm{\Sigma }^{}(u){\displaystyle \frac{\lambda }{4}}\mathrm{\Sigma }(u)=0,`$ (3)
with the IR and UV boundary conditions

$`\left[p^2{\displaystyle \frac{d\mathrm{\Sigma }}{dp}}\right]_{p=\mu }=0,\left[p{\displaystyle \frac{d\mathrm{\Sigma }}{dp}}+\mathrm{\Sigma }\right]_{p=\alpha }=0,`$ (4)
where $`u=p^2/\mathrm{\Sigma }_0^2`$. The general solution of Eq. (3) can be written as

$`{\displaystyle \frac{\mathrm{\Sigma }(p)}{\mathrm{\Sigma }_0}}=c_1u_1(p)+c_2u_2(p),`$ (5)
where we choose

$`u_1(p)=F({\displaystyle \frac{1+i\nu }{4}},{\displaystyle \frac{1i\nu }{4}};{\displaystyle \frac{3}{2}};{\displaystyle \frac{p^2}{\mathrm{\Sigma }_0^2}}),`$
$`u_2(p)=\left({\displaystyle \frac{p}{\mathrm{\Sigma }_0}}\right)^{(1+i\nu )/2}`$
$`\times `$ $`F({\displaystyle \frac{1i\nu }{4}},{\displaystyle \frac{1+i\nu }{4}};1+{\displaystyle \frac{i\nu }{2}};{\displaystyle \frac{\mathrm{\Sigma }_0^2}{p^2}})+\text{c.c.},`$ (7)
as a particular pair of independent solutions of Eq. (3), with $`\nu =\sqrt{4\lambda 1}`$, and where $`F`$ is the hypergeometric function. If the infrared scale is set equal to zero ($`\mu =0`$), then the IR boundary condition gives the constraint $`c_2=0`$, since the function $`u_2`$ is too singular at $`p=0`$. Then, the UV boundary condition fixes the relationship between the dynamical mass $`\mathrm{\Sigma }_0`$, $`\nu `$, and $`\alpha `$. For nonzero $`\mu `$, the boundary conditions (4) determine the mass spectrum, giving rise to the equation

$`f(\mathrm{\Sigma }_0/\alpha ,\mu /\alpha ,\nu )=A_1B_2A_2B_1=0,`$ (8)
where

$`A_i\left[p{\displaystyle \frac{du_i}{dp}}+u_i\right]_{p=\alpha },B_i\left[p^2{\displaystyle \frac{du_i}{dp}}\right]_{p=\mu }.`$ (9)
By making use of the formulas for differentiating the hypergeometric function erdelyi  the real functions $`A_i`$ and $`B_i`$ are expressed as

$`A_1`$ $`=`$ $`F({\displaystyle \frac{1+i\nu }{4}},{\displaystyle \frac{1i\nu }{4}};{\displaystyle \frac{1}{2}};{\displaystyle \frac{\alpha ^2}{\mathrm{\Sigma }_0^2}}),`$ (10)
$`A_2`$ $`=`$ $`\left({\displaystyle \frac{\alpha }{\mathrm{\Sigma }_0}}\right)^{(1+i\nu )/2}{\displaystyle \frac{(1i\nu )}{2}}`$
$`\times `$ $`F(1{\displaystyle \frac{1i\nu }{4}},{\displaystyle \frac{1+i\nu }{4}};1+{\displaystyle \frac{i\nu }{2}};{\displaystyle \frac{\mathrm{\Sigma }_0^2}{\alpha ^2}})+\text{c.c.}`$
and

$`B_1`$ $`=`$ $`\mu {\displaystyle \frac{\mu ^2}{\mathrm{\Sigma }_0^2}}{\displaystyle \frac{(1+\nu ^2)}{12}}`$ (12)
$`\times `$ $`F(1+{\displaystyle \frac{1+i\nu }{4}},1+{\displaystyle \frac{1i\nu }{4}};{\displaystyle \frac{5}{2}};{\displaystyle \frac{\mu ^2}{\mathrm{\Sigma }_0^2}}),`$
$`B_2`$ $`=`$ $`\mu \left({\displaystyle \frac{\mu }{\mathrm{\Sigma }_0}}\right)^{(1+i\nu )/2}{\displaystyle \frac{(1+i\nu )}{2}}`$ (13)
$`\times `$ $`F(1+{\displaystyle \frac{1+i\nu }{4}},{\displaystyle \frac{1i\nu }{4}};1+{\displaystyle \frac{i\nu }{2}};{\displaystyle \frac{\mathrm{\Sigma }_0^2}{\mu ^2}})+\text{c.c.}`$
$`=`$ $`\mu \left({\displaystyle \frac{\mu }{\mathrm{\Sigma }_0}}\right)^{1/2}\sqrt{1+\nu ^2}`$
$`\times `$ $`\left[\text{Re}(F(\nu ))\mathrm{cos}\beta _2+\text{Im}(F(\nu ))\mathrm{sin}\beta _2\right],`$
where $`\beta _2=(\nu /2)\mathrm{ln}(\mu /\mathrm{\Sigma }_0)\mathrm{tan}^1\nu ,`$ and $`F(\nu )`$ is shorthand notation for the hypergeometric function:

$`F(\nu )=F(1+{\displaystyle \frac{1+i\nu }{4}},{\displaystyle \frac{1i\nu }{4}};1+{\displaystyle \frac{i\nu }{2}};{\displaystyle \frac{\mathrm{\Sigma }_0^2}{\mu ^2}}).`$ (14)
Since we study the critical behavior, we can always assume that $`\mathrm{\Sigma }_0\alpha `$, and therefore we can use the asymptotic expressions for $`A_1`$ and $`A_2`$:

$`A_1`$ $``$ $`\left({\displaystyle \frac{\alpha }{\mathrm{\Sigma }_0}}\right)^{1/2}\left[{\displaystyle \frac{\mathrm{\Gamma }(1/2)\mathrm{\Gamma }(i\nu /2)}{\mathrm{\Gamma }^2((1+i\nu )/4)}}\left({\displaystyle \frac{\alpha }{\mathrm{\Sigma }_0}}\right)^{i\nu /2}+\text{c.c.}\right]`$ (15)
$`=`$ $`\left({\displaystyle \frac{\alpha }{\mathrm{\Sigma }_0}}\right)^{1/2}|c(\nu )|\sqrt{1+\nu ^2}\mathrm{cos}\left(\alpha _2+\theta \right)`$
and

$`A_2`$ $``$ $`\left({\displaystyle \frac{\alpha }{\mathrm{\Sigma }_0}}\right)^{1/2}\sqrt{1+\nu ^2}\mathrm{cos}\alpha _2,`$ (16)
where

$`\alpha _2={\displaystyle \frac{\nu }{2}}\mathrm{ln}{\displaystyle \frac{\alpha }{\mathrm{\Sigma }_0}}+\mathrm{tan}^1\nu ,\theta =\mathrm{arg}c(\nu ),`$
$`c(\nu )={\displaystyle \frac{\mathrm{\Gamma }(3/2)\mathrm{\Gamma }(i\nu /2)}{\mathrm{\Gamma }((1+i\nu )/4)\mathrm{\Gamma }((5+i\nu )/4)}}.`$
By making use of (2.10.2) of erdelyi , we can write $`B_1`$ in a form similar to $`B_2`$:

$`B_1`$ $`=`$ $`\mu \left({\displaystyle \frac{\mu }{\mathrm{\Sigma }_0}}\right)^{(1+i\nu )/2}{\displaystyle \frac{(1+i\nu )}{2}}c(\nu )`$ (17)
$`\times `$ $`F(1+{\displaystyle \frac{1+i\nu }{4}},{\displaystyle \frac{1i\nu }{4}};1+{\displaystyle \frac{i\nu }{2}};{\displaystyle \frac{\mathrm{\Sigma }_0^2}{\mu ^2}})+\text{c.c.},`$
$`=`$ $`\mu \left({\displaystyle \frac{\mu }{\mathrm{\Sigma }_0}}\right)^{1/2}\sqrt{1+\nu ^2}|c(\nu )|[\text{Re}(F(\nu ))`$
$`\times `$ $`\mathrm{cos}(\beta _2+\theta )+\text{Im}(F(\nu ))\mathrm{sin}(\beta _2+\theta )].`$
Finally, using the above expressions (12)–(17), the gap equation (8) can be expressed as

$`f`$ $`=`$ $`\mu |c|(1+\nu ^2)\left({\displaystyle \frac{\alpha \mu }{\mathrm{\Sigma }_0^2}}\right)^{1/2}\mathrm{sin}\theta [\text{Re}(F(\nu ))`$ (18)
$`\times `$ $`\mathrm{sin}(\beta _2\alpha _2)\text{Im}(F(\nu ))\mathrm{cos}(\beta _2\alpha _2)].`$
One can convince oneself that for $`\mu \mathrm{\Sigma }_0`$ the last equation is reduced to

$$\mathrm{cos}\left(\frac{\nu }{2}\mathrm{ln}\frac{\alpha }{\mathrm{\Sigma }_0}+\theta +\mathrm{tan}^1\nu \right)=0,$$
(19)
which gives the well-known result for the dynamical mass exhibiting the essential singularity at $`\nu 0`$:

$$\mathrm{\Sigma }_0=\alpha \mathrm{exp}\left[\frac{2(\pi n+\pi /2\theta \mathrm{tan}^1\nu )}{\nu }\right].$$
(20)
We recall that only the solution with $`n=1`$ corresponds to the stable ground state. Note also that for $`N`$ close to $`N_c`$ ($`\nu 0`$) the dynamically generated mass $`\mathrm{\Sigma }_0`$ is much less than the scale $`\alpha `$, providing a hierarchy of mass scales in the model under consideration.

On the other hand, for $`\mathrm{\Sigma }_0\mu `$, we can use the power expansion of $`F(\nu )`$ to get the equation for the dynamical mass $`\mathrm{\Sigma }_0`$. By expanding $`F(\nu )`$ of Eq. (14) for $`\mathrm{\Sigma }_0\mu `$, we obtain

$`\text{Re}(F(\nu ))=1+{\displaystyle \frac{\mathrm{\Sigma }_0^2}{\mu ^2}}\rho (\nu ^2)\mathrm{cos}\varphi ,`$ (21)
$`\text{Im}(F(\nu ))={\displaystyle \frac{\mathrm{\Sigma }_0^2}{\mu ^2}}\rho (\nu ^2)\mathrm{sin}\varphi ,`$ (22)
where

$`\rho (\nu ^2)={\displaystyle \frac{1}{8}}\sqrt{{\displaystyle \frac{(25+\nu ^2)(1+\nu ^2)}{4+\nu ^2}}},`$ (23)
$`\varphi =\mathrm{tan}^1{\displaystyle \frac{\nu (13+\nu ^2)}{102\nu ^2}}.`$ (24)
Thus Eq. (18) reduces to

$`{\displaystyle \frac{\mathrm{\Sigma }_0^2}{\mu ^2}}={\displaystyle \frac{1}{\rho (\nu ^2)}}{\displaystyle \frac{\mathrm{sin}\left(\frac{\nu }{2}\mathrm{ln}\frac{\alpha }{\mu }+2\mathrm{tan}^1\nu \right)}{\mathrm{sin}\left(\frac{\nu }{2}\mathrm{ln}\frac{\alpha }{\mu }+2\mathrm{tan}^1\nu +\varphi \right)}}.`$ (25)
Since, for $`0<\nu <1`$ $`\varphi `$, is negative, we find that a nontrivial or real solution for $`\mathrm{\Sigma }_0`$ arises when $`\nu `$ exceeds the critical value $`\nu _c`$ determined by the equation

$`{\displaystyle \frac{\nu _c}{2}}\mathrm{ln}{\displaystyle \frac{\alpha }{\mu }}+2\mathrm{tan}^1\nu _c=\pi ,`$ (26)
so that the $`\mathrm{sin}`$ in the numerator of Eq. (25) becomes negative (for $`\nu >\nu _c`$). Note that the form of the gap equation (18) is different in two regions $`\mu \mathrm{\Sigma }_0`$ and $`\mu \mathrm{\Sigma }_0`$: while in the first one ($`\mu \mathrm{\Sigma }_0`$) we observe oscillations in the mass variable, in the second one ($`\mu \mathrm{\Sigma }_0`$) such oscillations disappear. This is reflected also in the character of the mass dependence on the coupling constant \[compare Eqs. (20) and (29) below\]. In general, it can be shown that Eq. (18) has $`n`$ nontrivial solutions where the number $`n`$ is given by

$`n=\left[{\displaystyle \frac{\nu }{2\pi }}\mathrm{ln}{\displaystyle \frac{\alpha }{\mu }}+{\displaystyle \frac{2}{\pi }}\mathrm{tan}^1\nu \right],`$ (27)
and the symbol $`[C]`$ denotes the integer part of the number $`C`$. For small $`\nu _c`$, Eq. (26) gives

$`\nu _c={\displaystyle \frac{\pi }{2+(1/2)\mathrm{ln}(\alpha /\mu )}}.`$ (28)
A similar functional dependence of $`N_c`$ was guessed in Ref. pefr03  by fitting their numerical study of the gap equation with a massive photon.

Near the critical point, we find from Eq. (25) the mean-field scaling law for $`\mathrm{\Sigma }_0`$,

$`{\displaystyle \frac{\mathrm{\Sigma }_0^2}{\mu ^2}}=\sigma (\nu \nu _c),`$ (29)
where the positive “nonuniversal” constant of proportionality, $`\sigma `$, is

$`\sigma ={\displaystyle \frac{1}{\rho (\nu _c^2)}}\left[{\displaystyle \frac{(\pi 2\mathrm{tan}^1\nu _c)}{\nu _c}}+{\displaystyle \frac{2}{\nu _c^2+1}}\right]{\displaystyle \frac{1}{\mathrm{sin}(\pi +\varphi _c)}},`$
with $`\varphi _c`$ the value of $`\varphi `$ at $`\nu =\nu _c`$.

We now compare the analytically obtained critical value of $`\nu _c`$ (Eq. (26)) or $`N_c`$ as a function of the ratio $`\pi e^2/\mu `$ ($`\alpha =Ne^2/8`$) with the numerical solution of the integral equation (1) which was obtained by using numerical methods described in Ref. guhare96 . Basically, starting with some initial guess for $`\mathrm{\Sigma }(p)`$, the integral equation is iterated using splines on a logarithmic momentum scale, until sufficient precision is reached. Both the analytical solution for Eq. (2) and numerical solution for Eq. (1) are depicted in Fig. 1. The analytical solution of Eq. (2) agrees reasonable well with the numerical solution of Eq. (1) for values of $`N1`$. However, Eq. (2) (hence Eq. (26)) provides a rather poor description of dynamical symmetry breaking for small values of $`N`$ ( actually it does not allow quenched limit in contrast to Eq.(1)). The bending of the curve for the analytic solution for $`0N1`$ is due to the linear dependence of the effective UV cutoff $`\alpha `$ on $`N`$. Furthermore, the numerical solution shows the existence of a critical ratio $`\pi e^2/\mu 40`$ for $`N=0`$ (i.e., quenched QED3) below which no chiral symmetry breaking occurs.

From Fig. 1 it can be extracted that in order to find chiral symmetry breaking for $`N=2`$ at least a ratio $`\pi e^2/\mu 5\times 10^3`$ is required. Note that in order to draw decisive conclusions on chiral symmetry breaking on a lattice, the volume of the lattice $`L^3`$ must be large not just in dimensionless units but also compared to any dynamically generated correlations in the system. In particular, the physical volume $`(L/\beta )^3`$ is required to be large lattice2 . Typical lattice simulations use at the most a ratio $`L/\beta 50100`$ and appear to lie outside the region for dynamical symmetry breaking. Indeed, by associating $`La\pi /\mu `$ and $`e^21/(\beta a)`$, we have $`L/\beta \pi e^2/\mu `$; thus, dynamical symmetry breaking near $`N=2`$ requires a ratio $`L/\beta 5\times 10^3`$. This explains the absence of any signs of chiral symmetry breaking in lattice simulations with $`N=2`$. For $`N=1`$ the ratio is $`L/\beta 240`$ for the numerical solution of Eq. (1); therefore, we believe that lattice simulations with a single fermion flavor are crucial in establishing whether there is a critical number $`N_c`$ in QED3 and how this value depends on the size $`L`$ of the lattice. Moreover, as is evident from Eq. (29), the scaling near the critical point is expected to be of the mean-field type and such a scaling can also be checked using lattice simulations.

We are grateful to V.A. Miransky for fruitful discussions and comments and S. Hands for valuable correspondence. We thank A.H. Hams for help in the numerical programs. V.P.G. acknowledges support from the Natural Sciences and Engineering Research Council of Canada. The research of V.P.G. has been also supported in part by the SCOPES-projects 7 IP 062607 and 7 UKPJ062150.00/1 of Swiss NSF and by NSF Grant No. PHY-0070986.