Postsynaptic Potential Energy as Determinant of Synaptic Plasticity

Metabolic energy can be used as a unifying principle to control neuronal activity. However, whether and how metabolic energy alone can determine the outcome of synaptic plasticity remains unclear. This study proposes a computational model of synaptic plasticity that is completely determined by energy. A simple quantitative relationship between synaptic plasticity and postsynaptic potential energy is established. Synaptic weight is directly proportional to the difference between the baseline potential energy and the suprathreshold potential energy and is constrained by the maximum energy supply. Results show that the energy constraint improves the performance of synaptic plasticity and avoids setting the hard boundary of synaptic weights. With the same set of model parameters, our model can reproduce several classical experiments in homo- and heterosynaptic plasticity. The proposed model can explain the interaction mechanism of Hebbian and homeostatic plasticity at the cellular level. Homeostatic synaptic plasticity at different time scales coexists. Homeostatic plasticity operating on a long time scale is caused by heterosynaptic plasticity and, on the same time scale as Hebbian synaptic plasticity, is caused by the constraint of energy supply.


INTRODUCTION
Although the brain accounts for only 2% of body mass, it consumes 20% of the resting metabolic energy produced by the whole body (Attwell and Laughlin, 2001;Harris et al., 2012). Within the brain, neurons utilize 75-80% of this energy, and the remainder is used by the neighboring glial cells. Housekeeping tasks use 25% of the total neuronal energy. Maintaining resting membrane potential (15%), firing action potentials (16%), and synaptic transmission (44%) compose the energetically most expensive processes (Harris et al., 2012;Howarth et al., 2012). Thus, the majority of energy used by neurons is locally consumed at the synapse. In addition to the energetic costs of neural computation and transmission, experimental evidence indicates that synaptic plasticity is metabolically demanding (Mery and Kawecki, 2005;Jaumann et al., 2013;Placais and Preat, 2013;Placais et al., 2017). The energy cost of synaptic plasticity is estimated based on the neurophysiological and proteomic data of rat brains depending on the level of protein phosphorylation; this cost constitutes a small fraction of the energy used for fast excitatory synaptic transmission, which is typically 4.0-11.2% (Karbowski, 2019). However, the quantitative relationship between the changes in synaptic weights (potentiation or depression) and energy consumption remains unclear.
Considering the consistency of corresponding experiments, a large number of synaptic plasticity models have been established. These models are mainly biophysical models based on calcium hypothesis (Shouval et al., 2002;Graupner and Brunel, 2012) and phenomenological models based on pre-and post-synaptic spikes or voltage (Bienenstock et al., 1982;Pfister and Gerstner, 2006;Clopath et al., 2010). Although these models are successful in experimental reproduction, they ignore the role of metabolic energy. Growing evidence suggests that metabolic energy may be a unifying principle governing neuronal activities (Laughlin, 2001;Niven and Laughlin, 2008;Hasenstaub et al., 2010;Yu and Yu, 2017), thereby naturally leading people to focus on the relationship between metabolic energy and synaptic plasticity in recent years. Sacramento et al. (2015) showed that unbalanced synaptic plasticity rules can lead to sparse connectivity and energy-efficient computation. Li and van Rossum (2020) assumed that the metabolic energy for every modification of a synaptic weight is proportional to the amount of change, regardless of whether this is positive or negative. They proposed a synaptic caching algorithm based on this assumption. The proposed algorithm can enhance energy efficiency manifold by precisely balancing labile forms of synaptic plasticity with many stable forms. However, energy is expressed by synaptic weights in these studies. Whether synaptic plasticity can be fully quantified by energy remains unclear.
Potential energy is stored in transmembrane ion gradients. When postsynaptic neurons are stimulated by external stimuli (such as synaptic input), the changes in gating state, channel conductance, and current are driven by the energy stored in the Na + and K + gradients, and no adenosine triphosphate (ATP) is consumed in this process. These gradients and stored potential energy are partially depleted and must be actively restored. The change in postsynaptic potential energy indirectly reflects the consumption or supply of metabolic energy because the active recovery of potential energy needs ATP. In this study, we express the postsynaptic potential energy as the integral of the product of postsynaptic membrane potential and the postsynaptic membrane current density on stimulation time. The potential energy with membrane potential lower than a certain threshold is called subthreshold potential energy. The part with membrane potential greater than the threshold is called suprathreshold potential energy. The baseline potential energy is the result of downscaling the amplitude of the subthreshold potential energy. The synaptic weights are expressed by a simple linear relationship between the subthreshold potential energy and the suprathreshold potential energy and are constrained by the energy supply. The simulation results show that the model can reproduce a series of classic synaptic plasticity experiments, which indicate that our model is feasible.

Construction of Synaptic Plasticity Model
Our model uses postsynaptic potential energy to express the change in synaptic weights. Postsynaptic potential energy P is the integral of the product of postsynaptic membrane potential V m and postsynaptic membrane current density I m to stimulation time t, that is, P = V m I m dt. To explain more clearly how the model works, we divide the operation process of the model into four stages at each time step ( Figure 1A). In stage (1), if no stimulation occurs, the postsynaptic neuron is in a resting state. At this time, postsynaptic potential energy P is the resting state potential energy P rest . Taking the resting state potential energy as the reference point of potential energy, let P = P rest = 0. In stage (2), stimulation causes potential energy P to deviate from the resting state potential energy. Potential energy P after stimulation is separated into two parts. The first part is called subthreshold potential energy P sub , and its membrane potential V m is less than the threshold potential V th . The second part is called suprathreshold potential energy P sup , and its membrane potential V m is greater than V th . Thus, P = P sub + P sup . In stage (3), the subthreshold potential energy is multiplied by a constant A r between zero and one, which is called the baseline coefficient, and the baseline potential energy P bas is obtained. Thus, P bas = A r P sub . In stage (4), we assume that the change in synaptic weights is proportional to the difference between the baseline potential energy and the suprathreshold potential energy and then test the rationality of this hypothesis by comparing it with a series of synaptic plasticity experimental results. Therefore, the change in synaptic weight W is expressed as follows: where A is a positive constant, which represents the linear transformation between potential energy and synaptic weight and is called the amplitude coefficient (equivalent to the learning rate). The synaptic weights increase when the baseline potential energy is greater than the suprathreshold potential energy ( Figure 1B). The synaptic weights decrease when the baseline potential energy is less than the suprathreshold potential energy ( Figure 1C). The synaptic weights do not change when they are the same. We assume that the amplitude of postsynaptic potential energy P cannot be greater than the maximum energy supply S because the change in potential energy is constrained by energy supply. Unless otherwise specified, the energy supply in this study represents the maximum energy supply that can be provided. We call the maximum potential energy P max where its amplitude is the same as energy supply S, but the sign is consistent with potential energy P. Then, P max = S sign (P), where sign is the sign function, and |P| ≤ S and |P| ≤ |P max |. Given that the dynamic characteristics of energy supply are unclear, we propose a simple formula for calculating the energy supply with time where t is the time of stimulation, and τ is the time constant of energy supply which is usually much larger than the membrane time constant (Attwell and Laughlin, 2001). R is equivalent to the total energy supply of synapses per unit time, which is a constant and is called the rate of energy supply. S 0 is the minimum energy supply to maintain the normal function of neurons, which is a constant greater than zero and equivalent to the energy supply in the resting state. The formula represents the maximum energy that can be provided at time t by multiplying Rt of the energy supply linearly increasing with time and a damping factor S damp = e −t/τ , which decreases exponentially with time. S = S 0 at t = 0 FIGURE 1 | Model schematic. The height of the bars represents the amplitude of different potential energy (P, P sub , P sup , etc.). P is the change value relative to P rest . Because we assume that P rest is 0, we call P the potential energy. (A) Mathematical expression of synaptic plasticity Model; (B) calculation of synaptic weights when the baseline potential energy is greater than the suprathreshold potential energy. In the resting state, the potential energy of the postsynaptic membrane is assumed to be zero. The potential energy P of the postsynaptic membrane is lower than that of the resting state. Thus, P < P rest is negative. Subthreshold potential energy P sub and suprathreshold potential energy P sup are negative. Baseline coefficient A r is greater than zero, so baseline potential energy P bas has the same sign as subthreshold potential energy P sub and is also negative. Here, P bas > P sup , and the synaptic weight increases; (C) calculation of synaptic weights when P bas < P sup . Similar to that in (B), the synaptic weight decreases at this time because P bas < P sup ; (D) calculation of synaptic weight when the amplitude of potential energy P exceeds energy supply S at every time step. The amplitude of P is adjusted to the maximum potential energy P max , which is the same as S. Assuming that the membrane potential is greater than the threshold potential, the amplitude of P sup is reduced. The amplitude of P sub is reduced if the membrane potential is less than the threshold potential. Thus, P sub + P sup = P = P max .
or t → ∞, the energy supply is minimum, and the energy supply is maximum at t = τ . Some researches showed that the relationship between Na/K pump current and intracellular sodium ion concentration and ATP concentration can be expressed by the Hill equation (Figure 4 in Despa and Bers, 2003;Figures 2, 10 in Glitsch, 2001). Therefore, many scholars take the intracellular sodium ion concentration as the measurement of energy consumption (Hasenstaub et al., 2010). The concentration of sodium ions in stimulated neurons increases from the concentration in the resting state to a peak (at the end of stimulation) and then returns to the concentration in the resting state due to the action of the ion pump. The changes in postsynaptic membrane potential and ionic current caused by spikes in neurons also have the characteristics of first rising and then returning to the original position with time, which is usually described in the exponential form of t e −t/τ (Gerstner, 1995;Bohte et al., 2002). Given this, we also use this exponential form to describe the kinetics of intracellular sodium ions as an expression of energy consumption or energy supply. A similar exponential expression has also been successfully applied in the constraints of metabolic energy on the synaptic connection (Yuan et al., 2018).
Potential energy P is adjusted to P max if its amplitude exceeds energy supply S ( Figure 1D) so that its amplitude is equal to energy supply S. This adjustment results in a decrease in the amplitude of suprathreshold potential energy P sup or subthreshold potential energy P sub so that their sum is equal to P max . The subthreshold potential energy is adjusted, and the suprathreshold potential energy remains unchanged when the membrane potential at time t is smaller than the threshold potential. The suprathreshold potential energy is adjusted, and the subthreshold potential energy remains unchanged when the membrane potential at time t is greater than the threshold potential. The adjustment of subthreshold potential energy causes the amplitude of baseline potential energy P bas to change in the same proportion (as shown in the next section). The specific implementations of Equations (1) and (2)

Determination of Model Parameters
In accordance with Equations (1) and (2), the model includes six parameters, namely, amplitude coefficient A, baseline coefficient A r , threshold potential V th , energy supply rate R, minimum energy supply S 0 , and time constant of energy supply τ . For the frequency-dependent pairing protocol used by Sjöström et al. (2001), we chose the final parameter by trial and error. We simulated and tested 27 sets of parameters R = 100, 150,200, τ = 1, 2, 3 and S 0 = 10, 20, and 30, respectively. By comparing the simulation results with the results of synaptic plasticity experiment, the parameter set which is most consistent with the experimental results was selected: A = 0.02, A r = 0.2, V th = −60 mV, R = 175 fJ/(µm 2 s), τ = 2 s, S 0 = 25 fJ/µm 2 . In the next section, we introduce several classical experimental protocols of synaptic plasticity and illustrate how our model reproduces these different experimental results with the same set of parameters through potential energy and energy supply.

Reproduction of the Experimental Results of Homosynaptic Plasticity
Nineteen distal and proximal compartments (magenta, Figure 2A) were simulated in the basal dendrites of the L5 pyramidal neuron model. We followed two different experimental protocols on homosynaptic plasticity to compare our model with the experimental data. The first protocol was the classical spike-timing-dependent plasticity (Markram et al., 1997;Bi and Poo, 1998;Sjöström et al., 2001). Each distal and proximal compartments were connected to one synapse. Postsynaptic spikes were induced by the injection of 1 nA and 3 ms current pulses into the soma of postsynaptic neurons. The initial synaptic weights were set to 0.5. For the study of spike frequency dependence, pairs of pre-post ( Figure 2B) or post-pre ( Figure 2C) spikes separated by 10 ms were repeated five times at different frequencies of 5 Hz up to 50 Hz with steps of 5 Hz and for 0.1 Hz. For the study of spike timing dependence ( Figure 3A), pairs of pre-post or post-pre spikes at 20 Hz were repeated five times for different time intervals t (1, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, and 20 ms). The computational results (weights and potential energy) repeated five times were multiplied by a scaling factor of 12 (60/5) to mimic 60 pairs of presynaptic and postsynaptic spikes in the experimental protocols of Sjöström et al. (2001). The study of Sjöström et al. (2001) focused on the weight change as a function of the frequency for a fixed t in this pairing protocol. The second protocol was synaptic afferent where only presynaptic spikes were induced (Bliss and Lømo, 1973;Dudek and Bear, 1992). To test our model on a consistent set of data, we took the measurements of Dudek and Bear (1992) in this study because sufficient quantitative information can be found in their study. Each distal and proximal compartments were connected to one synapse. The activation of the synapse connected to each compartment consisted of 20 pulses delivered by a Poisson process at input frequencies ranging from 1 to 50 Hz (1, 3, 5, 10, 20, 30, 40, and 50 Hz) ( Figure 3B).
For the frequency-dependent pairing protocol without energy supply constraints, the amplitude of postsynaptic potential energy decreased with the increase in spike frequency (blue solid lines; top; Figures 2B,C); this finding is consistent with the calculation results of the relationship between metabolic energy and frequency in neurons (Yi et al., 2016). With the decrease in spike frequency, the time of stimulation increased gradually. Increasing cases were found, where the amplitude of postsynaptic potential energy without energy constraint (named as unconstrained energy) exceeded that of the maximum potential energy P max (blue shaded area under the red dotted line; top; Figures 2B,C). The unconstrained energy with amplitude greater than that of P max was adjusted to the same as P max to obtain the postsynaptic potential energy with energy supply constraints (named as constrained energy) (black solid lines; top; Figures 2B,C). The amplitudes of unconstrained energy in all postsynaptic compartments were greater than the amplitudes of P max when the frequency was <10 Hz. At this time, the amplitudes adjusted were the largest, resulting in the overlap between the constrained energy and P max . The adjustment of potential energy reduced the amplitude For the time-dependent pairing protocol (Figure 3A), the stimulation time of different pairing time intervals is the same because the spike frequency was fixed at 20 Hz. Thus, the maximum potential energy does not change with the time interval (red dashed lines; top; Figure 3A). Similar to the analysis in the previous section, the unconstrained energy with amplitude greater than P max was adjusted the same as P max . The amplitude of constrained energy (black solid lines; top; Figure 3A) was smaller than that of the unconstrained energy (blue solid lines; top; Figure 3A). These adjustments led to the corresponding changes in the baseline potential energy and suprathreshold potential energy (middle; Figure 3A) and made the synaptic plasticity constrained by energy supply more consistent with the experimental results than that without energy supply constraint (bottom; Figure 3A).
In the synaptic afferent protocol (Figure 3B), increasing cases were found, where the amplitude of postsynaptic potential energy without energy constraint exceeded that of P max , when the input frequency was <5 Hz. The adjustment of this unconstrained energy to the maximum potential energy led to the gradual approaching and overlapping of the constrained energy and the maximum potential energy (top; Figure 3B). If the amplitude of unconstrained energy is less than that of maximum potential energy, the constrained and unconstrained variables (i.e., potential energy, baseline and suprathreshold potential energy, and weights) should overlap because it is unnecessary to adjust the potential energy. However, the constrained and unconstrained variables did not overlap when the frequency was >5 Hz although the amplitude of unconstrained energy (blue solid lines; top; Figure 3B) was less than that of the maximum potential energy (red dashed lines; top; Figure 3B). This phenomenon was because all the calculation results in this study were the values at the end of stimulation. However, the adjustment of potential energy was conducted at every moment from the beginning to the end of stimulation (section Methods). These results indicated that before the end of stimulation, the unconstrained energy was adjusted several times because the amplitude exceeded that of the maximum potential energy. During the stimulation, the main adjustment was reflected in the suprathreshold potential energy if the input frequency was >5 Hz. However, the baseline potential energy was mainly adjusted if the input frequency was <5 Hz (middle, Figure 3B). The computational results showed that our model can quantitatively reproduce the results of the synaptic afferent experiment (bottom, Figure 3B).

Reproduction of Mexican Cap-Like Heterosynaptic Long-Term Depression
Heterosynaptic plasticity can be induced at synapses that are inactive during the induction of homosynaptic plasticity (Chistiakova et al., 2015;Zenke et al., 2017). High-frequency afferent tetanization induces a Mexican hat-like profile of response amplitude changes: Homosynaptic long-term potentiation (LTP) at stimulated inputs is surrounded by heterosynaptic long-term depression (LTD) (White et al., 1990;Royer and Paré, 2003). Each compartment on the thin basal branch of the L5 pyramidal neuron model is connected to two synapses. The initial synaptic weights were set to 0.5. To reproduce the heterosynaptic LTD, we used a similar protocol to that of Royer and Paré (2003). Homosynaptic LTP was induced with high-frequency stimuli (HFS) at the synapses connected to each distal or proximal compartment. HFS consisted of four series of 10 trains separated by 0.3 s, where each train consisted of 10 shocks (Poisson process) at 100 Hz. Each basal dendrite was divided into seven sites to compare with the experimental results. Site 0 was a compartment corresponding to homosynaptic LTP (magenta dots; Figure 4). Other compartments connected by heterosynapses were divided into six sites in accordance with the distance from site 0 (cyan dots; Figure 4). The value in each site (such as potential energy, weight, etc.) was the average for the values of all compartments in the site.
Homosynaptic LTP was induced in the synapses at the stimulation site, and heterosynaptic LTD (bottom, Figures 4B,D) occurred in the non-activated synapses of the same branch when HFS was performed at the proximal ( Figure 4A) or distal ( Figure 4C) in a thin branch. The statistical results for proximal and distal stimulation in all basal branches ( Figure 4E) showed that the heterosynaptic LTD decreased with the increase in distance from the stimulus site whereas homosynaptic LTP was induced. This Mexican hat-like heterosynaptic plasticity was in good agreement with the experimental results (bottom, Figure 4F). For heterosynaptic sites (sites 1-6), the amplitudes of unconstrained energy were always lower than that of the maximum potential energy. Therefore, the maximum potential energy did not have a constraining effect, which results in the overlap of constrained and unconstrained potential energy (top; Figures 4B,D,F). The homosynaptic LTP increases unlimitedly if it was not constrained by energy. This condition was because the unconstrained energy with amplitude greater than the maximum potential energy was not adjusted, and the difference between the baseline and suprathreshold potential energy was extremely large (middle; Figures 4B,D,F). The energy constraint made the postsynaptic potential energy equal to the maximum potential energy, thereby reducing the difference between the baseline and suprathreshold potential energy and controlling the homosynaptic LTP in a biologically reasonable range (bottom; Figures 4B,D,F). Although the sites of hetero LTD were more than that of homosynaptic LTP, the amplitude of homosynaptic LTP was larger than that of heterosynaptic LTD by comparing homosynaptic LTP and heterosynaptic LTD. As shown in the top panel of Figures 4B,D,F, the difference in the maximum potential energy was due to the different signs of the corresponding potential energy, and their energy supply was the same.

DISCUSSION
We presented a computational model of synaptic plasticity completely determined by energy and established a simple quantitative relationship between synaptic plasticity and postsynaptic potential energy. The synaptic weight is directly proportional to the difference between the baseline and suprathreshold potential energy and is constrained by the maximum energy supply. Considering that the dynamic characteristics of energy supply are unclear, we proposed a simple dynamic equation of energy supply and provided the upper limit of the amplitude of postsynaptic potential energy.  A,B), but for distal stimulation; (E,F) all stimulated sites (E) and the corresponding computational results (F). The lines (solid and dashed) and the shaded regions are the mean and SD, respectively, and overall proximal and distal compartments (magenta) shown in (E). Dot line with errors (bottom, F) represents the experimental data from Royer and Paré (2003). For clarity, the SDs of the potential energy with constraints (black lines, top), baseline and suprathreshold energy (middle), and weights without constraints (blue lines, bottom) are not shown.
The constraint of energy supply improves the performance of synaptic plasticity and avoids setting the hard boundary of synaptic weights. In the classical frequency-dependent pairing protocol, six parameters of the model were determined by trial and error. With such a set of parameters, our model reproduced several experimental results of homosynaptic plasticity and the Mexican hat-like heterosynaptic LTD, showing that our model can unify the homo-and heterosynaptic plasticity.

Quantitative Relationship Between Synaptic Plasticity and Metabolic Energy
The absolute value of potential energy can be regarded as the metabolic energy consumed because it restores postsynaptic potential energy to the resting state. Our model assumes that the synaptic weight is proportional to the difference between the baseline and suprathreshold potential energy and is constrained by the maximum energy supply. The model can reproduce a series of experimental results of homo-and heterosynaptic plasticity, which shows that our hypothesis is feasible. Can metabolic energy replace potential energy to express this linear relationship? If the potential energy, baseline, and suprathreshold potential energy are all negative, the absolute value of the baseline potential energy is called the baseline metabolic energy. The absolute value of the suprathreshold potential energy is called the suprathreshold metabolic energy. If the baseline or suprathreshold potential energy is a constant that does not change with time, then Li and van Rossum's (2020) hypothesis that metabolic energy is directly proportional to the change in synaptic weight is correct. The linear relationship between synaptic weights and metabolic energy is only valid in some cases.

Interaction Mechanism Between Hebbian and Homeostatic Synaptic Plasticity
At present, the biological significance of Hebbian synaptic plasticity (positive feedback) and homeostatic synaptic plasticity (negative feedback) remains controversial. Specifically, how these opposing forms of plasticity that share common downstream mechanisms work in the same networks, neurons, and synapses remain unclear (Turrigiano et al., 1998;Feldman, 2002;Turrigiano and Nelson, 2004;Swanwick et al., 2006;Rannals and Kapur, 2011). In recent years, these conditions have been discussed extensively by leading experts in the field (Vitureira and Goda, 2013;Fox and Stryker, 2017;Keck et al., 2017;Yee et al., 2017). One view is that homeostatic plasticity operates on a long time scale and does not interfere with synaptic changes induced by Hebbian plasticity (Turrigiano, 2012;Tononi and Cirelli, 2014;Hengen et al., 2016). Another view is that Hebbian and homeostatic synaptic mechanisms may be parallel; thus, they can interfere with each other in the same synaptic subset (Desai et al., 2002;Kim and Tsien, 2008;Keck et al., 2011;Vlachos et al., 2013;Frank, 2014;Li et al., 2014). Based on observing the existence of fast and input-specific homeostatic mechanisms, a signal pathway-based model was proposed to adjust the balance between Hebbian and homeostatic synaptic plasticity. The model was used to explain the interaction mechanism between Hebbian and homeostatic plasticity on the same time scale (Galanis and Vlachos, 2020).
Our model can integrate the two different viewpoints and give a unified explanation. The homeostatic synaptic plasticity at different time scales coexists. First, we believe that the homeostatic plasticity operating on a long time scale is caused by heterosynaptic plasticity. Our simulation showed that the amplitude of heterosynaptic LTD is extremely small, especially under high-frequency tetanic stimulation (bottom; Figures 4B,D,F). This condition indicates that the changes in heterosynaptic plasticity caused by normal neural activities are extremely small or difficult to confirm under the Hebbian time scale. Experiencing a longer stimulation than the Hebbian time scale is necessary before these changes can be evident. This heterosynaptic plasticity accumulated over a long period becomes homeostatic plasticity on a long time scale. Second, homeostatic synaptic plasticity on the same time scale as Hebbian synaptic plasticity is caused by the constraint of energy supply. The synaptic strength does not increase continuously under high-frequency stimulation nor does it decrease unlimitedly under low-frequency stimulation due to the constrained energy (bottom; Figures 2, 3B). On the contrary, the amplitude of synaptic enhancement or inhibition is reduced to match the energy supply because the energy supply gradually decreases after the time of stimulation is greater than its time constant (Equation 2). The final result is a homeostatic synaptic plasticity parallel to the Hebbian time scale. We propose a unified mechanism for the interaction between Hebbian and homeostatic synaptic plasticity based on the above analysis. The homeostatic homo-and heterosynaptic plasticity coexist with homo-and heterosynaptic plasticity. The time scale of homeostatic homosynaptic plasticity is the same as that of homosynaptic plasticity (i.e., Hebbian synaptic plasticity), which is rapid and input-specific and is caused by the limitation of energy supply. The homeostatic heterosynaptic plasticity has a long time scale, which is caused by the long-term accumulation of heterosynaptic plasticity. Although we do not fully understand the molecular mechanism of heterosynaptic plasticity and energy supply and the actual energy supply dynamics, we believe that the analysis of this mechanism from the cellular level is still valuable.

Limitations of Our Approach
First, our synaptic plasticity model can reproduce a series of classical synaptic plasticity experiments using a detailed biophysical model of a single pyramidal neuron. Although the feasibility of the model can be confirmed, further examining the consistency between the model and the experimental results under many stimulation protocols is necessary. Second, in our model parameters, the threshold potential V m of −60 mV corresponds to the leakage potential E L (section Methods) in the biophysical model of neurons. The baseline coefficient A r of 0.2 indicates that the baseline potential energy is one-fifth of the subthreshold potential energy. Whether the two parameters are universal for different neuron models and their biophysical significance remains unclear. Third, the proposed dynamic equation of energy supply (Equation 2) is not supported by experimental data. Can the equation be used as a theoretical prediction to guide future experiments? Can different dynamic equations supported by experiments achieve the same effect of energy constraint in our model? These questions are worthy of further exploration. Finally, the interaction between Hebbian and homeostatic plasticity for large-scale neural networks has an important influence on the learning and memory ability of neural networks. We did not study the validity and scalability of the model in the neural network environment, especially in largescale neural networks, which is an important direction of our future work.

Model of Neuron and Synapse
All simulations in this study were conducted on Brian 2 neuron simulator in Python (Goodman and Brette, 2009). We used the model and parameters of the biophysical neurons and synapses developed by Bono and Clopath (2017) on Brian 2. The implementation of the pyramidal neuron and synaptic model completely adopted the codes posted by Bono and Clapath on ModelDB (https://senselab.med.yale.edu/modeldb/). Given that L5 and L2/3 pyramidal neuron models are found to have similar results, we only conducted simulation studies on synaptic plasticity in L5 pyramidal neurons (Bono and Clopath, 2017). L5 pyramidal neurons are composed of a spherical soma, an axon, and many dendrite branches, with a total of 1,181 compartments. The leakage potential E L and resting potential of each compartment are −60 and −69 mV, respectively. Following codes of the neuron and synaptic model of Bono and Clopath (https://senselab.med.yale.edu/modeldb/), we can obtain the membrane potential V m and membrane current density I m of each compartment under any stimulation protocol. The unit of membrane potential is mV, and the unit of membrane current density is ampere/meter 2 , abbreviated as A/m 2 . We chose the unit of membrane current density as pA/µm 2 equivalent to A/m 2 because the geometric size of neurons in Brian 2 is usually expressed in µm.

Synaptic Plasticity
The variables in Brian 2 simulator are usually expressed and calculated in the form of differential equations. Thus, our synaptic plasticity model needs to be realized in the form of differential equations.
The energy supply (Equation 2) is calculated using two differential equations. The differential equation of exponential decay factor S damp = e −t/τ is expressed as dS damp /dt = -S damp /τ with an initial value of one. The part of energy supply that increases linearly with time, S lin = R t, is expressed as dS lin /dt = R, and the initial value is zero. Therefore, the energy supply of t is S = S damp S lin + S 0 .
Postsynaptic potential energy P is expressed as the integration of postsynaptic unit membrane power to time and is constrained by energy supply. The differential expression is as follows: where sign (·) is a symbolic function. The value of the function is −1 when the parameter is negative and is 1 when the parameter is positive. The value of this function is 0 when it is equal to 0. The unit of P is fJ/µm 2 , that is, 10 −15 J/µm 2 , and the initial value is 0. The unit of S and S 0 is the same as that of P, and the unit of R is fJ/(µm 2 s). Following the definition and considering the limitation of energy supply, the differential forms of baseline potential energy and suprathreshold potential energy are expressed in V m and I m as follows: where Θ(·) is the Heaviside step function. The function value is zero when the parameter is negative; otherwise, it is one. The initial values of baseline potential energy P bas and suprathreshold potential energy P sup are zero. After substituting Equations (4) and (5) into Equation (1), the differential form of synaptic weights can be expressed completely as follows: where The pseudocodes of our synaptic plasticity model is described in Python as follows: ###################################################### ###### Codes of parameter setting and variable initialization for synaptic plasticity model if "data from Sjöström": # in fig2b and fig3a A scale = 60/5 else: A scale = 1 A = 0.02 # amplitude coefficient A r = 0.2 # baseline coefficient V th = −60 * mV # threshold potential R = 175 # dimensionless-energy-supply rate τ = 2 * second # time constant of energy supply S 0 = 25 # dimensionless-minimum-energy supply W = 0.5 # initial weights of all synapses with energy supply constraints S damp = 1 # initial S damp of all synapses S lin = 0 # initial S lin of all synapses P = 0, P bas = 0, P sup = 0 # initial P, P bas , and P sup of all synapses with energy supply constraints W' = 0.5 # initial weights of all synapses without energy supply constraints P' = 0, P' bas = 0, P' sup = 0 # initial P, P bas , and P sup of all synapses without energy supply constraints ###### The main codes of our model dS damp /dt = -A scale * S damp /τ # differential form of equation S damp = e −t/τ dS lin /dt = A scale * R/second # differential form of equation S lin = R t dP/dt = A scale * sign(S damp * S lin +S 0 -abs(P)) * V m * I m1 /mV /second # codes of Equation 3, dP/dt = V m I m sign(S -|P|), I m1 is dimensionless I m dW/dt = A scale * sign(S damp * S lin +S 0 -abs(P)) * A * (A r * (V m <V th ) -(V m ≥V th )) * V m * I m1 /mV /mV /second # codes of Equation 6, dW/dt = A V m I m Φ(V m -V th ) sign(S -|P|) dP bas /dt = A scale * sign(S damp * S lin +S 0 -abs(P)) * A r * (V m <V th ) * V m * I m1 /mV /mV /second # codes of Equation 4, dP bas /dt = A r V m I m Θ(V th -V m ) sign(S -|P|) dP sup /dt = A scale * sign(S damp * S lin +S 0 -abs(P)) * (V m ≥V th ) * V m * I m1 /mV /mV /second # codes of Equation 5, dP sup /dt = V m I m Θ(V m -V th ) sign(S -|P|) dP ′ /dt = A scale * V m * I m1 /mV /second # codes of Equation 3 without energy supply constraints, dP/dt = V m I m dW ′ /dt = A scale * A * (A r * (V m <V th ) -(V m ≥V th )) * V m * I m1 /mV /mV /second # codes of Equation 6 without energy supply constraints, dW ′ /dt = A V m I m Φ(V m -V th ) dP ′ sub /dt = A scale * A r * (V m <V th ) * V m * I m1 /mV /mV /second # codes of Equation 4 without energy supply constraints, dP ′ bas /dt = A r V m I m Θ(V th -V m ) dP ′ sup /dt = A scale * (V m ≥V th ) * V m * I m1 /mV /mV /second # codes of Equation 4 without energy supply constraints, dP' sup /dt = V m I m Θ(V m -V th ) ######################################################

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found at: https://github.com/HWChenCSU/energydeterminant-plasticity.

AUTHOR CONTRIBUTIONS
HC designed the study, designed and implemented the models, and wrote the manuscript. LX analyzed the model and participated in discussions. YW and HZ discussed the results and commented on the manuscript. All authors contributed to the article and approved the submitted version.