It was previously reported, that temperature may significantly influence neural dynamics on the different levels of brain function. Thus, in computational neuroscience, it would be useful to make models scalable for a wide range of various brain temperatures. However, lack of experimental data and an absence of temperature-dependent analytical models of synaptic conductance does not allow to include temperature effects at the multi-neuron modeling level. In this paper, we propose a first step to deal with this problem: A new analytical model of AMPA-type synaptic conductance, which is able to incorporate temperature effects in low-frequency stimulations. It was constructed based on Markov model description of AMPA receptor kinetics using the set of coupled ODEs. The closed-form solution for the set of differential equations was found using uncoupling assumption (introduced in the paper) with few simplifications motivated both from experimental data and from Monte Carlo simulation of synaptic transmission. The model may be used for computationally efficient and biologically accurate implementation of temperature effects on AMPA receptor conductance in large-scale neural network simulations. As a result, it may open a wide range of new possibilities for researching the influence of temperature on certain aspects of brain functioning. Keywords Temperature · AMPA-receptor · Modeling of synapses · Analytical modeling · Uncoupling set of ODEs 1 Introduction understanding of temperature effects on different levels of brain function may be useful in developing more sophisti- From a medical perspective, it has been suggested that tight cated methods of treatment for different neurological disor- control of brain temperature in patients, suffering during a ders which are sensitive to temperature (including hot-water epilepsy (Kowacs et al. 2005), autism (Helt et al. 2008)or post-traumatic period is highly recommended (Shigemori et al. 2012). However, despite the fact that the techniques brain injury (Mrozek et al. 2012; Dietrich 1992). of the control of the brain temperature have been devel- On the level of single neurons, multiple effects of oped, direct mechanisms of the influence of temperature on temperature on brain function have been observed. The most neural dynamics are still uncertain (Badjatia 2009). Better important from the perspective of neural dynamics are: 1) Temperature influences membrane resting potential (Hodgkin and Huxley 1952; Buzatu 2009). 2) Temperature affects ion channels dynamics (Hille Action Editor: Upinder Singh Bhalla 2001; Sterratt 2015). 3) Temperature affects synaptic transmission (Asztely Dominik S. Kufel et al. 1997; Weight and Erulkar 1976;Schiffand email@example.com Somjen 1985). The Polish Children’s Fund, Pasteura 5a, Temperature effects on membrane resting potentials (the 02-093 Warsaw, Poland Goldman-Hodgkin-Katz equation) and on ion-channel Department of Physics & Astronomy, University College dynamics (e.g. Hodgkin and Huxley (1952)) are now rela- London, 0 Gower St, WC1E 6BT, London, UK tively well-characterized. But the influence of temperature Faculty of Mathematics, Physics and Computer Science, on synaptic transmission has proven to be more difficult to Institute of Computer Science, Department model (e.g. De Schutter et al. (2009)). This may be because of Neuroinformatics, Maria Curie-Sklodowska University in Lublin, Akademicka 9, 20-033 Lublin, Poland of the various processes involved in synaptic transmission, 380 J Comput Neurosci (2018) 44:379–391 such as presynaptic release of neurotransmitter; dynamics if temperature coefficients in one kinetic scheme were of a vesicle pore; diffusion of neurotransmitter; binding of found, they would be invalid for other schemes (unless the neurotransmitter; and kinetics of postsynaptic receptors. one finds a way to link different kinetic schemes, The influence of temperature on each of this processes is which is currently not possible apart from linking very different, and they combine to the significant modification simple kinetic models (Shelley and Magleby 2008)). of the synaptic transmission (Fuxe et al. 2005). In fact, both of the approaches described above are similar, Therefore, generally, in the creation of biologically-real as amplitude and time constants in phenomenological models in computational neuroscience, it is useful to make modeling (under certain assumptions) may be interpreted them easily scalable for different temperatures. This is as a combination of different kinetic rates (Destexhe et al. especially important because some of the neurobiological 1994b). experiments - for example, in vitro studies - are usually Generally, the problem of including temperature effects conducted at temperatures lower than physiological ones. in synapse modeling is complex and yet not well- Therefore, a better knowledge about the temperature characterized. Both the first and second approaches dependence of synaptic transmission may be crucial for are hard to generalize for different phenomenological linking in vitro and in vivo studies. Furthermore, an optimal functions describing synaptic conductance or for different way of including a full description of temperature effects in microphysiological kinetic schemes. Including temperature neural simulations may allow computational brain research effects, require additional neurobiological experiments, to address new areas: for example, investigation of the which does not allow previously developed models to be influence of temperature on such neurological disorders like easily scalable for a wide range of brain temperatures. In hot-water epilepsy, cerebral brain injury or autism. this paper, a novel approach to the problem of including It is possible to tackle the problem of including temperature effects on modeling synapses is proposed. On temperature effects on synapses from different perspectives: the basis of previous experimental and numerical research, we construct assumptions for a new analytical model (1) A first approach is to include some multiplicative to include temperature effects in modeling the kinetics coefficients accounting for the influence of tempera- of the α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic ture in phenomenological synapse models e.g. alpha acid (AMPA) receptor. First, using Monte Carlo simulation function, dual-exponential functions, single exponen- and Markov modeling we propose simplifications of an tial functions. In this approach, it is required to experimental kinetic scheme (Postlethwaite et al. 2007) multiply all of these time constants and amplitudes to allow for a closed form solution of the set of ODEs of phenomenological functions by some (probably describing this problem. Second, we introduce the concept different) factors associated with temperature. These of uncoupling of the differential equation system describing coefficients would mimic the increase in the speed AMPA receptor kinetics. Third, after solving the resulting of chemical reactions with temperature (according set of differential equations, we compare results using the to Arrhenius equation). However, there is a signif- constructed model with numerical and experimental data. icant problem related to this approach. The values Finally, we suggest that our model provides a simple way to of temperature multiplicative factors are not known mimic temperature effects in neural dynamics simulations apriori (for an arbitrarily chosen kinetic scheme). at low frequencies, regardless of the phenomenological These values would have to be obtained for each function used to describe AMPA synaptic conductance. study separately by performing additional neurobio- logical experiments at different temperatures, which is usually not possible for in vitro research. 2 Methods and model (2) A second approach is to model synapses on the micro- physiological level - to investigate temperature effects In this section, we construct a new, analytical model of on the kinetics of synaptic receptor proteins, with con- the conductance of AMPA-type synapse for low-frequency formation dynamics described by kinetic schemes. To stimulations using uncoupling assumption of set of ODEs include temperature effects, it is required to multi- along with simplifications from experimental and numerical ply all of the kinetic rate constants between differ- data. ent conformational states by coefficients dependent on temperature. This approach was previously taken experimentally (Postlethwaite et al. 2007;Caisetal. Monte Carlo simulation of synaptic transmission Some of 2008). Nonetheless, the possible problem is that all the assumptions of the analytical model presented below of the temperature coefficients (which scale rate con- are based on the analysis of the data from Monte Carlo stants) are specific for given kinetic scheme. So, even simulation of the synaptic transmission. The simulation was J Comput Neurosci (2018) 44:379–391 381 Table 1 Parameters of the AMPA receptor model (see kinetic scheme analytical model are presented in the Table 1. The code in and text) 2 MCell is available. Parameter Description Value Analytical model Our model is based on the following k Agonist binding rate 10 [1/molar · 1/s] 3 numerical and experimental findings: k Agonist unbinding rate 8 · 10 [1/s] k Channel opening rate 20 · 10 [1/s] (1) Acceleration in postsynaptic AMPA receptor kinetics k Channel closing rate 10 · 10 [1/s] is the predominant effect of increased temperature k Desensitization rate 4 · 10 [1/s] on altered synaptic responses at low frequencies k Resensitization rate 15 [1/s] ((Postlethwaite et al. 2007)). With this assumption, −4 A Amplitude of glutamate 7.48 · 10 [molar] modeling of temperature effects on synapses was concentration simplified by considering temperature effects only on ω Decay time constant of 2471 [1/s] AMPA receptor kinetics, rather than also on modified glutamate concentration presynaptic release and/or neurotransmitter diffusion Q Coefficient of temper- 2.4 dynamics. Furthermore, we assumed that to include ature dependence of temperature effects on AMPA receptor kinetics it kinetic rates is sufficient to multiply all of the rate constants (k ,k ,k ,k ,k ) by a single temperature coefficient b o c d r Q (Postlethwaite et al. 2007). constructed based on the assumptions and parameters of (2) As suggested by Postlethwaite et al. (2007), temper- Postlethwaite et al. (2007) - code of their original simulation ature effects are mediated by driving AMPARs to is available using the MCell simulator (Stiles et al. 1996). higher sub-conductance states. To include higher sub- The most important assumptions of the simulation are as conductance states in an analytical model of AMPA follows: receptors, a few simplifications of the complex 13 1. The geometry used in the Monte Carlo simulation refers state, 30 transitions kinetic scheme of Postlethwaite to the AMPA-type synapse at the calyx of Held of a et al. (2007)(Scheme 1inFig. 1) were made. Scheme rat. The morphological data was given by Satzler ¨ et al. 1 was re-written into the simplified form of Scheme (2002). In the simulation, presynaptic and postsynaptic 2. This form uses the symmetry of states and transi- terminals were separated by the synaptic cleft of tions in Scheme 1. Based on uncoupling of equations, 28nm. Postsynaptic terminal consisted Postsynaptic described in point (4) below, this symmetry is evident Density (PSD) with an area of 0.32 micrometer on in that one can divide Scheme 1 into five orders of 0.32 micrometer, populated by 80 AMPA receptors. sub-conducting states. Additionally, four neighboring PSDs were included - (3) For simplification all of the state transitions except separated by 317 nm. the transition from closed to bound states in AMPA 2. Vesicle from which glutamate was released was a cube receptor kinetics (considering the single mesh of with the volume equal to the volume given by Satzler ¨ Scheme 2inFig. 2)wereassumedtobeMarkov et al. (2002) and connected to the synaptic cleft by a models: time and voltage independent and dependent gradually opening (with exponential dynamics) fusion only on the occupancy of neighboring states as was pore. Vesicle was released at variable locations above previously proposed by Destexhe et al. (1994b). a central postsynaptic density (PSD). Each vesicle Additional simplification was as follows: modifications contained 6000 glutamate molecules. in the directions of the transitions with rates k and k were r c 3. The diffusion rate of the glutamate was assumed to be made to rewrite Scheme 1 to Scheme 2 (these transitions are 6 2 equal to 3· 10 cm /s for the receptor kinetic parameters diagonal in Scheme 2 vs vertical in Scheme 1). This change detailed in Table 1. simplified the dynamics of each closed (C) state by reducing −6 All of the simulations used 1 · 10 s time step. For its coupling with the adjacent open and desensitized states more detailed discussion about the assumptions of the of the AMPA receptor. This approach allows for a much Monte Carlo simulation see Postlethwaite et al. (2007). https://senselab.med.yale.edu/modeldb/ShowModel.cshtml? All standard parameters used in the simulation and further model=239072 This approach is analogous to the way of including temperature effects on voltage-gated ion-channels (Hodgkin and Huxley 1952)and is motivated by the Arrhenius equation. Q10 is the multiplicative factor https://senselab.med.yale.edu/modeldb/showModel.cshtml? for increasing a given speed of reaction upon a 10 degree Celsius model=85981 increase of temperature. 382 J Comput Neurosci (2018) 44:379–391 easier analytical solution of the differential equations des- in the i-th order sub-conductance bound state (C ) cribing the simplified kinetic scheme. is perceived as 1 [generalization of assumption (3)]. Separately, the 1st order (and symmetrically, with However, we do introduce a scaling function to additional assumptions described below, the i-th order) of differentiate the absolute values of fractions in these the AMPAR kinetic states was assumed to behave as a single C states. For each i-th order of conductance λ (t ) i i mesh, according to Scheme 3 in Fig. 3. It was assumed that function is introduced. λ (t ) scales relative fraction of C0 1. This is the case when very few receptors bind channels in each state to an absolute (scaled identically glutamate so that nearly all receptors remain in form C . to all of the orders of the kinetic scheme) fraction. This is motivated by comparing the number of channels in This method allows us to uncouple set of twelve different states in the more detailed Monte Carlo simulation coupled differential equations with a complex formu- described above (see Fig. 4). Therefore, the fraction of lation to set of 4 pairs of differential equations (cou- channels in C state is considered always to be 1. In a pled only in pairs, rather than between different orders further analogous simplification, each state in the i-th mesh, of sub-conductance). Using this approach, we are able with i = 1, 2, 3, 4, is likewise assumed to remain at to include higher order sub-conductance states and, as 1. This assumption is necessary for an analytical model a result, find an analytical solution of the problem. because otherwise, the system of differential equations is (5) Glutamate binding was assumed to be independent, too complex to be solved analytically. Only the open states similarly to model by Robert and Howe (2003) of the AMPA receptor (O ,O ,O ,O ) contribute to the This departure from Postlethwaite et al. (2007)was 1 2 3 4 synaptic conductance (Smith et al. 2000). motivated by a disproportionate increase of the complexity (in comparison to the gain in accuracy) of (4) To describe the kinetic scheme with m states using an analytical solution of a set of differential equations differential equations it is necessary to write m − 1 when assuming cooperative binding. coupled differential equations, which complexity is (6) Glutamate concentration was assumed to be time- proportional to the number of transitions between dependent according to a single exponential decay different states (Destexhe et al. 1994b). With the function (Scimemi and Beato 2009), with parameters purpose to simplify this description we introduce the fitted to a simulation of glutamate concentration in assumption of uncoupling. In mathematical terms, the synaptic cleft (and consequently at the PSD) assumption of uncoupling the differential equations for in the Monte Carlo model of synaptic transmission. each single mesh may be written, for i-th order, as: Therefore, the function describing the binding rate, from closed to bound states, has a form dependent on k (t )(x +x ) k y +k z +k (x +x ) (1) b i−1 i c i+1 r i+1 u i+1 i the concentration of glutamate at the PSD: −ωt Particulary, for 1st order we obtain: k (t ) = k Ae (3) b b where the parameters ω [1/s] and A [molar] were fitted k (t )(x + x ) k y + k z + k (x + x ) (2) b 0 1 c 2 r 2 u 2 1 from averaging glutamate concentration (in the cleft) in Monte Carlo simulation of synaptic transmission. where x , z , y are the fractions of channels in a state i i i Therefore, the function describing the binding rate, from C , D , O respectively. i i i closed to bound states, has a form dependent on concentra- We see that from the perspective of the (i+1)-th order sub-conductance state, the fraction of channels tion of glutamate at the PSD: with a time constant ω and Fig. 1 Scheme 1 - modified kinetic scheme model by Postlethwaite et al. (2007), independent binding was assumed (see below) and no transitions between desensitized states (with minor influence on accuracy of results) J Comput Neurosci (2018) 44:379–391 383 Fig. 2 Scheme 2. Kinetic scheme used for construction of an analytical model consists five orders of subconductance (index numbers of states: 0,1,2,3,4) and four meshes (colored triangles) amplitude A fitted to the numerical data from the Monte (2007), and motivated by previous experimental work by Carlo simulation. Smith et al. (2000). This assumption suggests that we may αt A dual-exponential function [in the form of y = (Ae − break down the problem of finding the synaptic conductance βt Be )] was also used as a fitting method for glutamate into finding the sum of four functions (one for each order of concentration, but the resulting, more complex, system of sub-conductance, see Fig. 2). equations could not be solved analytically, and this increase Combining the above assumptions simplifies the original in complexity was not compensated for by a substantial complex kinetic scheme, containing 13 states and multiple gain in accuracy of the model in comparison to the single- transitions (described by 12 coupled differential equations). exponential function. Modification of directions of transitions (from O and D ) using the uncoupling concept splits the corresponding (7) According to experimental data, AMPA receptors system of differential equations to four, one-side dependent are tetramers (Rosenmund et al. 1998). Conductance and specular between themselves (see Fig. 2), single meshes of AMPA receptor can be described as a sum of as in Fig. 3. Furthermore, by assuming the fraction of conductances of all orders of subconductance multiplied channels in state C for i-th order of scheme to be equal to i−1 by different constants for different orders of states in 1 we are able to solve a coupled pair of differential equations kinetic scheme as suggested by Sahara and Takahashi (2001): n=4 g(t) = g a y (t ) (4) 4 i i i=1 where, g is a peak conductance of a channel in 4-fold bound state, n is a number of orders in kinetic scheme and a , y are scalling factor and fraction of open channels in i i the i-th state respectively. Conductances of different orders of states in kinetic scheme were set as fractions of the peak conductance at the 4-fold bound state O (O : a = 0.1, O : a = 0.4, O : a = 0.7, 4 1 1 2 2 3 3 O : a = 1.0) as was proposed by Postlethwaite et al. 4 4 Departure from model based on neurotransmitter concentration occurring as a pulse (described by Dirac Delta at t ) proposed by pulse Destexhe et al. (1994b) or Destexhe et al. (1994a) was motivated by availability of direct data of glutamate concentration on PSDs in Monte Carlo simulation and unit inconsistency problem (which is a result of modeling using Dirac Delta at t combined with the assumption Fig. 3 Scheme 3. Single mesh triangle from Scheme 2 described by an pulse about considering fraction of channels in C state as 1). independent pair of coupled differential equations 0 384 J Comput Neurosci (2018) 44:379–391 within every single mesh. Thus, simplification leads to 4 where S = k − ω, P = k + k + k − ω, R =−k + k + c d o u c d independent pairs of coupled differential equations. For the k + k . o u first order of conductance we know the solution explicitly As it is possible to be seen, the first-order approximating and for higher orders we use λ (t ) = x (t ) (which is the function is a sum of exponents (as suggested by Destexhe i i−1 solution in respect to the fraction of channels in state C et al. (1994b)). i−1 of pair of differential equation for (i-1)-th order), which The full solution for all orders of a kinetic scheme assures information flow from lower to higher orders of (Scheme 2) can be found in Appendix A. AMPAR sub-conductance. Making these assumptions, one obtains the general system of coupled linear ODEs, describing the 1st through 3 Results 4-th orders of kinetic states (Scheme 2): Analytical model of AMPA receptor We found that, for dx −ωt = k Ae λ (t ) − (k + k + k )x (t ) (5) kinetic rates fitted from the model of Postlethwaite et al. b i o u d i dt (2007), our analytical model is able to reproduce two of their key results from detailed Monte Carlo simulation. dy These results describe the dynamics of fractions of AMPAR = k x (t ) − k y (t ) (6) o i c i dt channels in different states (compare Fig. 5 here with Fig. (t ) =[O (t )] is a fraction of channels in an open 2B of Postlethwaite et al. (2007)) and the time courses of where y i i state of i-th order, x (t ) =[C (t )] is a fraction of channels AMPAR synaptic conductance at two different temperatures i i in a bound state of i-th order, λ (t ) is a function to convert (compare Fig. 6 here with Fig. 1A of Postlethwaite et al. (2007)). fraction of all channels to same absolute scale (not only relative for each order) - for the i-th order of Scheme 2 it equals to solution with respect to fraction of channels Fraction of channels Fractions of AMPAR channels in dif- in state C of two differential equations of (i-1)-th order: ferent states (Fig. 5) are both qualitatively and quantitatively λ (t ) = x (t ). Using this approach we may include higher i i−1 (after normalization: see Discussion) consistent with this sub-conductance states of AMPA receptor with an analytical obtained in Monte Carlo simulation (with a mean error 5% approach, due to the uncoupling of differential equations for the time courses of the open states). describing the kinetic scheme. Just after the t = 0, binding of glutamate leads AMPARs The above set of linear coupled ODEs [(5)and(6)] has to transition from unbound to bound states. The fraction a closed-form solution. For the first order, with boundary of channels in bound states is dependent on glutamate conditions y (0) = 0and x (0) = 0 the solution is: 1 1 concentration at PSDs and rate of unbinding in the AMPA Ak k Ak k Ak k receptor kinetic scheme. Excluding effects of unbinding, b o b o b o −ωt −(P +ω)t −k t y (t ) = e + e − e (7) with the diffusion of neurotransmitters in synaptic cleft SP RP RS Fig. 4 Fraction of AMPAR channels in unbound state C . Fraction of Fig. 5 Fraction of channels in different states as a function of time. channels decays from 1 to about 0.8 in 3ms, supporting our assumption “Open, scaled channels” refers the sum of the fractions of channels that the fraction of channels in state C in first mesh is far larger than in open states, with each fraction multiplied by its respective peak i−1 the fraction of channels in the other states conductance J Comput Neurosci (2018) 44:379–391 385 Scheme 1). Underestimation of the fraction of channels in desensitized states is due to the modification in directions of transitions for the first order of sub-conductance, smaller fraction of AMPARs is in a bound state. Resensitization (due to its low transition likelihood) has a minor influence on the results. Synaptic conductance The AMPAR conductance curve (Fig. 6) obtained from the analytical model is able to repro- duce (with 5% accuracy for the relative amplitude and peak time, which is within the experimental uncertainty range) the shape and scale of temperature effects on synaptic trans- mission (compared with Fig. 1B of Postlethwaite et al. (2007)). At 35 C both the rise and decay time constants of synaptic conductance are smaller in value (peak time is ◦ ◦ Fig. 6 Conductance curve of AMPA synapse in 25 C and 35 C. shorter). The peak conductance is larger (ratio about 1.25) ◦ ◦ Conductance in 35 C in comparison to 25 C has larger and quicker ◦ ◦ and is reached quicker in 35 C in comparison to 25 C. ◦ ◦ peak. Transition between 25 C and 35 C was achieved only by However, the analytical model predicts a too rapid rise- multiplication of all rate constants in the kinetic scheme by Q = 2.4 time of conductance in comparison to experimental data (see the time of peak on Fig. 6 and on Fig. 1B of considered as a random walk, the mean distance of a Postlethwaite et al. (2007)). This is due to the assumption single neurotransmitter from the location of release (vesicle (6) (see Methods). Namely, only a single exponential decay pore) should increase over time proportionally to N , time (from a peak value at t = 0) of glutamate concentration where N is the number of steps). Hence, the glutamate was used - which only roughly approximates reality (see concentration should decay with 1/ N , and the number Fig. 7). However, we did not include any other, more of bound states should increase proportionally to 1 − N . complex glutamate concentration functions as they did not However, as time passes some of the AMPAR channels allow for closed-form analytical solution to the differential unbind neurotransmitters (transitioning from C to C 1 0 equations system including higher sub-conductance states. state), therefore eventually reaching equilibrium - so the It was found that modifying the diffusion rates of the function of bound state fraction (in time) should be close to neurotransmitters by the temperature coefficient cannot a ‘flattened’ 1 − N . increase the accuracy of the solution (in comparison to During synaptic transmission, AMPAR protein under- the experimental data) - as presented in Fig. 8. Q diff goes conformational changes. The rate of these changes is coefficient of diffusion (effectively multiplying diffusion proportional to temperature - the effect of temperature is coefficient of glutamate) with temperature, leads to quicker reflected in our analytical model by scaling all of the rate decay of glutamate concentration in the synaptic cleft constants of the kinetic scheme by a Q parameter. The continuous growth of fraction of channels in desensitized states can be attributed to that the resensitizing rate of reac- tion is about three orders of magnitude smaller than the rate of desensitization. Therefore, channels after entering, are unlikely to leave desensitized states - the fraction of channels in desensitized states slowly approaches the frac- tion of channels in all bound states. Generally, the analytical model slightly underestimates (about 9% of the difference between analytical model and numerical results) the fraction of AMPAR channels in states that are bound and desensi- tized. This may come from the assumption in the analytical model about directions of transitions away from open states, which go from O states to C states (rather than C ) i i−1 i and from D to C states (rather than C ). Thus, in the i i−1 i first order of sub-conductance states (see Fig. 2), some transitions from the open state go back to the unbound Fig. 7 Glutamate concentration obtained from Monte Carlo simulation ◦ −ωt state (rather than the first bound closed state as assumed in in 25 C and assumed fitting curve proportional to e 386 J Comput Neurosci (2018) 44:379–391 Fig. 8 Influence of different Q diffusion coefficients on synaptic conductance Fig. 9 Comparing ratio of peak amplitudes and peak times 35 Cand 25 C for different order approximations of higher sub-conductance states and quicker rise time and decay of synaptic conductance. However, increasing Q causes also a smaller AMPAR diff peak conductance than is observed experimentally. This comparison to the first order approximation), with only supports the previous conclusion of Postlethwaite et al. a minor influence on the ratio of times to peak (Fig. 9). (2007) for the predominant role of postsynaptic kinetics Furthermore, it was found that it is possible to achieve in mediating temperature effects on synapses. In turn, this the same dynamics of AMPAR synaptic conductance result may be important in the context of possible medical (qualitatively and quantitatively) by using the 3rd order applications. Namely, the creation of a drug capable of approximation (but not the 1st or 2nd order) and changing altering receptor kinetics may lead to successful prevention the fraction of the peak conductance at the 3rd order of adverse temperature effects on synapse dynamics. open state from 0.7 to 0.9. Therefore, we support one of The significance of including higher subconductance the results of Postlethwaite et al. (2007), which claims states in the analytical model was also investigated. It that higher temperature drives AMPARs to higher sub- turned out, that including higher and higher states of conductance states (rather than only increasing unitary subconductance leads to saturating increase of ratio (in subconductances) cause the increase in conductance peak ◦ ◦ 35 C relatively to 25 C) of synaptic conductance peak amplitude). amplitudes (for about 12% in comparison to first order approximation), with minor influence on ratio of times of peak (Fig. 9). Furthermore, it was found that it is 4 Discussion possible to achieve same dynamics of AMPAR synaptic conductance (qualitatively and quantitatively) by using 3rd Normalization of fraction of channels Without any correc- approximation (without including 4-th order) and changing tion, the model would not be able to fit experimental data the fraction of the peak conductance at the 4-fold bound quantitatively (because the sum of all of the fractions of state for the 3rd state from 0.7 to 0.9. Therefore, we channels in different states does not equal 1). The reasons suppor the conclusion by Postlethwaite et al. (2007), are the assumptions (3) and (4): the fraction of channels in which claims that higher temperature leads AMPARs to state C for the (i+1)-th order mesh equals 1. To resolve higher conducting states (thus increasing conductance peak this problem, we introduce a normalization constant, which amplitude). scales the analytical model output to agree with the results The significance of including higher sub-conductance of detailed Monte Carlo simulation. states in the analytical model was also investigated. First, we calculated the sum of all bound channels in Including higher states of sub-conductance leads to a the analytical model. The results are presented in Fig. 10. ◦ ◦ saturating increase (in 35 C relative to 25 C) of synaptic The sum of all bound channels firstly rapidly increases conductance peak amplitudes (an increase of 12%, in and then after a time of approximately 0.5 ms gradually decreases over time. This is due to resensitization, which For example, a drug for slowing down the kinetics in a state of the drives AMPARs to the C0 state (as assumed in Scheme 2). brain hyperthermia - or a drug for speeding up kinetics in case of hypothermia. Additionally as the concentration of the glutamate in the J Comput Neurosci (2018) 44:379–391 387 mEPSC peak amplitude was not successfully reproduced by an independent binding model for the kinetic scheme they proposed. Although it is not possible to research this relationship using our differential-equation based model which describes only the averaged dynamics of the system, we suggest that the kinetic model proposed here (kinetic scheme in Fig. 2) may reproduce some aspects of the biological variability observed experimentally. To verify this we have run three sets of Monte Carlo simulations (with assumptions discussed in the Methods section) - two for kinetic scheme proposed by Postlethwaite et al. (2007) (one with cooperative and one with independent binding) and one set with the kinetic scheme proposed in this paper (with independent binding). The rate constants of the models were fitted to reproduce average behavior for the cooperative binding of Postlethwaite et al. (2007). We found that for the Fig. 10 Fraction of bound AMPAR channels for different orders of kinetic Scheme 2 Scheme 2 (Fig. 2) with independent binding skewness = 0.62 and CV = 0.27 are close to the values resulting from kinetic scheme by Postlethwaite et al. (2007) with cleft is low, the binding rate is substantially lower than cooperative binding (skewness = 0.53, CV = 0.31) unbinding rate, so no new channels get bound. unlikely to kinetic scheme by Postlethwaite et al. (2007) Second, we investigated what is the absolute number of with independent binding(skewness = 0.25, CV = 0.19). the bound channels in the Monte Carlo simulation over time. The absolute number of unbound channels (easily convertible for the absolute number of bound channels as Temperature dependence of AMPA receptor conductance they sum up to one) is presented in Fig. 4. Finally, we Our simulations show that, for a single synapse, increased compared the absolute number of bound channels (in Monte temperature causes larger peak amplitude of AMPAR con- Carlo simulation) at the peak of the synaptic conductance ductance, which is also achieved quicker (the conductance (0.5ms after stimulation) with the sum of all bound channels peaks faster in time). However, these results are not easily (in the analytical model) at the respective peak (0.3ms after interpretable in a more general context. One of the examples stimulation). A fraction of 0.13 of the channels is bound at is predicting an influence of temperature-modified synap- the peak in the Monte Carlo simulation, which corresponds tic conductance on temporal summation of signals across to ≈ 1 according to the analytical model. Therefore, to morphology of a neuron. As a single synaptic conductance quantitatively reproduce Monte Carlo results to a good has a larger peak amplitude in higher temperature, a smaller level of approximation, we only have to multiply all of the number of EPSPs should elicit an action potential at the fractions of channels in the analytical model by 0.13. The higher temperature. However, rise and decay time constants normalization factor is weakly dependent on the order of of the AMPAR conductance function are quicker at the subconductance assumed - as the 1st order contributes by higher temperature, so the effective summation window of far the most to the sum of all bound channels. the different signals should be shorter and therefore they should sum up less efficiently. These opposing features of AMPAR synaptic conductance hinder a simple description Glutamate binding model An analytical model showed of compounded temperature effects on multi-synaptic net- that, with the assumptions we made, it is possible to works: it is difficult to strictly predict (because we do not reproduce experimental and averaged numerical results of know the relative importance of amplitude and time con- the Monte Carlo simulation using independent glutamate stants of synaptic conductance curve) how temperature may binding. However, considering the complexity of biological influence temporal summation of the synaptic signal. In systems and the role of variability in their behavior, we the future, more detailed study on this topic may allow do not argue that independent binding is a biological us to investigate temperature effects from the level of sin- reality, but rather a reliable approximation for simulating gle synapses to that of large neural networks, which may aspects of AMPAR dynamics and temperature response. In help in better understanding of complex and paradoxical Postlethwaite et al. (2007), using Monte Carlo simulation, field interactions in brain imposed by temperature (Ander- sen and Moser 1995). Additionally, detailed investigation it was shown that variability of both rise and decay time constants of AMPAR conductance as a function of the of the influence of temperature on field potentials may be 388 J Comput Neurosci (2018) 44:379–391 important in the context of temperature-sensitive epilepsy AMPA type synapse, capable of simulating temperature (Traub and Wong 1982). effects, has few potential applications. First, our analytical model was validated with experimental data and it may be implemented with high accuracy and efficiency in Uncoupling assumption accuracy To further test the accu- large neural network simulations, which (when combined racy of our uncoupling assumptions of differential equa- with previous studies about temperature-dependence of tions in the analytical model, we carried out additional the conductance of voltage-gated ion channels) may open Monte Carlo simulation of synaptic transmission for the new possibilities of researching temperature influences on kinetic scheme proposed here, by comparing the dynam- neural dynamics in computational neuroscience. Second, ics of Scheme 2 without uncoupling, to the dynamics with due to the generality of this model, it is weakly dependent uncoupling (Scheme 3). This comparison illustrated that the on the kinetic scheme or phenomenological method of uncoupling assumption is fulfilled with different accuracy synaptic conductance modeling we have chosen. The for each order (Fig. 11). model provides simple theoretical linking between some The assumption is best fulfilled for the first order of models created in one temperature to any other (in kinetic scheme (see Fig. 2) and the error associated with this a reasonable physiological range), by using the model order of sub-conductance is < 0.5%. For the second and developed here as a ‘linking bridge’, without performing fourth orders, the error is around 15%. The worst accuracy additional experiments (however its accuracy has to be is for the third order, with an error level of roughly 40% at further carefully tested). the peak of synaptic transmission. However, from (1), we This theoretical linking could be done as follows: may see that the accuracy of the uncoupling assumption is dependent also on a glutamate concentration: the higher the 1. To some synaptic conductance curve, we fit (with glutamate concentration is at PSDs, the better the accuracy free parameters being rate constants and glutamate of the assumption. Therefore these error percentages should concentration constants A, ω (with the constraints for be divided by a factor of about 1.8, as the glutamate their values being in physiological range) of the model concentration function is not able to correctly capture the developed here. dynamics of glutamate concentration in the synaptic cleft 2. In fitted model, to determine dynamics at a different (see Fig. 7). temperature, we multiply all of the AMPAR kinetic rate constants by an appropriate temperature dependent New modeling method of temperature eﬀects on AMPA factor (Q = 2.4 for our simulations) and hence create receptor The creation of an analytical model for the a new synaptic conductance curve. Fig. 11 Uncoupling assumption validation for different orders of kinetic scheme. ‘Left’ and ‘Right’ denote for numerical values of left and right side of equation (1) respectively J Comput Neurosci (2018) 44:379–391 389 Fig. 12 Scheme 4 3. Finally, we fit desired phenomenological model to the not known apriori (as described earlier in the Introduction). synaptic conductance curve we achieved in the previous Thus, we will find the temperature coefficients by fitting the step (see also Scheme 4 on Fig. 12). rate constants of the simple models from Schemes 5, 6 and 7 to the analytical model (developed here) - this process is This approach shall be insensitive to possible fitting described in more detail in the previous paragraph. The full parameter degeneracy, due to the use of the common solutions of the coupled linear ODEs describing Schemes multiplication Q = 2.4 factor. 5,6 and 7 are presented in Appendix B. The model may be efficiently implemented in NEURON For kinetic models from Scheme 5 and 6 solutions (Hines and Carnevale 1997) using NMODL (Hines and may be discussed together as r and r rate constants 2 3 Carnevale 2000), making it ready for easy implementation are indistinguishable from each other (so we can reduce in neural network simulations. them to one constant). Models from Scheme 5 and 6 were unable to capture the dynamics of the AMPA receptor Comparison of the performance of the model with the in different temperatures, as relatively large fitting error simple AMPA kinetic models Set of the simple AMPA could hinder the subtle temperature dependence of the receptor kinetic models have been proposed previously by amplitude and time courses of synaptic conductance. For Destexhe et al. (1994b). Here, we will analyze three of kinetic model from Scheme 7, it was possible to obtain the AMPA kinetic models proposed, which are shown in accurate fit of the data in 25 C. However, for the best fit Schemes 5, 6 and 7 (Fig. 13). Employing the assumptions in 35 C the model was underfitting the amplitude of the (1),(2),(5) and (6) (see Methods) we may solve the set synaptic conductance, which is crucial in the investigation of differential equations describing the kinetics of the of temperature dependence. The solution of the kinetic AMPA receptor for a given kinetic model and include model from Scheme 7 needs additional constraints on the influence of temperature by multiplication of the rate parameters during the fitting process (as there are square constants by Q coefficients (yet not assuming that all roots in the solution). Moreover, the Q parameters used in of the coefficients are equal to each other). However, this model do not have a physical meaning, as (for the best- temperature coefficients for the arbitrary kinetic scheme are fit) some of them were found to be lower than 1 (so they should decrease the speed of the conformational changes of AMPARs). Therefore, without the more careful investigation of the simpler kinetic models, we propose that the analytical model developed here is a more convenient and accurate way to include temperature effects on AMPA receptor conductance. 5 Conclusion In the present study, an analytical model of an AMPA- type synapse including temperature effects was created. The model was developed on the bases of previous Marko- vian models describing the kinetics of the AMPA receptor and was simplified by uncoupling of the differential equa- tion system, and by kinetic scheme modifications motivated by Monte Carlo simulation of synaptic transmission. This method may be used to make simple models of synaptic conductance easily-scalable for any temperature and may provide simple theoretical linking of neurobiological mea- surements (involving AMPA-type synapse) conducted in Fig. 13 Kinetic schemes 5,6,7 (Destexhe et al. 1994b)usedinthe comparison of performance with the model developed in this paper different temperatures. Due to its accuracy (in comparison 390 J Comput Neurosci (2018) 44:379–391 to experimental data) and efficiency, this model may be Additionally, we would like to thank Themistoklis Mavrogordatos and Dr. Paul Bartlett for the critical review of the language used in the used in large neural network simulations. This opens new paper. possibilities to research various temperature effects on neu- ral dynamics in large-scale multi-neuron experiments and Compliance with Ethical Standards simulations. It may provide a theoretical basis for better understanding of different neurological disorders associated Conﬂict of interests The authors declare that they have no conflict of with sub- and super- physiological temperatures. In con- interest. junction with some previously proposed models (Wazn ˙ y Open Access This article is distributed under the terms of the and Wojcik 2014), this approach may shed some light on Creative Commons Attribution 4.0 International License (http:// understanding paradoxical temperature influence in serious creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give neurological disorders like autism spectrum disorder (Helt appropriate credit to the original author(s) and the source, provide a et al. 2008), which will be a scope of our interest in the link to the Creative Commons license, and indicate if changes were forthcoming future. made. Acknowledgements We would like to thank Prof. Daniel Wojcik, Polish Children’s Fund and S.Staszic High School in Lublin for help in organizing an internship in the Nencki Institute of Polish Academy Appendix A of Sciences. We would like to thank Department of Neuroinformatics of Maria Curie-Sklodowska University in Lublin for lending the The solution of all order of differential equation system computational power. Special thanks also to Dr. Matthias Hennig and describing AMPA receptor kinetics: Mr. Piotr Kononowicz, who helped in the development of the work. tω (kc −P)t kc t t(P +S+ω) t(kc +ω) kc t e S e e e (k − G) e (ω − S) e (ω − S) −t(k +ω) 2 −t(k +P +2ω) c c y (t ) = 0.1Ae k − + + k + 0.4A e k − − + scaled o b o R PR PS PR(P − ω)(S − ω) R(P − ω)(S − ω)ω P(S − ω)ω(R + ω) t(P +2ω) e (ω − P) + k R(P − ω)(S − ω)(R + ω) t(P +S+ω) t(k +2ω) k t c c 2e (k − G)(−G + S − ω)(S − G) e (−G + S − ω)(ω − P )(S − G) e (k − G)(S − G) c c 3 −t(k +P +3ω) +0.35A e k − − + 2 2 PR(k − 3ω)(P − 2ω)(P − ω)(R + ω)(R + 2ω) R(P − 2ω)(P − ω)ω (R + ω)(R + 2ω) PRω (R + ω)(R + 2ω) t(P +3ω) t(k +ω) 2e 2e (k − G)(−G + S − ω) − − k R(k − 3ω)(R + ω)(R + 2ω) R(P − ω)ω (R + ω)(R + 2ω) t(kc +3ω) t(P +S+ω) e (S −G)(−G+S − ω)(ω−P)(3ω−G)(−G+k −3ω) e (k −G)(S −G)(−G+S −ω)(−G+k − 3ω) c c c 4 −t(k +P +4ω) +A e k − + 6R(P − 3ω)(P − 2ω)(P − ω)ω (R + ω)(R + 2ω)(R + 3ω) PR(k − 4ω)(P −3ω)(P −2ω)(P −ω)(R +ω)(R + 2ω)(R + 3ω) t(k +2ω) e (k − G)(−G + S − ω)(−G + k − 3ω) c c 2R(P − 2ω)ω (R + ω)(R + 2ω)(R + 3ω) t(kc +ω) t(P +4ω) kc t e (k − G)(S − G)(−G + k − 3ω) e e (k − G)(S − G)(−G + S − ω) c c c − − + k 3 3 2R(P − ω)ω (R + ω)(R + 2ω)(R + 3ω) R(k − 4ω)(R + ω)(R + 2ω)(R + 3ω) 6PRω (R + ω)(R + 2ω)(R + 3ω) where, S = k − ω, R =−k + k + k + k , P = where B is a constant. c c d o u k + k + k − ω, G = k + k + k Solution of the system of two coupled ODEs described d o u d o u by kinetic model of Destexhe et al. (1994b) (Scheme 6): Appendix B t(−(r +r )) t(r +r −ω) 2 3 2 3 Br e e − 1 y(t) = r + r − ω Solution of the system of two coupled ODEs described by 2 3 kinetic model of Destexhe et al. (1994b) (Scheme 5): Solution of the system of two coupled ODEs described −r t t(r −ω) 2 2 Br e e − 1 y(t) = by kinetic model of Destexhe et al. (1994b) (Scheme 7): r − ω 2 2 2 − t (r +r +r ) −4r r +r +r +r 2 2 t (r +r ) +2(r −r )r +r 2 3 4 2 4 2 3 4 2 t ((r +r ) +2(r −r )r +r ) 2 3 3 2 4 2 3 3 2 4 4 4 y(t) = Br e −1 + e r (r − ω)− r −1 + e 1 2 4 4 √ √ 2 2 2 2 t (r +r ) +2(r −r )r +r 0.5t(r +r +r + (r +r ) +2(r −r )r +r )−2ω) 2 3 3 2 4 4 2 3 4 2 3 3 2 4 4 2 2 + 1 + e − 2e (r + r ) + 2(r − r )r + r (r − ω) 2 3 3 2 4 4 4 2 2 t (r +r ) +(r −r )r +r 2 2 3 3 2 4 4 2 − −1 + e r (r − ω) 2 −4r r + (r + r + r ) r r − (r + r + r )ω + ω 3 4 2 4 2 3 4 2 4 2 3 4 J Comput Neurosci (2018) 44:379–391 391 References underlies temperature-dependent changes in synaptic strength at the rat calyx of Held. The Journal of Physiology, 579(1), 69–84. Robert, A., & Howe, J.R. (2003). How AMPA receptor desensitization Andersen, P., & Moser, E.I. (1995). Brain temperature and hippocam- depends on receptor occupancy. The Journal of Neuroscience, pal function. Hippocampus, 5(6), 491–498. 23(3), 847–858. Asztely, F., Erdemli, G., Kullmann, D. M. (1997). Extrasynaptic glu- Rosenmund, C., Stern-Bach, Y., Stevens, C. (1998). The tetrameric tamate spillover in the hippocampus: dependence on temperature structure of a glutamate receptor channel. Science, 280, 1596– and the role of active glutamate uptake. Neuron, 18, 281–293. Badjatia, N. (2009). Hyperthermia and fever control in brain injury. Sahara, Y., & Takahashi, T. (2001). Quantal components of the Critical care medicine, 37(7), S250–S257. excitatory postsynaptic currents at a rat central auditory synapse. Buzatu, S. (2009). The temperature-induced changes in membrane Journal of Physiology, 536, 189–197. potential. Biology Forum/Rivista di Biologia, 102(2), 199–217. Satzler, ¨ K., Sohl, ¨ L.F., Bollmann, J.H., Borst, J.G.G., Frotscher, M., Cais, O., Sedlacek, M., Horak, M., Dittert, I., Vyklicky, L. et al. (2002). Three-dimensional reconstruction of a calyx of (2008). Temperature dependence of NR1/NR2B NMDA receptor Held and its postsynaptic principal neuron in the medial nucleus channels. Neuroscience, 151(2), 428–438. of the trapezoid body. Journal of Neuroscience, 22(24), 10567– Destexhe, A., Mainen, Z.F., Sejnowski, T.J. (1994a). An efficient method for computing synaptic conductances based on a kinetic Schiff, S.J., & Somjen, G.G. (1985). The effects of temperature on model of receptor binding. Neural Computation, 6(1), 14–18. synaptic transmission in hippocampal tissue slices. Brain research, Destexhe, A., Mainen, Z.F., Sejnowski, T.J. (1994b). Synthesis 345(2), 279–284. of models for excitable membranes, synaptic transmission and Scimemi, A., & Beato, M. (2009). Determining the neurotransmitter neuromodulation using a common kinetic formalism. Journal of concentration profile at active synapses. Molecular neurobiology, Computational Neuroscience, 1(3), 195–230. 40(3), 289–306. Dietrich, W.D. (1992). The importance of brain temperature in cerebral Shelley, C., & Magleby, K.L. (2008). Linking exponential components injury. Journal of Neurotrauma, 9, S475–85. to kinetic states in Markov models for single-channel gating. Fuxe, K., Rivera, A., Jacobsen, K.X., Hoistad, M., Leo, G., Horvath, Journal of Eneral Physiology, 132(2), 295–312. T.L., et al. (2005). Dynamics of volume transmission in the brain. Shigemori, M., Abe, T., Aruga, T., Ogawa, T., Okudera, H., Ono, Focus on catecholamine and opioid peptide communication and J., et al. (2012). Guidelines for the Management of Severe Head the role of uncoupling protein. Journal of Neural Transmission, Injury, 2nd Edition Guidelines from the Guidelines Committee 112, 65–76. on the Management of Severe Head Injury, the Japan Society Helt, M., Kelley, E., Kinsbourne, M., Pandey, J., Boorstein, H., of Neurotraumatology. Neurologia Medico-Chirurgica, 52(1), 1– Herbert, M., Fein, D. (2008). Can children with autism recover? If so, how Neuropsychology Review, 18(4), 339–366. De Schutter, E., Roth, A., van Rossum, M.C. (2009). 6. Computational Hille, B. (2001). Ion channels of excitable membranes Vol. 507. modeling methods for neuroscientists, (pp. 139–160). Cambridge: Sinauer: Sunderland. The MIT Press. Hines, M., & Carnevale, N. (1997). The NEURON simulation Smith, T., Wang, L., Howe, J. (2000). Heterogeneous conductance environment. Neural Computation, 9, 1179–1209. levels of native AMPA receptors. Journal of Neuroscience, 20, Hines, M., & Carnevale, N. (2000). Expanding NEURON’s repertoire 2073–2085. of mechanisms with NMODL. Neural Computation, 12, 995– Sterratt, D.C. (2015). Q10: the effect of temperature on ion channel kinetics. In: Encyclopedia of Computational Neuroscience, pp. Hodgkin, A.L., & Huxley, A.F. (1952). A quantitative description of 2551–2552. membrane current and its application to conduction and excitation Stiles, J.R., Van Helden, D., Bartol, T.M., Salpeter, E.E., Salpeter, in nerve. The Journal of Physiology, 117(4), 500. M.M. (1996). Miniature endplate current rise times <100 s Kowacs, P.A., Marchioro, I.J., Silva Jr, E.B.D., Rocha, S.F., Simao, ˜ from improved dual recordings can be modeled with passive C.A., Meneses, M.S. (2005). Hot-water epilepsy, warm-water acetylcholine diffusion from a synaptic vesicle. Proceedings of the epilepsy, or bathing epilepsy? Report of three cases and National Academy of Sciences, 93, 5747–5752. considerations regarding an old theme. Arquivos de Neuro- Traub, R.D., & Wong, R.K. (1982). Cellular mechanism of neuronal psiquiatria, 63(2B), 399–401. synchronization in epilepsy. Science, 216(4547), 745–747. Mrozek, S., Vardon, F., Geeraerts, T. (2012). Brain temperature: Wazn ˙ y, M., & Wojcik, G.M. (2014). Shifting spatial attention- physiology and pathophysiology after brain injury. Anesthesiology numerical model of Posner experiment. Neurocomputing, 135, research and practice, 2012, 13. 139–144. Postlethwaite, M., Hennig, M.H., Steinert, J.R., Graham, B.P., Weight, F.F., & Erulkar, S.D. (1976). Synaptic transmission and effects Forsythe, I.D. (2007). Acceleration of AMPA receptor kinetics of temperature at the squid giant synapse. Nature, 261, 720–722.
Journal of Computational Neuroscience
– Springer Journals
Published: May 11, 2018