Energetic substrate availability regulates synchronous activity in an excitatory neural network

Neural networks are required to meet significant metabolic demands associated with performing sophisticated computational tasks in the brain. The necessity for efficient transmission of information imposes stringent constraints on the metabolic pathways that can be used for energy generation at the synapse, and thus low availability of energetic substrates can reduce the efficacy of synaptic function. Here we study the effects of energetic substrate availability on global neural network behavior and find that glucose alone can sustain excitatory neurotransmission required to generate high-frequency synchronous bursting that emerges in culture. In contrast, obligatory oxidative energetic substrates such as lactate and pyruvate are unable to substitute for glucose, indicating that processes involving glucose metabolism form the primary energy-generating pathways supporting coordinated network activity. Our experimental results are discussed in the context of the role that metabolism plays in supporting the performance of individual synapses, including the relative contributions from postsynaptic responses, astrocytes, and presynaptic vesicle cycling. We propose a simple computational model for our excitatory cultures that accurately captures the inability of metabolically compromised synapses to sustain synchronous bursting when extracellular glucose is depleted.

Individual bursts are recorded per electrode as a list of individual burst timestamps and lengths, and the MC Rack output file is converted to ASCII format using the software MC Tools (Multichannel Systems).
The ASCII file output from MC Tools is used as an input for a custom-built synchronous burst detection algorithm written in the R programming language. Electrodes are pre-processed to detect outliers, which are identified and excluded from further analysis if they have a number of recorded individual bursts falling more than 1.5 times the interquartile range above the third quartile or below the first quartile, or zero. Time is binned into 100ms intervals and a logical vector spanning the length of the entire recording is initialized for each remaining electrode, containing value TRUE if an individual burst is recorded for that electrode during the given time interval and FALSE otherwise. An integer-valued burst profile vector, also spanning the length of the recording, is then constructed by summing the total number of electrodes individually bursting at each time interval. Run-length encoding of consecutive intervals within the burst profile vector where four or more electrodes are individually bursting simultaneously is used to mark the start and end of synchronous bursts, which at their peak typically involve 35-45 electrodes. Typical synchronous bursts last 600-800ms whilst individual bursts at electrodes usually have a length of 80-500ms. The total count of synchronous bursts over the course of a recording is divided by length of that recording (in mins) to calculate SBPM for further analysis. SBPM is capped at a maximum of 10 to avoid over-counting in hyper-active cultures.

Network and neural model
We employ the same models for the neural network and the neural dynamics as in the work of [2], which we recapitulate in this section, but modify the model of synaptic dynamics as described in the following section. As in the original paper, we simulate a network of 400 neurons, organized on a physical grid of 20x20, where the probability of connection is a decaying function of distance. The connectivity between neurons is defined by the binary connectivity matrix W, where W(i, i) = 0, d(i, j) returns the distance between the two units, and σ represents the connectivity decay length.
The neural model consists of a simplified Hodgkin-Huxley model described by: where V is the neuron's membrane potential, C is the membrane capacitance, and I app , I Na , I K , I L , I syn represent the external, sodium, potassium, leak, and total synaptic currents, respectively. ξ denotes a Wiener process, and σ ξ represents the standard deviation of the noise. In this simplified version of the model, only the dynamics of the potassium activation variable n are computed via: while the sodium activation variable m is assumed to instantaneously adapt to: The currents, are then calculated via: All parameter are presented in Table S1. No external current is applied, but frozen noise (constant in time but variable across units) of zero mean is provided to simulate heterogeneity in the neurons resting potentials (−65 ± 0.2mV). The calculation of the total synaptic current is presented after the introduction of the synaptic dynamics.

Metabolic-dependent synaptic dynamics
As in [2], we consider synapses endowed with synaptic facilitation, captured by facilitation variable x, whose dynamics evolve according to: where x increases after every action potential of the presynaptic neuron, when V is above the threshold T in the Heaviside function H, and decays to resting value X with time constant τ f . Beyond the use as a facilitation variable, functionally, x is generically employed in the model as a proxy for the level of recent neural activity, which becomes useful to evaluate both refractoriness and metabolic consumption, as we later describe.
In the original work of [2], three synaptic vesicle pools are simulated, representing the fraction of available free, docked, and recovering vesicles. Here we simulate only the fraction of docked vesicles, assuming that their recovery is simply modulated by the ratio of metabolic flux to neural activity. We have shown that, in absence of this metabolic coupling, by setting the other vesicle pools to a constant in the original model, the bursting behaviour remains unchanged, showing that the simulation of the dynamics of these two additional pools is not relevant for the purposes of the present work, and justifying this simplification. The docked vesicle pool (represented by y dock ), evolves according to: where represents the metabolic modulation function. We note here that we do not intend to capture the exact functional form of this coupling in nature (which is currently unknown), but simply to present a simplified model with the desired properties: the exponential ensures the sign of the function remains positive, becoming smaller as the ratio of recent past energy consumption to resource availability decreases. This ratio is here expressed as the proxy for recent neural activity x over the metabolic flux m φ , which we assume to be a function of glucose availability. High recent activity means here a larger resource depletion, and reduced capacity for y dock to recover, while a higher metabolic flux results in an increased capacity for vesicle recovery.

3/9
Finally, as in [2] we consider synapses to become refractory after persistent firing of the presynaptic neuron, when y dock reaches a minimal depletion variable y min dock , and stays refractory until this value has recovered above this minimal value and x has has also reached a reference value x re f . We note that the equations for x and for the synaptic pools in this work and in the original model are only expressed in terms of global properties of the presynaptic neuron. All synapses of the same presynaptic neuron follow the same dynamics and can be simulated as one. In particular, all synapses corresponding to the same presynaptic neuron will become refractory at the same time. Refractory synapses do not contribute to the total synaptic current each postsynaptic neuron receives. If a spike arrives at time t 0 to a non-refractory synapse j, it produces a current given by: which integrates the number of depleted docked vesicles and includes a proportionality constant K I . The two Heaviside functions ensure that currents are only produced during spikes and only if the fraction of docked vesicles has not dropped below a minimal value. The total current is then computed as: where a delay time t del = 1ms is employed.

Simulations
Simulations were performed in python, employing an adaptation of the original matlab code of [2]. The full code for the simulations here presented can be found at https://github.com/r-echeveste/meta_ reg_net.
Network dynamics were simulated for 100s, after an initial 5s transient. The evolution of the system was computed using a fourth-order Runge-Kutta method, with dt = 0.05ms. Simulations were conducted in two different metabolic flux conditions: a high metabolic flux case, corresponding to high levels of extracellular glucose, and a low metabolic flux case, corresponding to low levels of extracellular glucose. The parameters for both simulations are found in Table S1.