Neural field theory of synaptic plasticity

Journal of Theoretical Biology 285 (2011) 156–163

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Neural field theory of synaptic plasticity P.A. Robinson a,b, a b

School of Physics, University of Sydney, Sydney, NSW 2006, Australia Brain Dynamics Center, Sydney Medical School - Western, University of Sydney, Westmead, NSW 2145, Australia

a r t i c l e i n f o

a b s t r a c t

Article history: Received 9 January 2011 Received in revised form 16 June 2011 Accepted 17 June 2011 Available online 18 July 2011

Plasticity is crucial to neural development, learning, and memory. In the common in vivo situation where postsynaptic neural activity results from multiple presynaptic inputs, it is shown that a widely used class of correlation-dependent and spike-timing dependent plasticity rules can be written in a form that can be incorporated into neural field theory, which enables their system-level dynamics to be investigated. It is shown that the resulting plasticity dynamics depends strongly on the stimulus spectrum via overall system frequency responses. In the case of perturbations that are approximately linear, explicit formulas are found for the dynamics in terms of stimulus spectra via system transfer functions. The resulting theory is applied to a simple model system to reveal how collective effects, especially resonances, can drastically modify system-level plasticity dynamics from that implied by single-neuron analyses. The simplified model illustrates the potential relevance of these effects in applications to brain stimulation, synaptic homeostasis, and epilepsy. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Plasticity Modeling Neural field theory STDP Neural systems

1. Introduction Neural plasticity involves activity-driven changes to strengths of synapses between neurons that alter neural connectivity (Hebb, 1949; Abbott and Nelson, 2000; Dayan and Abbott, 2001). It is essential to development of connectivities for information processing, learning, and memory, with analogs in artificial neural nets (Abbott and Nelson, 2000; Dayan and Abbott, 2001). Large-scale modification of synaptic strengths may also be responsible for the phenomenon of kindling, in which repeated epileptic seizure induction lowers the seizure threshold (Engel and Pedley, 1997; Wasserman et al., 2008; Swartz, 2009), and for some effects of transcranial magnetic stimulation (TMS) and electroconvulsive therapy (Wasserman et al., 2008; Swartz, 2009), and in applications such as control of synchrony in distributed systems (Lubenov and Siapas, 2008). Plasticity is thus very important in neural systems, but understanding at the systems level has been impeded because no existing plasticity theory spans from synapses to large-scale dynamics. The core aims of this paper are thus to show how this gap can be bridged using neural field theory (NFT) and that the results imply that emergent systems-level effects are critical to understanding plasticity in the brain as a whole.

 Correspondence address: School of Physics, University of Sydney, Sydney, NSW 2006, Australia. E-mail address: [email protected]

0022-5193/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2011.06.023

Plasticity has long been studied on the basis that synapses are modified by neural activity. Hebb (1949) postulated that synapses strengthen when pre- and postsynaptic activity are positively correlated, implying a connection between them, as encapsulated in the dictum ‘‘neurons that fire together wire together.’’ Plasticity has since explained many phenomena, especially in development and learning (Abbott and Nelson, 2000; Dayan and Abbott, 2001). In correlation-dependent plasticity (CDP) synaptic strengthening depends only on the time difference between pre- and postsynaptic activity, regardless of order. Usually spike-timing-dependent plasticity (STDP) strengthens a synapse (termed long-term potentiation or LTP) when a postsynaptic spike (i.e., an action potential) follows a presynaptic one within a certain temporal window (Abbott and Nelson, 2000; Dayan and Abbott, 2001), consistent with causality, and weakens it if the order is reversed (termed long-term depression or LTD). Most analyses of plasticity have concentrated on spike-based interactions at individual synapses in cases where a postsynaptic spike is caused by one or a few presynaptic ones. In reality, a typical cortical neuron receives inputs from  4000 others each generating  10 spikes/s (Braitenberg and Schuz, 1991; Robinson et al., 2004; Deco et al., 2008), so each postsynaptic spike results from dynamics driven by  4000 presynaptic ones—i.e., 40 000 presynaptic spikes per second drive 10 postsynaptic ones. In reality, the effective ratio is somewhat less than this (but still very large) because the typical 10–20 ms time constant of the combined synaptic, dendritic, and soma-voltage responses means that it is only the  4002800 presynaptic spikes within this typical time window that determine the soma voltage and hence

P.A. Robinson / Journal of Theoretical Biology 285 (2011) 156–163

the postsynaptic firing response (Deco et al., 2008). Thus the commonly assumed picture of each postsynaptic spike being generated by one or a few presynaptic spike(s) can be very misleading when populations are involved—multineuron effects must also be incorporated, particularly network oscillations involved in brain rhythms, which correlate large-scale activity (Steriade et al., 1990; Nunez, 1995; Robinson et al., 2002; Robinson, 2003; Fries, 2005; Busza´ki, 2006). Work on rate-based plasticity and its relationship to singleneuron STDP shows that the former reproduces many features of the latter if the analysis preserves temporal correlations (Kempter et al., 1999; Abbott and Nelson, 2000; Dayan and Abbott, 2001; these references mostly concentrated on single neurons, but the key point that correlations must be preserved applies more widely). Hence, and because of the large-scale nature of many applications of plasticity mentioned above, we seek an approach that will enable rate-based plasticity rules to be analyzed in a way that incorporates systems-level effects. Neural field theory (NFT) provides a promising framework to study plasticity in neural populations and at the systems level because it averages over the dynamics of many individual neurons (see Deco et al., 2008 for a review of NFT) to track the large-scale dynamics of populations. Although it is might be thought that averaging would remove the ability to incorporate relative timings of pre- and postsynaptic activity and that a spikebased description would thus be essential to study any aspect of STDP, we show below that this is not true in general and that NFT can thus be used to study the average effects of rate-based STDP dynamics on neural populations. This is especially significant because NFT can easily include multiscale spatial and temporal features in the dynamics and has successfully predicted a wide variety of linear and nonlinear phenomena, including time series, spectra, correlations, evoked responses, and seizures (Beurle, 1956; Wilson and Cowan, 1973; Lopes da Silva et al., 1974; Freeman, 1975; Nunez, 1974, 1995; Jirsa and Haken, 1996; Rennie et al., 2002; Robinson, 2003, 2005; Robinson et al., 1997, 2002, 2004, 2008; Steyn-Ross et al., 1999, 2005; Breakspear et al., 2006; Hutt et al., 2003; Coombes, 2005; Deco et al., 2008), including analysis of an example of CDP as a mechanism for memory erasure in a purely cortical system (Steyn-Ross et al., 2005). Overall, the above points highlight that activity dependent plasticity is a crucial component of neuronal dynamics that underlies many aspects of self-organization in learning and memory. This paper analyzes the dynamics and frequency dependence of this plasticity by deriving an expression for the rate of change of expected (average) connection strengths in terms of the frequency content of external inputs to a neuronal system, and (in the linear limit) the system’s underlying transfer functions. This enables one to predict whether connection strengths will increase or decrease given the frequency content of the inputs to which the system is exposed. The analysis relies on noting that average plasticity responses depend upon correlations between presynaptic inputs and postsynaptic responses. This enables one to derive the expected correlation, given the transfer functions mapping input to responses. These transfer functions are estimated from established sets of differential equations describing neuronal responses to input that can be linearized around their fixed point (for some assumed connection strengths). One can then derive equations for the changes in expected connection strength in the frequency domain to see how connection strengths will change given inputs with different frequencies or spectra. The structure of this paper is as follows: In Section 2 we adopt a rate-based plasticity rule from the literature, generalize it to span both CDP and STDP, and briefly outline how it arises from spike-based rules while clarify the assumptions implicit in the

157

necessary averaging. One could equally well use a range of other rules that involve spike timings (Bienenstock et al., 1982; Abbott and Nelson, 2000; Dayan and Abbott, 2001; Izhikevich and Desai, 2003), but STDP and CDP are selected here because they are among the most basic and widely studied. It is next shown how to write our rule in terms of frequency responses in the case that activity changes can be approximated as perturbations to a steady state. The resulting expressions are then incorporated into NFT in the context of a simplified model system that includes the key biological effects of synaptic and dendritic dynamics, threshold firing responses, axonal propagation within and between neural populations, stimulus inputs, and internal feedbacks. For simplicity we restrict attention to cases where the activity can be approximated as spatially uniform, yielding a formulation that is similar to neural mass approximations, except that it retains characteristic propagation timescales. Exploration of these results in Section 3 shows that network interactions strongly modify plasticity relative to the single-synapse cases, sometimes even reversing its sign to convert LTP into LTD and vice versa. These changes are elucidated in the frequency domain, where it is shown that system-wide resonances can dominate the dynamics through the frequency-dependent correlations they induce between pre- and postsynaptic activity. Analysis of the model system also further illustrates that these effects have relevance to plasticity phenomena implicated in phenomena such as TMS, synaptic homeostasis, and seizure kindling, simplified versions of which serve as examples. Finally, Section 4 discusses the main findings and their relevance to applications to more detailed corticothalamic networks, also noting aspects that can be generalized to treat other plasticity rules and systems.

2. Theory In this section we set down the relevant plasticity rule from the literature and briefly outline its derivation and underlying assumptions. We then re-express it in a spectral form that gives new insights into its dynamics. In Section 2.2, NFT is then introduced and used to obtain equations for plasticity dynamics in a simple exemplar model system.

2.1. Plasticity rule Many rules proposed for change dsab of the strength sab of a synapse to a neuron of type a from type b (sab o0 for inhibition) depend on relative timings of pre- and postsynaptic spikes and have the form X ð1Þ ðdsab Þj p Hðt 0j ti Þ, i

(Abbott and Nelson, 2000; Dayan and Abbott, 2001; Izhikevich and Desai, 2003), where ti and t 0j are the times of particular pre- and postsynaptic spikes, i and j, and H is a window function (H depends on a and b more generally). Although other rules are possible, (1) encompasses the most commonly used CDP and STDP rules and exemplifies the issues discussed in Section 1. In the cases of interest here, each postsynaptic spike in neuron a results from multiple incoming spikes, and interacts with the effects of all presynaptic spikes at their afferent synapses to produce plastic changes. The sum in (1) can then be written as an integral over spikes approximated as delta functions: Z 1 ðdsab Þj p dtHðtÞfab ðtÞ, ð2Þ 1

158

fab ðtÞ ¼

P.A. Robinson / Journal of Theoretical Biology 285 (2011) 156–163

X

dðtti Þ,

ð3Þ

i

where t ¼ t 0j t and fab represents presynaptic spikes. To determine the mean rate of change of the synaptic strength sab we multiply (2) by the postsynaptic spike rate Qa ðt 0 Þ, and average over a time window of characteristic width T around t, where T is longer than the timescales of the plasticity window H and of pre- and postsynaptic interspike intervals (ISIs), but shorter than synaptic evolution—i.e., T is a few seconds. One implementation of this averaging is to take an exponential moving average over a time T; Robinson et al. (2008) also used this approach to obtain the postsynaptic firing rate Qa from individual spikes. In practical terms there is no contradiction between this limited averaging and the infinite bounds on the t integral in (2) because the rapid decrease of the window function HðtÞ with jtj effectively cuts off the integral in (2) at a few tens of ms. A more rigorous treatment would amend (2) to have finite bounds, but this would make negligible difference to the following results for the reasons just given. The above steps yield the following equation for the evolution of sab on timescales of order T or longer: Z 1 d/sab S ¼B dt/Qa ðt þ tÞHðtÞfab ðtÞS, ð4Þ dt 1 where B is a constant (SI unit: V s2 ) and angle brackets denote the average over the window of width T. It is important to note that we average the product Qa Hfab , not the separate factors, to preserve correlations between the rates Qa and fab . Because we are interested here long-term plasticity, rather than in phenomena that occur on very short timescales, we average over timescales longer than those of H and pre- and postsynaptic ISIs, (i) hundreds of incoming spikes at multiple synapses are typically involved in producing a given postsynaptic spike, as discussed in Section 1 and (ii) the relative timing t 0j ti is broadly distributed so the result does not depend significantly on individual spike timings. Hence, the time difference between pre- and postsynaptic spikes can be approximated as continuous in (4) without appreciably changing the result (this is equivalent to approximating an integral by a discrete sum of values of the integrand, multiplied by the time differences between those values). Averaging is thus equivalent to approximating both the spike trains fab and Qa as smoothed instantaneous quantities, which is valid a fortiori because we henceforth interpret (4) as being further averaged over many neurons to become the population mean. The quantity t is then the time difference between each change in presynaptic activity and any resulting postsynaptic activity change, as determined by neural dynamics (see Section 2.2). Eq. (4) differs from the response to constant average rates because /Qa ðt þ tÞfab ðtÞSa /Qa ðt þ tÞS/fab ðtÞS:

ð5Þ

We thus stress that, despite averaging, NFT can treat STDP effects, provided that the additional step of replacing the left of (5) by the right is not taken. In Section 2.2, we introduce NFT and use it to evaluate (4). Large-scale normal brain activity is characterized by a steady state plus perturbations, a view that has been borne out by numerous investigations, although pathological large-scale nonlinear dynamics can also occur and be treated by NFT—e.g., in seizures (Robinson et al., 2002; Breakspear et al., 2006; Deco et al., 2008). Note that this does not contradict the existence of highly nonlinear single- or few-neuron dynamics at the small scale; e.g., the dynamics of spike generation itself. In the steady state, /dsab =dtS ¼ 0 for self-consistency, and net plasticity is due to perturbations. Hence, in neural populations we reinterpret Qa and fab in (4) as perturbations from mean values

[denoted ð0Þ], which can be caused by naturally occurring or artificial stimuli, for example. We write these perturbations as ð0Þ dQa ðtÞ ¼ Qa ðtÞQað0Þ and dfab ðtÞ ¼ fab ðtÞfab and note that the time average of the perturbations is zero (any part that did not average to zero would be incorporated in the mean value). Writing the (perturbed) quantities on the right of (4) as Fourier transforms, and evaluating the average (i.e., the zero-frequency component) by integrating over t, gives Z 1 d/sab S B n ¼ dodQa ðoÞHn ðoÞdfab ðoÞ, ð6Þ dt 2p 1 where asterisks denote complex conjugates. Note that the ð0Þ product of the two zeroth-order terms Qað0Þ fab must give no net contribution in order to avoid contradiction with these being steady-state values. Similarly, terms linear in the perturbations average to zero, because the perturbations themselves do. At this point the perturbations have not been assumed to be small, although we make this approximation below in Section 2.2. Note that, in deriving (6), we formally extend the bounds of integration to infinity, but this makes little difference to the final results because the form of the window function HðoÞ and the other factors in (6) restricts contributions to moderate o (see Section 3 for further discussion). In evaluating (6) in what follows we employ the widely used form (Abbott and Nelson, 2000) HðtÞ ¼ A þ expðt=tP Þ, HðtÞ ¼ A expðt=tP Þ,

t 4 0,

ð7Þ

t r 0,

ð8Þ

where tP is the plasticity timescale [different values could be used in (7) and (8), but we set them equal for simplicity] and the A 7 are constants. This expression covers many possible forms of window, as discussed by Abbott and Nelson (2000), including all possible combinations of signs of CDP and STDP. Fourier transformation of (7) and (8) gives HðoÞ ¼

H0 þ iotP H1

ð9Þ

1 þ ðotP Þ2

with H0 ¼ ðA þ þA ÞtP and H1 ¼ ðA þ A ÞtP , where H1 ¼ 0 corresponds to CDP and H0 ¼ 0 yields STDP; for reference, the most common model of CDP, with HðtÞ ¼ dðtÞ, corresponds to A þ ¼ A ¼ 1=ð2tP Þ in the limit tP -0. Eq. (9) implies CDP tends to be more sensitive to the lowest-o contributions of the other factors in the integrand in (6) than STDP, which is preferentially driven by activity at o  1=tP . In nonlinear applications (e.g., seizures, ECT), (4) and (6) must be evaluated numerically, but other applications can involve perturbations small enough to use linear approximations that have yielded many experimentally verified NFT predictions, including steady states and second order quantities such as spectra and correlation functions (Robinson et al., 1997; Robinson, 2003; Deco et al., 2008). 2.2. Neural field theory and model system To evaluate (6) in terms of linear perturbations dQa and dfab of the rates, one must adopt a specific NFT and physical system in order to calculate these quantities. Here we use a recent NFT (Robinson, 2005), to analyze plasticity in the system in Fig. 1, which includes the key features of input fax to neurons a from

x

ax

a

aa

Fig. 1. Simplified model system with internal (a) and external (x) populations and a feedback loop showing pulse densities faa and fax .

P.A. Robinson / Journal of Theoretical Biology 285 (2011) 156–163

external neurons x, and internal feedback faa with delay t0. This idealized system incorporates plasticity dynamics at external and internal synapses has feedbacks and delays like those in real systems (Robinson, 2005; Deco et al., 2008), and is sufficient to identify the main types of emergent systems-level dynamics. For simplicity we omit spatial dependences. NFT averages neural properties over scales of  0:1 mm, containing many neurons (Robinson, 2005; Deco et al., 2008). The soma potential Va responds to afferent spikes via synaptic, dendritic, and soma dynamics, with an overall characteristic time constant 1=a. The resulting response approximately obeys  2 d aþ Va ðtÞ ¼ a2 ½naa faa ðtt0 Þ þ nax fax ðtÞ, ð10Þ dt where nab ¼ Nab sab , and Nab is the mean number of connections per cell of type a from cells of type b (Robinson et al., 1997; Robinson, 2005; Deco et al., 2008). Here we measure Va relative to resting. Action potentials are produced at the axonal hillock at a rate Qa when Va exceeds the firing threshold. However, this threshold is not identical for all neurons, owing to differences in their individual properties and environments, including fluctuating depolarization levels, thereby leading to its smearing out in a neural population. An approximate population average response function is then the sigmoid Qmax , ð11Þ 1 þexp½CfVa ðtÞyg=s pffiffiffi where C ¼ p= 3, Qmax is the maximum firing rate, and y and s are the mean threshold and its standard deviation (Wilson and Cowan, 1973; Freeman, 1975; Robinson, 2005; Deco et al., 2008; Marreiros et al., 2008). Mean spike rates in axons, faa and fax , propagate approximately as if governed by damped wave (Jirsa and Haken, 1996; Robinson et al., 1997; Robinson, 2005; Deco et al., 2008). In the purely temporal case considered here, the spatial dependence is not considered and one has  2 d gab þ fab ðtÞ ¼ g2ab Qb ðtÞ, ð12Þ dt Qa ¼

where gab ¼ vab =rab , vab is the axonal velocity, and rab is the characteristic axonal range, where b¼a for activity within population a and b¼x for external inputs. The left side of (12) describes the temporal decay of spatially uniform activity in the absence of regeneration due to neural firing, while the right side represents the activity that is the source of axonal spikes. Eq. (13) does not have a spatial dependence of activity, so no spatial derivatives appear, as they would in a general neural field equation (which would have position as an additional argument), but it is not quite equivalent to a neural mass model, which would also have gab ¼ vab =rab ¼ 1 for all populations, corresponding to the signal transit time across the system being negligible. We keep finite gaa to enable the effects of delays within population a to be retained but, since we do not model sensory systems explicitly here, we set fax ¼ Qx (i.e., gax ¼ 1), which involves no loss of generality because any form of fax can still be specified. Eq. (10) implies that the soma voltage Va [and hence the firing rate Qa, via (11)] lags the terms on the right that drive it by the characteristic time 1=a required for a soma response. Hence, this NFT equation encapsulates phase differences between the rates Qa, faa , and fax that reflect relative timings of neural inputs and outputs and is closely related to spike timings in individual neurons (Robinson et al., 2008). Similarly (10) and (12) express the relative timings imposed by propagation within (timescale 1=g) and between (timescale t0) populations (Robinson et al., 2002; Robinson, 2005). Thus NFT retains the information on

159

relative timings required to evaluate (6) when many neurons are involved. The apparent paradox of how correlations in spike timings can be explored by a theory that does not represent spikes explicitly can be further resolved by the simple observation that when rates of spiking Qa and fab are both high (perhaps with a time offset due to dynamics), spikes are closer together in time and their timings are positively correlated. This is analogous to the situation pertaining to pressure waves in a gas—when pressure is high, molecules can be inferred to be close together, and to collide more often, even when they are not represented explicitly in the theory (e.g., hydrodynamics) that describes the waves. The NFT equations (10) and (12), with (11), constitute a closed set of equations for the voltages and firing rates. These have steady states obtained by setting all time derivatives to zero, eliminating Va, and solving the resulting equation ðs=CÞ ln½Qmax =Qað0Þ 1 þ ynaa Qað0Þ nax Qxð0Þ ¼ 0:

ð13Þ

In a more general corticothalamic model, the relevant generalization of (13) has been found to have a stable steady-state low-Qað0Þ solution at low Qxð0Þ , which is interpreted as representing the baseline of normal activity, and which yields firing rates in accord with experiment (Robinson et al., 1997, 2002, 2004; Deco et al., 2008). Eq. (13) has a similar low-firing rate solution and we interpret in the same way. The next step in evaluating (6) is to express the perturbations in rates dQa and dfab in terms of the external steady-state signal Qxð0Þ and its perturbations dQx . We linearize relative to the system steady state discussed in the previous paragraph by setting ra ¼ dQ a =dV a (evaluated at Qað0Þ ) and approximating the synaptic strengths sab as being constant on the timescale of fluctuations in dQa and dfab —i.e., by making an adiabatic approximation. The linear transfer functions dQa =dQx and dfab =dQx can then be evaluated straightforwardly by eliminating other perturbation variables from the linearized equations. These steps yield (omitting arguments o of some functions for brevity) Z d/sab S B 1 ¼ doKab ð1þ o2 tP2 Þ1 PðoÞ, ð14Þ dt p 0 where Kaa ¼ Re½ðH0 iotP H1 ÞGnaa ,

ð15Þ

n n Þ=Xax , Kax ¼ Re½ðH0 iotP H1 Þð1Xaa

ð16Þ

Gaa ¼ ð1io=gaa Þ2 ,

ð17Þ

Xaa ¼ Gaa Leiot0 Gaa ,

ð18Þ

Xax ¼ LGax ,

ð19Þ

L ¼ ð1io=aÞ2 ,

ð20Þ

  Xax dQx ðoÞ2  , PðoÞ ¼  1Xaa 

ð21Þ

PðoÞ is the spectrum of dQa , Gab ¼ ra nab , and dQx ðoÞ in (21) is the Fourier component of the external signal dQx . We find Kaa ¼ ½H0 ð1o2 =g2aa Þ2o2 tP H1 =gaa =ð1 þ o2 =g2aa Þ2 ,

ð22Þ

Kax ¼ ½H0 ð1o2 =a2 Þ þ 2o2 tP H1 =a=Gax ,

ð23Þ

where (23) is written for Gaa ¼ 0, the general result being too algebraically complicated to yield much insight beyond what can be gained from the form (16). Eq. (14) is the key result of the paper. It expresses the rate of change of expected connection strengths sab ð ¼ saa ,sax Þ in terms of an integral over the frequency spectrum of external inputs,

160

P.A. Robinson / Journal of Theoretical Biology 285 (2011) 156–163

have substantial effects, modifying plasticity dynamics from the forms found in single- or few-neuron analyses. Fig. 2(a) illustrates the integrands in (14) for gain jGaa j 5 1 (so there is negligible feedback via the loop shown in Fig. 1) vs. frequency f ¼ o=ð2pÞ. Fig. 2(b) shows the corresponding integrals up to an upper frequency bound fmax as a function of fmax . We interpret these results by first noting that (9) implies that the window function jHðoÞj peaks at o  1=tP (where tP is the plasticity time scale), which is also a low-pass cutoff, while jGaa j and jLj in (17) and (21) peak at o ¼ 0, with cutoffs set by the axonal rate gaa and the dendritic rate a, respectively. Hence, the integrands decrease rapidly as o increases beyond the lesser of oc ¼ minfa, gaa ,1=tP g and the bandwidth of dQx (  1 in Fig. 2)—i.e., beyond the lowest of the relevant low-pass cutoffs. This is seen in Fig. 2(a) where CDP integrands peak at low f and have a rapid fall off above about 5 Hz. Similarly, the magnitudes of the STDP integrands fall off steeply above 10–20 Hz. Eqs. (22) and (23) predict sign changes in the factors Kaa and Kax at o ¼ gaa , a, respectively, and these are confirmed by the results in Fig. 2(a). These sign changes occur when synaptodendritic delays cause the postsynaptic response to a sinusoidal wave peak to lie just before the next peak, reversing their effective order. In STDP Fig. 2(a) shows that jKaa j and jKax j increase from 0 at f¼0 with Kax 40 and Kaa o 0 for all f. Positive Kax accords with single-neuron analyses, since feedback is negligible for gains

modulated by the factors Kab (b ¼ x,a), which incorporate the system transfer functions. By assuming simple forms for the spectral profile of inputs we can use (14) to estimate the ensuing change in connection strengths.

3. Results We now show how (14) implies the emergence of systemslevel effects on plasticity that are especially clear in the frequency domain. In Sections 3.2–3.4 we then briefly outline how some of these effects potentially apply to analysis of plasticity in applications to TMS, sleep homeostasis, and seizure kindling, as illustrations of their wider importance. We stress that detailed verification of these possible applications will require analysis of more complex and realistic corticothalamic networks, so the present discussion aims only to demonstrate their potential relevance. 3.1. General features Equations (14), (22) and (23) reveal that plasticity dynamics can have strong dependences on stimulus frequency. In particular, the possibility exists that the denominator in (21) may vanish or become very small, meaning that system-level resonances may

0.008

0.3 aa

0.006 0.2

0.002

Integral

aa

0.004

Integrand

ax

ax

ax

0.1

0.0

0.000

aa

ax

aa -0.002 0

10

20

30

-0.1

40

0

10

20

f (Hz)

30

40

fmax (Hz)

0.02 0.03

aa

aa 0.02

ax

ax

Integrand

Integrand

0.01

0.00

-0.01

aa

0.01 ax

ax

0.00 aa -0.01

-0.02 0

5

10 f (Hz)

15

20

0

5

10 f (Hz)

15

20

Fig. 2. (a) Integrands in (14) vs. f and (b) integrals up to fmax vs. fmax for CDP (solid) and STDP (dots). Both have jQ ðoÞj ¼ 1, a ¼ 80 s1 , gaa ¼ 120 s1 , tP ¼ 0.01 s, t0 ¼ 0, Gaa ¼ 0, Gaa ¼ 0:01, Gax ¼ 1, A þ ¼ 1, and A ¼ 1 (CDP) or A ¼ 1 (STDP). (c) and (d) are as for (a) but with STDP multiplied by 3, t0 ¼ 0.1 s, Gaa ¼ 0:95 in (c), and Gaa ¼ 0:7 in (d). Temporal parameters have typical values from plasticity literature and from Deco et al. (2008).

P.A. Robinson / Journal of Theoretical Biology 285 (2011) 156–163

o  mp=ðt0 þ 2=aÞ,

ð24Þ

where the integer m is even for Gaa 40 and odd for Gaa o 0, in accord with theory (c.f., Robinson et al., 2002). Furthermore we see that they can change the effective order of Qa and fab , thereby reversing the signs of the integrands. We find numerically that this effect strengthens if jGaa j is larger, implying stronger feedback. The effect is also stronger if the plasticity timescale and delay satisfy tP 5 t0 , implying sharper temporal localization of feedbacks. When specialized to the present system, previous corticothalamic linear stability analyses imply that it is unstable for gains Gaa 41 (o ¼ 0 instability via saddle-node bifurcation) and Gaa t1 [3.5 Hz instability at m¼1 in (24) via a supercritical Hopf bifurcation], and is stable for 1 tGaa o 1 (Robinson et al., 2002; Breakspear et al., 2006). At Gaa ¼ 0:95, Fig. 2(c) shows that the m¼ 0, 2 resonances (0, 7 Hz) greatly enhance CDP, both absolutely and relative to STDP, the latter being little affected by the m ¼0 resonance at o ¼ 0. In the case of negative feedback with Gaa ¼ 0:9 [see Fig. 2(d)] the m ¼1 resonance has similar, weaker effects. The feedback loop thus causes every synapse to be affected by activity in every population, as seen in TMS experiments that examine the effects on other neural areas or populations of pulses applied to a specific area or population (Wasserman et al., 2008, Chapter 12). 3.2. Frequency dependences and TMS Fig. 3 shows regimes of LTP and LTD for Kaa and Kax from (15) and (16) vs. o and H1 =H0 , the latter ratio measuring the importance of STDP relative to CDP. It is seen that Kaa and Kax yield very different dynamics, depending on o and the plasticity parameters, H0 and H1. Both move from LTD to LTP as o increases, provided H0 40 and H1 t 0:5 for Kaa, and provided H0 4 0 and H1 \ 0:5 for Kax. This is consistent with data on neocortical neurons, which can display LTD or LTP (Abbott and Nelson, 2000; Dayan and Abbott, 2001; Tononi and Cirelli, 2003). The above frequency dependence also concurs with observations that show that TMS at 1 Hz induces LTD, causing temporary cortical lesions that are unresponsive to stimuli, whereas 5–20 Hz TMS causes LTP and increased excitability, with 5 Hz being marginal (Wasserman et al., 2008). Experimental frequency bounds on this crossover [via (22) and (23), as plotted in Fig. 3] are consistent with the ones predicted from estimated dendritic and axonal cortical parameters a, g ¼ 50150 s1 (Robinson et al., 2002, 2004), provided 1 tjH1 =H0 j (i.e., provided STDP exceeds CDP in magnitude). TMS thus constrains these parameters via NFT, which may thus help to provide new ways to quantify them and TMS effects more generally. 3.3. Sleep homeostasis Tononi and Cirelli (2003) hypothesized that synapses are plastically strengthened during wake, and weakened during slow wave sleep (SWS; dominated by activity at 0.5–4.5 Hz) to achieve daily homeostatic regulation of total synaptic number and

6 4 S

W

2 H1/H0

jGaa j 5 1, so postsynaptic activity Qa always follows presynaptic fax . Paradoxical values Kaa o0 occur because neither Qa nor faa can always be said to precede the other in the closed aa loop. Consistency with network propagation yields the overall net result that Qa causally precedes faa , even though faa also affects Qa. Resonances in the integrand in (14) occur where Xaa  1 [cf., (19)], so stimuli at these frequencies drive plasticity most effectively. For cases with significantly nonzero gain jGaa j, Fig. 2(c) and (d) show that these resonances occur at frequencies given by

161

Kax>0

0 -2

Kaa>0 W

S -4 -6 0

5

10 f (Hz)

15

20

Fig. 3. Regimes of LTP (above the solid curves for ax, below the dotted curves for aa) and LTD (the remaining range in each case) for H0 40 (LTP and LTD are reversed for H0 o 0). The pairs of curves shown indicate the most extreme possible positions of the transition curve for combinations of parameters in the ranges tP ¼ 0:0120:02 s, a ¼ 50280 s1 , and gaa ¼ 502120 s1 ; i.e., in each case the actual transition is predicted to occur at a single curve located somewhere between the pair of bounds shown. Vertical dashed and dot-dashed lines show the transition frequencies implied by TMS and sleep homeostasis, as discussed in Sections 3.2 and 3.3, respectively. Horizontal bars show examples of possible sleep–wake (S–W) paths that span the transition loci to allow homeostasis, as discussed in Section 3.3, where f is reinterpreted as an effective frequency of the whole EEG spectrum.

metabolic load. They also argued that this cyclic strengthening and pruning might assist memory consolidation [see also Diekelmann et al., 2011] and account for observed enhancements in task performance after sleep. Here we do not attempt to analyze the hypothesis of Tononi and Cirelli (2003) in detail, which would require study of a full corticothalamic model; rather, we show that our results for the highly simplified system of Fig. 1 are consistent with this picture. To do this, we integrate over o to find the net effect of the prevailing spectrum of activity fluctuations. We simplify interpretation by using the expression IðoÞ ¼

1 , ð1þ o2 =o20 Þð1 þ o2 =o2c Þ

ð25Þ

to approximate all the factors in (14) aside from the forms (15) and (16) of the Kab at low o. Substitution of (25) into (14) and evaluation of the resulting approximate integral then shows that pffiffiffiffiffiffiffiffiffiffiffiffiffi if o is replaced by oeff ¼ o0 oc , Fig. 3 also represents the behavior of the overall integral in (14) vs. feff ¼ oeff =ð2pÞ and the ratio of plasticity window parameters H1 =H0 (i.e., rather than showing the integrand vs. f and H1 =H0 , as labeled). To compare Fig. 3 with sleep and wake data, we note that typical bandwidths of EEG spectra imply o0  6 s1 in SWS, o0  20 s1 in wake, and oc ¼ 502100 s1 in both states (Niedermeyer and Lopes da Silva, 1999). So feff  2:524 Hz in SWS and 5–7 Hz in wake. Fig. 3 shows two examples of horizontal bars linking sleep (S) and wake (W) states, at H1 =H0 ¼ 7 3 for H0 40. The upper bar shows that the sleep–wake cycle can carry the ax synapses back and forth between LTP in wake and LTD in sleep, provided H1 is positive for those synapses. The lower bar shows that the same cycle can carry the aa synapses between the same regimes in the same states provided H1 is negative for those synapses (Abbott and Nelson, 2000, for example, noted that a range of possible plasticity windows can occur, corresponding to all combinations of signs of H1 and H0). The vertical dot-dashed

162

P.A. Robinson / Journal of Theoretical Biology 285 (2011) 156–163

lines in the figure show these effective frequency limits at all values of H1 =H0 , demonstrating that these bounds on the transition frequency from net LTD to net LTP agree with the theoretical transition frequency (i.e., they bracket the possible location of the theoretical transition frequency, which lies between the pairs of curved lines in the figure—see the caption for further explanation) provided 1 t jH1 =H0 j. This is the same bound as found from TMS considerations in Section 3.2. 3.4. Seizure kindling, TMS, and ECT The occurrence of LTP at f \5 Hz, discussed in Section 3.2, is also qualitatively consistent with kindling, where repeated seizure induction lowers the seizure threshold (Engel and Pedley, 1997). In previous work, the m ¼2 resonance from (24) in a more detailed corticothalamic model was identified with 10 Hz tonic-clonic seizures (Robinson et al., 2002; Breakspear et al., 2006; Deco et al., 2008). Induced 10 Hz seizures (which need nonlinear analysis, so present results are only qualitative) would thus plausibly increase jGaa j, leaving the system more prone to seizure induction thereafter. These seizure findings are also consistent with findings regarding ECT and TMS therapies for clinical depression, which likely involve plasticity (Wasserman et al., 2008; Swartz, 2009). TMS therapy is applied via repetitive TMS pulses, and is found to be most effective at an approximately 10 Hz repetition rate, while ECT induces therapeutic 10 Hz seizures via electrical stimulation. Both therapies plausibly increase synaptic strengths via plasticity, thereby raising brain activity from depressed levels. The present work implies that this optimal frequency arises because (i) stimuli must be above about 5 Hz to yield LTP, (ii) stimuli must be below about 15 Hz to escape low-pass filtering, and (iii) within the remaining 5–15 Hz range the corticothalamic resonance at the alpha frequency of approximately 10 Hz enhances plasticity, as illustrated in Fig. 2(c). These results have the interesting implication that the optimal frequency for TMS therapy should be close to each individual’s electroencephalographic alpha frequency.

4. Summary We have generalized a plasticity rule to include CDP and STDP and have expressed it in a spectral form suitable for evaluation using neural field theory. NFT analysis in the context of a simple exemplar system has elucidated systems-level emergent features of the plasticity equations and their potential for detailed application in a number of contexts. The main results are (i) It is shown that NFT can be used to evaluate the correlations between pre- and postsynaptic activity that are required by a class of plasticity rules that is sufficiently broad to cover the most commonly used types of CDP and STDP rules. This demonstrates that NFT can treat the average effects of STDP on neural populations, provided that inherent information on correlations between pre- and postsynaptic activity is not discarded. (ii) A spectral representation of plasticity is derived and its utility is demonstrated. These results reveal the strong roles of system feedbacks and resonances in plasticity, including reversing its sign under some circumstances (e.g., LTD becomes LTP and vice versa). Large-scale network effects must thus be considered when investigating plasticity, since real neurons are not isolated, but are embedded in networks, where they typically interact with thousands of others. The role found for system resonances is particularly important, as plasticity attributes new functional significance to them, and we speculate that this may help explain their prominence in real

brains. For example, Fries (2005) postulated that correlations induced by brain rhythms may be important in cognition, and the present results may imply similar links to memory formation. (iii) Applications of even our simplified model are consistent with key experimental results on brain stimulation via TMS, synaptic homeostasis, and seizure kindling, implying potential applicability to more realistic systems. In particular, the frequency of transition between LTD and LTP is correctly predicted to be circa 5 Hz, consistent with sleep and wake lying on opposite sides of the transition, as in the synaptic homeostasis hypothesis. The existence of an optimal frequency for TMS therapy near 10 Hz is also consistent with experiment, and NFT implies links to the mechanisms of alpha-rhythm generation and seizure kindling. The above results make the application of NFT to more detailed cortical, corticothalamic, and other systems a priority. As well as the applications to sleep, epilepsy, and brain stimulation mentioned above, future applications may also include studies of learning, memory, and TMS effects on cognition, where plasticity is also critical. The same techniques can also be applied to incorporate other types of plasticity rules into NFT.

Acknowledgments The author thanks J. Roberts for drawing Fig. 1. The Australian Research Council and the Westmead Millennium Institute supported this work. References Abbott, L.F., Nelson, S.B., 2000. Synaptic plasticity: taming the beast. Nat. Neurosci. Suppl. 3, 1178–1183. Beurle, R.L., 1956. Properties of a mass of cells capable of regenerating pulses. Philos. Trans. Roy. Soc. Lond. Ser. B 240, 55–94. Bienenstock, E.L., Cooper, L.N., Munro, P.W., 1982. Theory for the development of neuron selectivity: orientation specificity and binocular interaction in visual cortex. J. Neurosci. 2, 32–48. Braitenberg, V., Schuz, A., 1991. Anatomy of the Cortex: Statistics and Geometry. Springer, Berlin. Breakspear, M., Roberts, J.A., Terry, J.R., Rodrigues, S., Mahant, N., Robinson, P.A., 2006. A unifying explanation of primary generalized seizures through nonlinear brain modeling and bifurcation analysis. Cerebral Cortex 16, 1296–1313. Busza´ki, G., 2006. Rhythms of the Brain. Oxford University Press, Oxford. Coombes, S., 2005. Waves, bumps, and patterns in neural field theories. Biol. Cybern. 93, 91–108. Dayan, P., Abbott, L.F., 2001. Theoretical Neuroscience. MIT, Cambridge. Deco, G., Jirsa, V.K., Robinson, P.A., Breakspear, M., Friston, K., 2008. The dynamic brain: from spiking neurons to neural masses and cortical fields. Pub. Lib. Sci. Computat. Biol. 4,e1000092. ¨ Diekelmann, S., Buchel, Born, J., Rasch, B., 2011. Labile or stable: opposing consequences for memory when reactivated during waking and sleep. Nat. Neurosci. 14, 381–386. Engel, J., Pedley, T.A. (Eds.), 1997. Epilepsy: A Comprehensive Textbook. Lippincott-Raven, Philadelphia. Freeman, W.J., 1975. Mass Action in the Nervous System. Academic, New York. Fries, P., 2005. A mechanism for cognitive dynamics: neuronal communication through neuronal coherence. Trends Cogn. Sci. 9, 474–480. Hebb, D.O., 1949. The Organization of Behavior: A Neuropsychological Theory. Wiley, NY. Hutt, A., Bestehorn, M., Wennekers, T., 2003. Pattern formation in intracortical neuronal fields. Network: Comput. Neural Syst. 14, 351–368. Izhikevich, E.M., Desai, N.S., 2003. Relating STDP to BCM. Neural Comput. 15, 1511–1523. Jirsa, V.K., Haken, H., 1996. Phys. Rev. Lett. 77, 960–963. Kempter, R., Gerstner, W., van Hemmen, J.L., 1999. Hebbian learning and spiking neurons. Phys. Rev. E 59, 4498–4514. Lopes da Silva, F.H., Hoeks, A., Smits, H., Zetterberg, L.H., 1974. Model of brain rhythmic activity. The alpha-rhythm of the thalamus. Kybernetik 15, 27–37. Lubenov, E.V., Siapas, A.G., 2008. Decoupling through synchrony in neuronal circuits with propagation delays. Neuron 58, 118–131. Marreiros, A.C., Daunizeau, J., Kiebel, S.J., Friston, K.J., 2008. Population dynamics: variance and the sigmoid activation function. NeuroImage 42, 147–157. Niedermeyer, E., Lopes da Silva, F., 1999. Electroencephalography, fourth ed. Lippincott Williams & Wilkins, Philadelphia. Nunez, P.L., 1974. The brain wave equation: a model for the EEG. Math. Biosci. 21, 279–297.

P.A. Robinson / Journal of Theoretical Biology 285 (2011) 156–163

Nunez, P.L., 1995. Neocortical Dynamics and Human EEG Rhythms. Oxford University Press, New York. Rennie, C.J., Robinson, P.A., Wright, J.J., 2002. Unified neurophysical model of EEG spectra and evoked potentials. Biol. Cybern. 85, 457–471. Robinson, P.A., 2003. Neurophysical theory of coherence and correlations of electroencephalographic and electrocorticographic signals. J. Theor. Biol. 222, 163–175. Robinson, P.A., 2005. Propagator theory of brain dynamics. Phys. Rev. E 72, 011904. Robinson, P.A., Rennie, C.J., Rowe, D.L., 2002. Dynamics of large-scale brain activity in normal arousal states and epileptic seizures. Phys. Rev. E 65, 041924. Robinson, P.A., Rennie, C.J., Rowe, D.L., O’Connor, S.C., 2004. Estimation of multiscale neurophysiologic parameters by electroencephalographic means. Hum. Brain Mapp. 23, 53–72. Robinson, P.A., Rennie, C.J., Wright, J.J., 1997. Propagation and stability of waves of electrical activity in the cerebral cortex. Phys. Rev. E 56, 826–840. Robinson, P.A., Wu, H., Kim, J.W., 2008. Neural rate equations for bursting dynamics derived from conductance-based spike equations. J. Theor. Biol. 250, 663–672.

163

Steriade, M., Gloor, P., Llina´s, R.R., Lopes da Silva, F.H., Mesulam, M.M., 1990. Basic mechanisms of cerebral rhythmic activities. Electroenceph. Clin. Neurophysiol. 76, 481–508. Steyn-Ross, M.L., Steyn-Ross, D.A., Sleigh, J.W., Liley, D.T.J., 1999. Theoretical electroencephalogram stationary spectrum for a white-noise driven cortex: evidence for a general anesthetic-induced phase transition. Phys. Rev. E 60, 7299–7311. Steyn-Ross, M.L., Steyn-Ross, D.A., Sleigh, J.W., Wilson, M.T., Wilcocks, L.C., 2005. Proposed mechanism for learning and memory erasure in a white-noisedriven sleeping cortex. Phys. Rev. E 72, 1–11. Swartz, C.M., 2009. Electroconvulsive and Neuromodulation Therapies. Cambridge University Press, Cambridge. Tononi, G., Cirelli, C., 2003. Sleep and synaptic homeostasis: a hypothesis. Brain Res. Bull. 62, 143–150. Wasserman, E.M., Epstein, C.M., Ziemann, U., Walsh, V., Paus, T., Lisanby, S.H., 2008. The Oxford Handbook of Transcranial Stimulation. Oxford University Press, Oxford. Wilson, H.R., Cowan, J.D., 1973. A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik 13, 455–480.