Abstract
Introduction:
Computational models are valuable tools for understanding and studying a wide range of characteristics and mechanisms of the brain. Furthermore, they can also be exploited to explore biological neural networks from neuronal cultures. However, few of the current in silico approaches consider the energetic demand of neurons to sustain their electrophysiological functions, specifically their well-known oxygen-dependent firing.
Methods:
In this work, we introduce Digitoids, a computational platform which integrates a Hodgkin-Huxley-like model to describe the time-dependent oscillations of the neuronal membrane potential with oxygen dynamics in the culture environment. In Digitoids, neurons are connected to each other according to Small-World topologies observed in cell cultures, and oxygen consumption by cells is modeled as limited by diffusion through the culture medium. The oxygen consumed is used to fuel their basal metabolism and the activity of Na+-K+-ATP membrane pumps, thus it modulates neuronal firing.
Results:
Our simulations show that the characteristics of neuronal firing predicted throughout the network are related to oxygen availability. In addition, the average firing rate predicted by Digitoids is statistically similar to that measured in neuronal networks in vitro, further proving the relevance of this platform.
Dicussion:
Digitoids paves the way for a new generation of in silico models of neuronal networks, establishing the oxygen dependence of electrophysiological dynamics as a fundamental requirement to improve their physiological relevance.
1 Introduction
Exploring how neurons process and transmit information is crucial for advancing our knowledge of the brain. Along with the study of biological neural networks in cultures or in in vitro slices (; ; ; Van Pelt et al., 2005), computational, or in silico, models have been successfully exploited, e.g., to support the study of neuronal network modulation and delineate potential mechanisms underlying activity patterns (; Masquelier and Deco, 2013; Sukenik et al., 2021; Wen et al., 2022). Several model-based solutions for generating virtual representations of neural cells able to replicate the salient properties of experimentally observed behaviors have been proposed (Lonardoni et al., 2015). For instance, intuitive and easy to use simulators (e.g., BRIAN 2, NEST, NEURON) have been employed to simulate spiking neural network models (; ; Stimberg et al., 2019). Traditionally, they include mathematical descriptions of the single-neuron activity, ranging from simple phenomenological characterization of neuronal spiking (Izhikevich, 2003) to more complex, biophysical conductance-based simulations of ion fluxes between the intra and the extracellular space (), as well as models of cell-cell connections to replicate the neuronal network architecture (Markram et al., 2015; Masoli et al., 2022; Potjans and Diesmann, 2014). Some studies also incorporate more sophisticated models, e.g., including astrocytes via tripartite synapses (Lenk et al., 2020).
However, few of these approaches include energetic considerations, i.e., the dynamics of ATP hydrolysis (Kuznetsov, 2024; Wei et al., 2014). It is well-known that metabolism is involved in brain functionality: nutrients—and, in particular, oxygen (O2)—fuel brain specialized functions, determining the electrophysiological dynamics and brain plasticity, up to cognitive functions (Watts et al., 2018). More specifically, beyond the basic activities common to other cells (e.g., DNA and RNA synthesis), resource uptake in neurons is also dedicated to support spiking, because of the role of the Na+-K+-ATP pump in signal propagation (; Lennie, 2003). Since ATP dephosphorylation depends on the rate of O2 consumption, its dynamics can be monitored (; Özugur et al., 2020). O2 dependence is also crucial for in vitro slice preparations, requiring humid and well-oxygenated environment for their culturing (Sanchez-Vives et al., 2000). An analytical formulation describing O2-dependent firing was proposed by Wei and collaborators (Wei et al., 2014) to elucidate the mechanisms of seizure development and termination, as well as their interaction with energy metabolism. This model assumes that O2 variations depend on the diffusion from the bath solution and on the neuronal consumption rate for firing, but it does not consider that O2 can also be consumed for sustaining other metabolic functions of the cell (; Lennie, 2003).
The formulation proposed by Wei’s team was applied to brain tissue slices. However, O2 is also crucial in in vitro cultures: for example, in traditional monolayers, cells are inevitably exposed to different O2 levels when varying the amount of medium or the O2 boundary concentration (; ; Pacitti et al., 2019; Walsh et al., 2005). Starting from Wei et al.’s model, we have developed a computational platform able to mimic the in vitro electrophysiological behavior of neuronal cultures at the single-cell and network level. We refer to Digitoids as the digitalized versions of in vitro neuronal monolayers obtained from dissociated neurons, in which the dependence on O2 concentration of network dynamics is considered. As in vitro networks can have different culture conditions and layouts (; ; ; ), the platform is purposely designed to be modular, thus the user can generate Digitoids matching any type of in vitro neuronal network. Here we describe the theory and computational setup of Digitoids. For testing the performance and highlighting the crucial role of O2 in describing firing dynamics in neuronal cultures, we digitalized the layouts of neuron monolayers seeded on commercial micro-electrode arrays (MEAs). The O2-dependent model of firing and metabolism was implemented on digitalized networks to assess if a degree of similarity can be found between the Digitoids’ output and the corresponding experimental data from MEA recordings, comparing the predictivity of our platform to that of traditional models which neglect the dependence of firing activity on O2 supply. Albeit preliminary, these results highlight the significant role of O2 dynamics in network behaviors and thus the necessity of including energetic considerations while mathematically describing electrophysiological activity in cell cultures.
2 Materials and methods
2.1 Theory and outline of the computational platform
Figure 1. A shows the in vitro scenario simulated by the Digitoids. It is composed of a well seeded with neurons, supplied with a layer of culture medium of height h. The cells are assumed to be homogeneously distributed on the bottom of the well (at z = 0). Four phenomena occur in the system: i) O2 diffusion through the medium, ii) O2 consumption by neurons to fuel both basic cellular processes and electrophysiological activity, iii) neuron firing and iv) neural network dynamics, i.e., the transfer of electrical information via synaptic-mediated connections among cells. Considering the symmetry of the system, O2 diffusion can be assumed to occur only along the z axis and independently of the x and y directions (McMurtrey, 2016; Patterson and Mazurek, 2010; Place et al., 2017). Each neuron at z = 0 consumes O2 as described in the subsection 2.1.2 Single-neuron model, generating an axial concentration gradient and a consequent downward flux. Moreover, O2 diffusion and reaction can be simulated as “background dynamics”, given that their characteristic times are significantly longer than those of the electrophysiological phenomena occurring on the x,y plane, where the neuron monolayer lies (Table 1). Transfer information is mainly influenced by the strength and number of synaptic connections between neurons. Thus, the O2-dependent single-neuron dynamics can be decoupled from those of the network as a whole. As such, the network (Figure 1A) can be considered as the integration of modules describing the O2 consumption—depending on its downward diffusion—as well as the firing for a single neuron (Figure 1B), modulated through the extent of its in-plane connectivity. Under these assumptions, the single-neuron metabolic and electrophysiological activity can be determined at each time step according to the O2 concentration perceived by the cells at z = 0, which is in turn updated depending on the single-neuron consumption and allows computing the diffusive flux magnitude along the medium column. On this basis, the O2 concentration profile at the subsequent time step can be estimated and the process iterated over time.
TABLE 1
| Phenomenon | Characteristic time (s) | References |
| O2 diffusion in z direction | ∼ 0.6 | (Magliaro et al., 2019) |
| O2 diffusion on x-y plane | ∼ 6 104 | (Magliaro et al., 2019) |
| O2 consumption | ∼ 31 | (Magliaro et al., 2019) |
| Neuron firing | ∼ 10–3 | () |
| Synapses | ∼ 2 10–4 | (; Wang et al., 2010) |
Characteristic times of the phenomena involved in the single-neuron model: O2 diffusion—in the x direction and in the x-y plane, O2 consumption, neuron firing and synapses.
From the evaluation of such characteristic times, it was possible to assume that single-neuron dynamics are decoupled from the network ones.
FIGURE 1
The computational platform was developed in Matlab (version R2023b the Mathworks Inc., Boston Massachusetts), exploiting the Simulink toolbox.
2.1.1 Diffusion model
O2 diffusion through the culture medium is modeled as a one-dimensional phenomenon governed by the Fick’s second law:
where c (mol m–3) is the O2 concentration and D (m2 s–1) is the diffusion constant of O2 in the culture medium. Eq. (1) is solved using the finite difference method according to the initial and boundary conditions. Specifically, assuming that the well is initially filled with O2-saturated culture medium, the initial condition is c(z, 0) = c0, and the air-medium interface maintains a uniform and time-invariant O2 concentration, i.e., c(h,t) = c0. Note that, as neurons consume O2 by means of a surface reaction (i.e., they sink O2 as an outward flux through the well bottom), there is no volumetric reaction term to include in Eq. (1).
2.1.2 Single-neuron model
The single-neuron model describes the O2 consumption for maintaining both vital and electrophysiological functions and the O2-dependent dynamics in each cell of the network as a function of the current O2 availability [i.e., c(0,t)] as input.
Regarding the O2 consumption, we assume that 75% of the available O2 is devoted to fuel neuronal spiking activity (namely, cf = 0.75⋅c), and the remaining 25% (namely, cnf = 0.25⋅c) to sustain basic cell processes (; Lennie, 2003). The O2 consumption rate of the whole neuron network (R(c), in mol m–3 s–1) can be thus expressed as:
where Rnf(cnf) is the rate at which O2 is consumed for cellular and sub-cellular processes not directly linked to electrophysiological activity, and Rf(cf) is the O2 consumption contributing to neuron firing. Specifically, Rnf(cnf) can be formulated according to the Michaelis-Menten kinetics (; Magliaro et al., 2019):
where sOCR (mol s–1) is the maximal consumption rate of a single cell in the network, ρcells (m–3) is the volumetric cell density of the monolayer and km (mol m–3) is the Michaelis-Menten constant, i.e., the O2 concentration corresponding to half saturation of the consumption rate. On the other hand, Rf(cf) was described by Wei and co-workers (Wei et al., 2014) as:
where Ipump (mol m–3 s–1) is the transport rate of ions across the membrane and α (a.u.) is a conversion factor from pump transport rate to time variation of O2 concentration. Ipump is related to intracellular (subscript i) sodium and the extracellular (subscript o) potassium concentrations as in the following equation:
in which we assume that the rate ρ (mol m–3 s–1) at which the pumps transport ions across the membrane depends on the O2 concentration according to a sigmoidal function.
In Eq. (6), ρmax (mol m–3 s–1) is the maximal rate at which the pump operates, i.e., when the medium is fully oxygenated. Therefore, Ipump regulates the trans-membrane electrochemical gradient depending on the O2 availability, which thus influences the membrane potential V (mV) and the firing activity of the neuron. The Hodgkin-Huxley (HH) model is used to describe the dynamics of V (; ):
where C (μF cm–2) is the membrane capacitance, Iext (μA cm–2) is the external applied or synaptic current from other neurons, INa,IK,ICl (μA cm–2) are the sodium, potassium and chloride currents. The latter corresponds to a leakage current, as it is mainly represented by flux of Cl– ions ().
It is worth highlighting that—as in the traditional formulation of the HH model—the membrane potential V (Eq. 7) depends on the potassium and sodium currents INa and IK:
where m,p, and n are activation and inactivation variables (their description is given by Supplementary Eqs. 1–7) of voltage-gated ionic channels, whose values range from 0 to 1 and define the fraction of open and closed channels throughout the membrane. For the sake of simplicity, the non-voltage-sensitive leaks were not reported. As detailed in Eqs. (8) and (9), INa and IK are in turn functions of the reversal potential ENa and EK, respectively, given by the Nernst equation:
However, while in the HH model the intracellular concentration of sodium and the extracellular concentration of potassium are considered as constants, in this formulation they are modulated by Ipump, which is a function of the local O2 concentration, as described through Eqs. (5) and (6). Thus, Nernst potentials of sodium and potassium (Eqs. 10 and 11) vary with O2. All the dynamics describing neuronal functioning are here assumed to occur at 37°C, corresponding to the physiological temperature for eukaryotic cells. Intracellular sodium and extracellular potassium concentrations are in turn modulated by INa, IK and Ipump, as described by the following equations (Eqs. 12 and 13):
More details on the model are provided in the Supplementary Text 1 (Eqs. 8–10).
2.1.3 Connectivity model
The neuronal network is generated connecting the single neurons. In this study, we implemented the neuron-to-neuron coupling via chemical synapses (Roth and van Rossum, 2009). Thus, the membrane potential of the i-th neuron is described by the following equation.
in Eq. (14) is the synaptic current input to the post-synaptic neuron i and it is modeled as:
in which we assume that the i-th neuron receives inputs from N pre-synaptic neurons. aij is the coefficient describing the connection between vertices i and j of the adjacency matrix A, obtained through the Watts-Strogaz method (more details in the next Section and in Supplementary Table 2). (mV) is the reversal potential of the synapse for the j-th pre-synaptic neuron and can assume the following values according to the nature of the synaptic connection (; Wei et al., 2014).
The value of the synaptic conductance (μS cm–2) is modified every time the pre-synaptic neuron fires, i.e., every time Vi exceeds the threshold value of 0mV with a positive derivative. At each spike, there is a release of neurotransmitter into the synaptic cleft, thus the synaptic conductance over time is modeled as an exponential decay:
where t0 is the time at which the spike is fired by the pre-synaptic neuron, is the maximal conductance value and τsyn is the decay time constant, which assumes the following values (Wei et al., 2014).
The synaptic dynamics are implemented in the model by updating the value of the synaptic conductance as follows (; Roth and van Rossum, 2009):
where μS cm–2 is the intensity of the synaptic update, the same for both excitatory and inhibitory synapses. In Digitoids, 80% of neurons are excitatory and 20% inhibitory.
2.1.4 Network model
It has been observed that the structure of neuronal networks in both brain tissues and cellular monolayers can be described by Small-World (SW) graphs (; ; ). Specifically, a SW graph shows intermediate characteristics between a random and a regular graph, with dense clustering of neighboring vertices and short distances between pair of vertices. Indeed, in vivo chemical synapses typically facilitate the formation of dense local connections between neurons, thus giving rise to clusters, as well as of long-range connections allowing clusters of neurons to communicate (). Given a network composed of n vertices and m edges, it can be described by the metrics reported in Supplementary Table 2 (; Watts and Strogatz, 1998).
We generated SW neural networks in a purposely developed Simulink library, which describes the wiring information through an adjacency matrix A, usually used to represent inter-neuron connections (; Poli et al., 2015; Shefi et al., 2002). Starting from the number of vertices and edges and the metrics characterizing the networks, A can be obtained using the Watts-Strogatz method (). Each coefficient of the matrix describes the connectivity between vertex i and j. Specifically, aij = 1 if an edge exists from vertex i to vertex j, otherwise it is 0. We thus exploited such adjacency matrices to create connections between neurons, defined by chemical synapses (Eqs. 14–20). Both the O2 diffusion and the single-neuron models were integrated in the library, which allows defining: (i) the initial and boundary O2 concentration c0 at the air-medium interface, (ii) the height of the medium h, and (iii) the metabolic and firing parameters of the neuron.
2.2 Impact of oxygen on single-neuron activity
For assessing the influence of O2 availability on firing, the single-neuron model coupled with O2 diffusion was first computed using stepwise variations of both (i) the boundary concentration of O2 from 0.2 mol m–3 (i.e., the maximum available oxygen concentration in water) to 0.04 mol m–3 [i.e., the critical oxygen concentration for cell survival ()] and (ii) the culture medium height h from 0.1 to 3 mm, based on the conditions usually used for neuron electrophysiological recordings (; Negri et al., 2020; Scelfo et al., 2012). All the parameter combinations were simulated for 20 s (variable step solver “ode15s” by Simulink) and are summarized in Table 2.
TABLE 2
| h (mm) | c0 (mol m–3) |
| 3 | 0.2 |
| 2 | 0.16 |
| 1 | 0.12 |
| 0.5 | 0.08 |
| 0.1 | 0.04 |
Values of the parameters simulated in the single-neuron configuration.
Every combination of the two parameters – medium height h and boundary O2 concentration c0−was tested, for a total of 25 configurations in the single-neuron model, to assess and characterize the influence of these parameters in shaping the resulting electrophysiological activity.
2.3 Analysis of single-neuron membrane potential
To characterize how the shape of the spike trains and the single-spike waveforms are influenced by the different combinations of c0 and h–and, thus, by the overall O2 availability within the system—we defined two new metrics: the Aspect Ratio (AR, expressed in logarithmically-scaled mV s–1) and the Dissipation Rate (DR, expressed in s–1), defined as follows:
where ΔVmax is the peak-to-peak amplitude of the highest spike in the train, ttrain is the time duration of the train and α (in mV s–1) is the average value of the first derivative of the envelope of the peaks in the train. ttrain was expressed as the difference between the end and start times tend and tstart, identified as the time at which the first derivative of the signal is equal to 0 and the time at which the signal amplitude > −60 mV (; Wilson and Emerson, 2002), respectively. Figure 2A reports a typical spike train, and a graphical representation of the quantities used to calculate AR and DR.
FIGURE 2
We separately assessed the correlation of each of the three metrics–ttrain, AR and DR—with the boundary O2 concentration c0 and the medium height h by computing the non-parametric Spearman coefficient (significance level of 0.05).
The shape of single spikes was also evaluated, calculating the peak-to-peak amplitude (vpp = vmax−vmin, expressed in mV), rise rate (rr = (vmax−vstart)/(tmax−tstart), in mV s–1) and fall rate (fr = (vmax−vend)/(tmax−tend) in mV s–1), where tmax is calculated as the time corresponding to the maximum of the spike (Figure 2B; ; Zaitsev et al., 2012).
Finally, to describe the features of the spikes fired by single neurons as a function of the balance between diffusive O2 supply and its consumption by the neurons irrespective of the specific setup of the simulation, we exploited the Thiele Modulus, Φ2. Specifically, Φ2 is defined as the ratio between the characteristic diffusion (τd) and reaction (τr) times. Since metabolism and firing occur simultaneously in the neuron domain, the reaction dynamics is driven by the faster of the two phenomena. Given that the reaction is described by the sum of two rates (Eq. 2), Φ2 can be formulated as follows:
where τnf and τf indicate the characteristic times of basal and firing-related O2 consumption, respectively. Refer to Supplementary Text 2 for further details on the derivation of Eq. (23). The shape metrics of the spike trains–AR and DR−were then also evaluated as a function of Φ2.
2.4 Assessment of digitoids performance
2.4.1 Digitoids versus experimental data
Digitoids performance was evaluated using the experimental data presented in , following the pipeline shown in Figure 3. In ref. (), the authors describe the morphology and the electrophysiological activity of neuron networks in vitro. The network activity was recorded via a MEA, and the mean Firing Rate (mFR) as well as the event synchronization were extracted. The topological evolution of the networks was mapped to a network graph, where neurons are represented as vertices and their physical connections as edges, and the SW metrics were defined. Their experimental setup and SW metrics are detailed in Supplementary Tables 1, 2, and a more in-depth description of the experimental set-up and procedures is provided in Supplementary Text 4. In our work, measurements from Day In Vitro (DIV) 11 to DIV 16—i.e., when the network exhibits a SW layout () —were exploited, without the intention of mapping the temporal evolution of the in vitro neuronal cultures. Given the number of vertices and edges for those DIVs reported in Ballesteros-Esteban and co-workers and the metrics characterizing the networks, the adjacency matrix was obtained through the Watts-Strogatz method, setting the rewiring probability to 0.5 (Watts and Strogatz, 1998). Three SW graphs were obtained for each DIV considered. The outcoming connectivity models are sparse (i.e., the number of edges is less than the possible number of edges in the order of O(q), where q is the total number of vertices), with a mean edge density (defined in Supplementary Table 2) of 2.5%, in consistence with previously reported experimental values (; ).
FIGURE 3
The SW layouts and the adjacency matrices were used to generate the corresponding Digitoids. The layouts, along with their number of vertices and edges, are reported in Supplementary Table 3, while the model parameters are listed in Table 3. The same SW layouts were used to build in silico neuronal networks where the traditional HH model was implemented instead of the O2-dependent one, described in Section 2.1.2. The single-neuron description was obtained by imposing the membrane pump to work optimally, i.e., with fixed pump rate ρmax (Table 3). The current components of the model are the same of the single-neuron model (Section 2.1.2)—i.e, INa, IK and ICl−−consistently with the model developed by Wei and co-workers (Wei et al., 2014). The neurons in the computational network models are spontaneously active due to potassium concentration in the bath. All the network models were simulated for 20 s with the variable-step solver “ode15s” of Simulink, with maximal step size of 0.4.
TABLE 3
| Model parameter | Value | References |
| Diffusion constant (D) | (McMurtrey, 2016) | |
| Oxygen Consumption Rate per cell (sOCR) | () | |
| Cell density (ρcells) | () | |
| Michaelis-Menten constant (km) | () | |
| Conversion factor from pump current to oxygen concentration (α) | 0.17 | (Wei et al., 2014) |
| Conversion factor current to concentration (γ) | (Wei et al., 2014) | |
| Ratio to intra/extracellular volume (β) | 7 | (Wei et al., 2014) |
| Maximal Na-K pump rate (ρmax) | (Wei et al., 2014) | |
| Membrane capacitance (C) | 1μF/cm2 | (Wei et al., 2014) |
| Maximal sodium conductance (GNa) | 30mS/cm2 | (Wei et al., 2014) |
| Maximal potassium conductance (GK) | 25mS/cm2 | (Wei et al., 2014) |
| Reversal potential of synapses, | 0mV, if excitatory | (; Wei et al., 2014) |
| −80mV, if inhibitory | ||
| Time constant of synapses, τsyn | 4ms | (; Wei et al., 2014) |
| 8ms | ||
| Synaptic update, | 0.5μScm−2 | () |
Values of the parameters used in the model.
This table reports all the constants adopted in the model described in this work.
2.4.2 Impact of oxygen on network-level activity
Six Digitoids (with SW layout size described in Supplementary Table 3) were developed and simulated to explore the effects of O2 deprivation on the network activity. For this purpose, the six Digitoids were first simulated in normal oxygenation conditions for cell culture, i.e., considering a boundary concentration c0 = 0.2mM. Then, the same networks were simulated lowering c0 to 0.04mM, i.e., the threshold O2 concentration ensuring physiological cell functioning and survival ().
2.5 Statistical analysis
Statistical analyses were performed using GraphPad Prism 8 (GraphPad Software, Boston, Massachusetts United States) to identify any significant differences between the mFR of the computational models and the experimental data. Thus, firstly, the distributions of mFR of the O2-dependent firing in Digitoids, the mFR experimentally measured in cultured neurons and the mFR values from the traditional HH model were tested for normality, by adopting the Shapiro-Wilk test (α = 0.05). Since the distributions were not Gaussian, the non-parametric Kruskal-Wallis test was used (α = 0.05). To compare mFR and event synchronization between the Digitoids simulated in normal (i.e., c0 = 0.2mM) and O2 deprivation (i.e., c0 = 0.04mM) conditions, the Mann-Whitney test was instead adopted (α = 0.05).
3 Results
3.1 Dependence of firing on oxygen availability
Figure 4 shows examples of the outcome of the Digitoids, i.e., the neuron membrane potential and the O2 concentration at the cell level (z = 0) taken over a time window of 20 s for different values of the boundary O2 concentration. As expected, the plots indicate that single neurons exhibit an O2-dependent firing, with reduced activity when the local concentration decreases. Indeed, when the neuron fires, the Na+-K+-ATP pump is activated, thus O2 is consumed (Eqs. 2–4), and its concentration at z = 0 decreases. Longer spike trains are generated if O2 availability is high.
FIGURE 4
For what concerns the sensitivity of the shape metrics to the parameters c0 and h, plots are reported in Supplementary Figures 5, 6. Specifically, Supplementary Figure 5 graphically depicts the dependence of ttrain, AR and DR (Eqs. 21 and 22) on c0 for each of the tested medium heights, while Supplementary Figure 6 reports their dependence on h parametrized with respect to c0. From the visual analysis of the plots, a monotonic relation can be identified between the parameters c0 and h–which set the availability of O2 over time to the neuron—and the train metrics. This suggests that the neuron is able to fire longer trains of action potentials when the O2 availability in the system is not a limiting factor, i.e., with highest c0 and lowest h. Furthermore, the Spearman correlation coefficient r was computed to provide a quantitative means of such dependencies. Numerical values of r are reported in Supplementary Tables 4–9 together with corresponding p-values. All the metrics display significant correlation with c0, while they significantly correlate with the medium height only when boundary O2 is maximal. Indeed, the single-neuron output is more sensitive to growing medium heights when O2 availability is not limited yet by reduced air saturation, that is c0 < 0.2 mM. Otherwise, supply constraints due to the increased diffusive path do not significantly affect the duration of spike trains.
Moreover, single spikes were identified for each combination of h and c0; for each spike, vpp, rr and fr were calculated, and their trend over time are shown in Supplementary Figures 1–3. At the beginning of the simulation (i.e., when O2 availability is high), vpp values appear independent of h and c0 (see Supplementary Figure 1). Then, vpp decreases over time with a rate depending on c0. In particular, we observed that the rate at which vpp decreases at the end of the spike train is higher for the lower values of c0. This is reported in Figure 5A, where the slope of vpp (dVpp) over the last four time points considered in the simulation is shown to better highlight the dependence on the different values of h and c0.
FIGURE 5
Figure 5B depicts ttrain as a function of Φ2. Notably, ttrain is sensitive to the level of O2 available to the neuron, as reported in , since it decreases with higher Φ2 (that is with lower c0 and higher h). This implies that firing is a diffusion-limited phenomenon, which is suppressed when it cannot be energetically sustained due to O2 depletion (Nieber, 1999; Pires Monteiro et al., 2021; Santiago et al., 2023). Moreover, the dispersion of ttrain values becomes narrower with increasing Φ2, indicating that the firing threshold is governed by O2 availability, which is in turn increasingly limited by diffusion as h increases.
The same trends with respect to Φ2 are observable for DR and AR−see Supplementary Text 3 and Supplementary Figure 4 for details.
3.2 Performance of the digitoids
Figure 6 shows the mFR obtained from: (i) the experimental in vitro recordings reported in ; (ii) the output from the corresponding Digitoids; (iii) the firing activity of the network with the same topological layout but the traditional formulation of electrophysiology according to the HH model. For all the DIV considered, no statistically significant differences were found between the mFR of Digitoids and the corresponding experimental data. On the other hand, the values of mFR of the traditional O2-independent HH model were significantly different if compared to both the in vitro observations and Digitoids predictions. The associated p-values are reported in Supplementary Table 10.
FIGURE 6
Further, the whole-network effect of O2 deprivation on mFR predicted by Digitoids is shown in Figure 7A. When accounting for reduced O2 availability (c0 = 0.04mM), Digitoids coherently predicted significantly lower mFR than that obtained for c0 = 0.2mM. The event synchronization was also evaluated in such conditions. Also in this case, O2 deprivation lowers the predicted synchronization values, with significant differences with respect to values predicted by the Digitoids with normoxic conditions (Figure 7B).
FIGURE 7
Supplementary Figure 7 depicts an example of the event synchronization calculated from one of the simulated Digitoids.
4 Discussion
O2 levels are crucial to neuronal function in vitro: they significantly affect viability, oxidative stress and mitochondrial function (Zhu et al., 2012). However, the influence of O2 on in vitro electrophysiological behavior is often neglected. In this work, we developed a computational platform—Digitoids—able to replicate a neuronal network in vitro. Digitoids embeds a model of neuron firing in which the O2 dynamics of diffusion and consumption are introduced and coupled with ionic transport across the cell membrane. The novelty of the proposed model resides in the coupling of O2 diffusion and consumption dynamics with neuronal electrical activity. Thanks to this approach, different culture conditions and layouts can be replicated obtaining descriptions of O2-dependent activity tailored on the specific system under study.
To demonstrate the importance of O2 in neuron firing, we computed different metrics of the spike train as well as of single spikes and assessed their dependency on O2 availability. Overall, the observed trends confirm that the electrophysiological behavior of single neurons is modulated by O2 supply. These results are supported by the significant correlation between the train metrics and the boundary O2 concentration, c0. Interestingly, neuron firing was found to be less sensitive to O2 fluctuations in conditions of limited resource availability (i.e., for high Φ 2 values). Indeed, reduced—or even non-significant—correlation coefficients of spike train characteristics with medium height are found when boundary O2 does not correspond to air saturation (i.e., c = 0.2 mM).
This behavior can be explained considering that reduced O2 hinders the homeostatic maintenance of ion concentrations between the intra and extracellular environments, which is responsible for sustaining the electrical activity of the neuron, as reported for both brain slices and in vitro cultures exposed to hypoxia (; ; Pires Monteiro et al., 2021; Spong et al., 2016; Zanelli et al., 2015). Under these conditions, the Na+-K+-ATP pump lacks sufficient resources to fuel ion transport, and thus firing decreases or even ceases (Nieber, 1999). The preliminary results obtained simulating O2 deprivation at the network level corroborate this evidence, suggesting that cells reduce their electrical activity and synchronization than in standard oxygenation at both the single-neuron and whole-network scale. These results are consistent with studies which reported reduced firing rate of cultured neurons when exposed to hypoxia (; ).
As a first preliminary assessment of the goodness of Digitoids predictions, we compared the simulated firing rate to that measured in neuronal networks seeded on commercial MEAs. No statistically significant differences were found between the experimentally measured mFRs and those predicted by Digitoids. Additionally, we compared the mFRs observed in vitro to predictions by the classic HH model applied to the same network layouts. The results are significantly different, highlighting that the mutual influence between local O2 concentration and the ion pump activity affects electrophysiological dynamics, as also captured by the analyses performed on the single-neuron output. It is worth highlighting that for Digitoids and experimental data the mFR is much lower than in the traditional HH model. In the latter, the initial values of simulation parameters (described in Section 2.4) indeed induce neurons to fire longer trains of APs, which are not limited by reduced O2 availability. Including O2 dynamics instead allows Digitoids to mimic its potential deprivation due to neuronal uptake, hindering the cross-membrane transport of ions as the energetic demand of membrane pumps cannot be satisfied.
The platform is designed to be modular and adaptable to different culture conditions by tuning the cell metabolic and electrophysiological parameters. More complex models of neuronal cultures, (e.g., co-cultures) can be developed by adding different single-cell blocks to the Simulink library to mimic other neural phenotypes. Following the approach described in , three-dimensional (3D) neuronal constructs can also be built by overlaying monolayers on each other. Thus, Digitoids can be extended to the simulation of neurospheres and cerebral organoids (Poli et al., 2019), supporting the investigation of their biophysical mechanisms. These constructs are particularly susceptible to O2 availability, as its depletion can lead to the formation of non-viable cores, hindering the development of mature traits and of a 3D neural network (Poli et al., 2019).
It is important to note that the experimental data used for comparison were derived from insect neurons, while the model parameters are typical of mammalian neurons. Nevertheless, the fundamental mechanisms underlying spike generation are similar across different species (Spong et al., 2016), and invertebrates are widely used to advance our understanding of more complex organisms (Newcomb et al., 2023; Sattelle and Buckingham, 2006). A more in-depth validation of our platform would require parallel recordings of electrophysiological and O2 dynamics in in vitro neurons. Specifically, perturbations will be added in the model and the predicted output will be compared to an experimental setting where the same perturbation is introduced (e.g., incubator O2 level drop). Furthermore, future effort will be carried out to integrate in the model also the O2 demand of synaptic activity (). For what concerns the parameters specific of the electrophysiological model, they will be tuned to better fit the recorded electrophysiology of in vitro neurons.
To proof the feasibility of using Digitoids, we exploited topological and electrophysiological data acquired on low density cultured networks. Thus, the simulated networks involve a relatively limited number of neurons and connections. To further expand the relevance of this work, larger networks can be developed and simulated. Larger-sized Digitoids can be developed with the same approach described in this work (Section 2.4.1) by defining a bioinspired (i.e., based on biologically observed features) adjacency matrix to layout the spatial distribution of single neurons and their connections within the Simulink framework.
In conclusion, this work represents a promising first step towards creating “digital twins” of in vitro neuronal networks. The approach implemented in Digitoids can be exploited for gaining important insights into brain pathophysiology. As an example, stroke and ischemia are characterized by low O2 levels, which lead to cognitive decline, neuronal damage and cell death (Klein Gunnewiek et al., 2020; Radenkovic et al., 2024; Voogd et al., 2023). In addition, neurodegeneration is known to be intimately linked to mitochondrial—and thus bioenergetic—dysfunctions (). Hence, Digitoids hold the potential to support, or even replace, primary neuronal cultures, as they are cost-effective, have a longer lifespan and allow high-throughput experiments that would be unfeasible in vitro (Velasco et al., 2020). Ongoing efforts include further model validation through detailed O2 and electrophysiological measurements, as well as expanding the model to include additional modules for different neuron types and 3D networks.
Statements
Data availability statement
The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding authors.
Ethics statement
The manuscript presents research on animals that do not require ethical approval for their study.
Author contributions
RF: Data curation, Formal Analysis, Software, Visualization, Writing – original draft. EB: Data curation, Formal Analysis, Visualization, Writing – review & editing. AA: Supervision, Writing – review & editing. CM: Conceptualization, Funding acquisition, Methodology, Project administration, Supervision, Writing – review & editing.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This work has received funding from the NAP project (HORIZON-EIC-2022-PathfinderOpen) under grant agreement No. 101099310.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The authors declare that no Generative AI was used in the creation of this manuscript.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fninf.2025.1549916/full#supplementary-material
References
1
Al-AniA.TomsD.KondroD.ThundathilJ.YuY.UngrinM. (2018). Oxygenation in cell culture: Critical parameters for reproducibility are routinely not reported.PLoS One13:e0204269. 10.1371/journal.pone.0204269
2
AntonelloP. C.VarleyT. F.BeggsJ.PorcionattoM.SpornsO.FaberJ. (2022). Self-organization of in vitro neuronal assemblies drives to complex network topology.eLife11:e74921. 10.7554/ELIFE.74921
3
AshwinP.CoombesS.NicksR. (2016). Mathematical frameworks for oscillatory network dynamics in neuroscience.J. Math. Neurosci.6:92. 10.1186/S13408-015-0033-6
4
AttwellD.LaughlinS. B. (2001). An energy budget for signaling in the grey matter of the brain.J. Cereb. Blood Flow Metab.211133–1145. 10.1097/00004647-200110000-00001
5
Ballesteros-EstebanL. M.LeyvaI.AlmendralJ. A.Sendiña-NadalI. (2023). Self-organization and evolution of structure and function in cultured neuronal networks.Chaos Solitons Fractals173:113764. 10.1016/j.chaos.2023.113764
6
BassettD. S.BullmoreE. (2006). Small-world brain networks.Neuroscientist12512–523. 10.1177/1073858406293182
7
BeanB. P. (2007). The action potential in mammalian central neurons.Nat. Rev. Neurosci.8451–465. 10.1038/nrn2148
8
BergerE.MagliaroC.PacziaN.MonzelA. S.AntonyP.LinsterC. L.et al (2018). Millifluidic culture improves human midbrain organoid vitality and differentiation.Lab Chip183172–3183. 10.1039/C8LC00206A
9
BettencourtL. M. A.StephensG. J.HamM. I.GrossG. W. (2007). Functional structure of cortical neuronal networks grown in vitro.Phys. Rev. E - Statist. Nonlinear Soft Matter Phys.75:021915. 10.1103/PhysRevE.75.021915
10
BorgesF. S.ProtacheviczP. R.SouzaD. L. M.BittencourtC. F.GabrickE. C.BentivoglioL. E.et al (2023). The roles of potassium and calcium currents in the bistable firing transition.Brain Sci.13:1347. 10.3390/brainsci13091347
11
BrissonC. D.LukewichM. K.AndrewR. D. (2013). A distinct boundary between the higher brain’s susceptibility to ischemia and the lower brain’s resistance.PLoS One8:e79589. 10.1371/journal.pone.0079589
12
BroselS.GrotheB.KunzL. (2018). An auditory brainstem nucleus as a model system for neuronal metabolic demands.Eur. J. Neurosci.47222–235. 10.1111/ejn.13789
13
Bustamante-BarrientosF. A.Luque-CamposN.ArayaM. J.Lara-BarbaE.de SolminihacJ.PradenasC.et al (2023). Mitochondrial dysfunction in neurodegenerative disorders: Potential therapeutic application of mitochondrial transfer to central nervous system-residing cells.J. Trans. Med.21:613. 10.1186/s12967-023-04493-w
14
CallegariF.BrofigaM.MassobrioP. (2023). Modeling the three-dimensional connectivity of in vitro cortical ensembles coupled to micro-electrode arrays.PLoS Comp. Biol.19:e1010825. 10.1371/journal.pcbi.1010825
15
ChenY. W.ZhangL. F.HuangJ. P. (2007). The Watts–Strogatz network model developed by including degree distribution: Theory and computer simulation.J. Phys. Math. Theoret.408237–8246. 10.1088/1751-8113/40/29/003
16
ChiappaloneM.PasqualeV.FregaM. (2019). “Preface: In vitro neuronal networks,” in In vitro neuronal networks: From culturing methods to neuro-technological applications. advances in neurobiology, edsChiappaloneM.PasqualeV.FregaM. (Berlin: Springer). 10.1111/jnc.16121
17
CompteA.Sanchez-VivesM. V.McCormickD. A.WangX.-J. (2003). Cellular and network mechanisms of slow oscillatory activity (<1 Hz) and wave propagations in a cortical network model.J. Neurophysiol.892707–2725. 10.1152/jn.00845.2002
18
de Santos-SierraD.Sendiña-NadalI.LeyvaI.AlmendralJ. A.AnavaS.AyaliA.et al (2014). Emergence of small-world anatomical networks in self-organizing clustered neuronal cultures.PLoS One9:e85828. 10.1371/journal.pone.0085828
19
Di FlorioM.IyerV.RajhansA.BuccelliS.ChiappaloneM. (2022). “Model-based online implementation of spike detection algorithms for neuroengineering applications,” in Proceedings of the 2022 44th annual international conference of the IEEE Engineering in Medicine & Biology Society (EMBC), (Piscataway, NJ: IEEE). 10.1109/EMBC48229.2022.9871444
20
DoornN.van HugteE. J. H.CiptasariU.MordeltA.MeijerH. G. E.SchubertD.et al (2023). An in silico and in vitro human neuronal network model reveals cellular mechanisms beyond NaV1. 1 underlying Dravet syndrome.Stem Cell Rep.181686–1700. 10.1016/j.stemcr.2023.06.003
21
DownesJ. H.HammondM. W.XydasD.SpencerM. C.BecerraV. M.WarwickK.et al (2012). Emergence of a small-world functional network in cultured neurons.PLoS Comp. Biol.8:e1002522. 10.1371/JOURNAL.PCBI.1002522
22
Emre KapucuF.VinogradovA.HyvärinenT.Ylä-OutinenL.NarkilahtiS. (2022). Comparative microelectrode array data of the functional development of hPSC-derived and rat neuronal networks.Sci. Data9:120. 10.1038/s41597-022-01242-4
23
Faria-PereiraA.MoraisV. A. (2022). Synapses: The Brain’s energy-demanding sites.Int. J. Mol. Sci.23:3627. 10.3390/ijms23073627
24
FiskumV.SandvigA.SandvigI. (2021). Silencing of activity during hypoxia improves functional outcomes in motor neuron networks in vitro.Front. Integrat. Neurosci.15:792863. 10.3389/fnint.2021.792863
25
GewaltigM.-O.DiesmannM. (2007). NEST (NEural Simulation Tool).Scholarpedia2:1430. 10.4249/SCHOLARPEDIA.1430
26
GhaderiP.MaratebH. R.SafariM.-S. (2018). Electrophysiological profiling of neocortical neural subtypes: A semi-supervised method applied to in vivo whole-cell patch-clamp data.Front. Neurosci.12:823. 10.3389/fnins.2018.00823
27
GordonJ.AminiS. (2021). General overview of neuronal cell culture.Methods Mol. Biol.23111–8. 10.1007/978-1-0716-1437-2_1
28
HinesM. L.CarnevaleN. T. (2001). Neuron: A tool for neuroscientists.Neuroscientist7123–135. 10.1177/107385840100700207
29
HodgkinA. L.HuxleyA. F. (1952a). A quantitative description of membrane current and its application to conduction and excitation in nerve.J. Physiol.117500–544. 10.1113/JPHYSIOL.1952.SP004764
30
HodgkinA. L.HuxleyA. F. (1952b). Currents carried by sodium and potassium ions through the membrane of the giant axon of Loligo.J. Physiol.116449–472. 10.1113/JPHYSIOL.1952.SP004717
31
HodgkinA. L.HuxleyA. F.KatzB. (1952). Measurement of current- voltage relations in the membrane of the giant axon of Loligo.J. Physiol.116424–448. 10.1113/JPHYSIOL.1952.SP004716
32
HofmeijerJ.MulderA. T. B.FarinhaA. C.van PuttenM. J. A. M.le FeberJ. (2014). Mild hypoxia affects synaptic connectivity in cultured neuronal networks.Brain Res.1557180–189. 10.1016/j.brainres.2014.02.027
33
HuchzermeyerC.BerndtN.Holzhü TterH.-G.KannO. (2013). Oxygen consumption rates during three different neuronal activity states in the hippocampal CA3 network.J. Cereb. Blood Flow Metab.33263–271. 10.1038/jcbfm.2012.165
34
HumpelC. (2015). Organotypic brain slice cultures: A review.Neuroscience30586–98. 10.1016/J.NEUROSCIENCE.2015.07.086
35
HumphriesM. D.GurneyK. (2008). Network “small-world-ness”: A quantitative method for determining canonical network equivalence.PLoS One3:e0002051. 10.1371/journal.pone.0002051
36
HyvärinenT.HyysaloA.KapucuE.AarnosL.VinogradovA.EglenS. J.et al (2019). Functional characterization of human pluripotent stem cell-derived cortical networks differentiated on laminin-521 substrate: Comparison to rat cortical cultures.Sci. Rep.9:17125. 10.1038/s41598-019-53647-8
37
IzhikevichE. M. (2003). Simple model of spiking neurons.IEEE Trans. Neural Netw.141569–1572. 10.1109/TNN.2003.820440
38
Klein GunnewiekT. M.Van HugteE. J. H.FregaM.GuardiaG. S.ForemanK.PannemanD.et al (2020). m.3243A > G-Induced mitochondrial dysfunction impairs human neuronal development and reduces neuronal network activity and synchronicity.Cell Rep.31:107538. 10.1016/j.celrep.2020.107538
39
KuznetsovA. V. (2024). Effects of time-dependent adenosine triphosphate consumption caused by neuron firing on adenosine triphosphate concentrations in synaptic boutons containing and lacking a stationary mitochondrion.J. Biomechan. Eng.146:111002. 10.1115/1.4065743
40
LenkK.SatuvuoriE.LallouetteJ.Ladrón-de-GuevaraA.BerryH.HyttinenJ. A. K. (2020). A computational model of interactions between neuronal and astrocytic networks: The role of astrocytes in the stability of the neuronal firing rate.Front. Comp. Neurosci.13:92. 10.3389/fncom.2019.00092
41
LennieP. (2003). The cost of cortical computation.Curr. Biol.13493–497. 10.1016/S0960-9822(03)00135-0
42
LonardoniD.Di MarcoS.AminH.BerdondiniL.NieusT. (2015). A computational model of cell culture dynamics: The role of connectivity and synaptic receptors in the appearance of synchronized bursting events.BMC Neurosci.16:177. 10.1186/1471-2202-16-S1-P177
43
MagliaroC.RinaldoA.AhluwaliaA. (2019). Allometric scaling of physiologically-relevant organoids.Sci. Rep.9:11890. 10.1038/s41598-019-48347-2
44
MarkramH.MullerE.RamaswamyS.ReimannM. W.AbdellahM.SanchezC. A.et al (2015). Reconstruction and simulation of neocortical microcircuitry.Cell163456–492. 10.1016/j.cell.2015.09.029
45
MasoliS.RizzaM. F.TognolinaM.PrestoriF.D’AngeloE. (2022). Computational models of neurotransmission at cerebellar synapses unveil the impact on network computation.Front. Comp. Neurosci.16:1006989. 10.3389/fncom.2022.1006989
46
MasquelierT.DecoG. (2013). Network bursting dynamics in excitatory cortical neuron cultures results from the combination of different adaptive mechanism.PLoS One8:e75824. 10.1371/journal.pone.0075824
47
McMurtreyR. J. (2016). Analytic models of oxygen and nutrient diffusion, metabolism dynamics, and architecture optimization in three-dimensional tissue constructs with applications and insights in cerebral organoids.Tissue Eng. Part C Methods22221–249. 10.1089/ten.TEC.2015.0375
48
NegriJ.MenonV.Young-PearseT. L. (2020). Assessment of spontaneous neuronal activity in vitro using multi-well multi-electrode arrays: Implications for assay development.Eneuro7:ENEURO.0080-19.2019. 10.1523/ENEURO.0080-19.2019.
49
NewcombJ. M.ToddK.BuhlE. (2023). Editorial: Invertebrate neurophysiology—of currents, cells, and circuits.Front. Neurosci.17:1303574. 10.3389/fnins.2023.1303574
50
NieberK. (1999). Hypoxia and neuronal function under in vitro conditions.Pharmacol. Therapeut.8271–86. 10.1016/S0163-7258(98)00061-8
51
ÖzugurS.KunzL.StrakaH. (2020). Relationship between oxygen consumption and neuronal activity in a defined neural circuit.BMC Biol.18:76. 10.1186/s12915-020-00811-6
52
PacittiD.PrivolizziR.BaxB. E. (2019). Organs to cells and cells to organoids: The evolution of in vitro central nervous system modelling.Front. Cell. Neurosci.13:129. 10.3389/fncel.2019.00129
53
PattersonM. S.MazurekE. (2010). Calculation of cellular oxygen concentration for photodynamic therapy in vitro.Methods Mol. Biol.635195–205. 10.1007/978-1-60761-697-9_14
54
Pires MonteiroS.VoogdE.MuzziL.De VecchisG.MossinkB.LeversM.et al (2021). Neuroprotective effect of hypoxic preconditioning and neuronal activation in a in vitro human model of the ischemic penumbra.J. Neural Eng.18:036016. 10.1088/1741-2552/abe68a
55
PlaceT. L.DomannF. E.CaseA. J. (2017). Limitations of oxygen delivery to cells in culture: An underappreciated problem in basic and translational research.Free Radical Biol. Med.113311–322. 10.1016/J.FREERADBIOMED.2017.10.003
56
PoliD.MagliaroC.AhluwaliaA. (2019). Experimental and computational methods for the study of cerebral organoids: A review.Front. Neurosci.13:162. 10.3389/fnins.2019.00162
57
PoliD.PastoreV. P.MassobrioP. (2015). Functional connectivity in in vitro neuronal assemblies.Front. Neural Circuits9:57. 10.3389/fncir.2015.00057
58
PotjansT. C.DiesmannM. (2014). The cell-type specific cortical microcircuit: Relating structure and activity in a full-scale spiking network model.Cereb. Cortex24785–806. 10.1093/cercor/bhs358
59
RadenkovicS.BudhrajaR.Klein-GunnewiekT.KingA. T.BhatiaT. N.LigezkaA. N.et al (2024). Neural and metabolic dysregulation in PMM2-deficient human in vitro neural models.Cell Rep.43:113883. 10.1016/j.celrep.2024.113883
60
RothA.van RossumM. C. W. (2009). Modeling synapses.Comp. Model. Methods Neuroscient.6, 139–160. 10.7551/mitpress/7543.003.0008
61
Sanchez-VivesM. V.NowakL. G.McCormickD. A. (2000). Cellular mechanisms of long-lasting adaptation in visual cortical neurons in vitro.J. Neurosci.204286–4299. 10.1523/JNEUROSCI.20-11-04286.2000
62
SantiagoJ.KreutzerJ.BossinkE.KallioP.le FeberJ. (2023). Oxygen gradient generator to improve in vitro modeling of ischemic stroke.Front. Neurosci.17:1110083. 10.3389/fnins.2023.1110083
63
SattelleD. B.BuckinghamS. D. (2006). Invertebrate studies and their ongoing contributions to neuroscience.Invertebrate Neurosci.61–3. 10.1007/s10158-005-0014-7
64
ScelfoB.PolitiM.RenieroF.PalosaariT.WhelanM.ZaldívarJ.-M. (2012). Application of multielectrode array (MEA) chips for the evaluation of mixtures neurotoxicity.Toxicology299172–183. 10.1016/j.tox.2012.05.020
65
ShefiO.GoldingI.SegevR.Ben-JacobE.AyaliA. (2002). Morphological characterization of in vitro neuronal networks.Phys. Rev. E66:021905. 10.1103/PhysRevE.66.021905
66
SpongK. E.AndrewR. D.RobertsonR. M. (2016). Mechanisms of spreading depolarization in vertebrate and insect central nervous systems.J. Neurophysiol.1161117–1127. 10.1152/jn.00352.2016
67
StimbergM.BretteR.GoodmanD. F. M. (2019). Brian 2, an intuitive and efficient neural simulator.eLife8:e47314. 10.7554/ELIFE.47314
68
SukenikN.VinogradovO.WeinrebE.SegalM.LevinaA.MosesE. (2021). Neuronal circuits overcome imbalance in excitation and inhibition by adjusting connection numbers.Proc. Natl. Acad. Sci. U. S. A.118:e2018459118. 10.1073/pnas.2018459118
69
Van PeltJ.VajdaI.WoltersP. S.CornerM. A.RamakersG. J. A. (2005). Dynamics and plasticity in developing neuronal networks in vitro.Prog. Brain Res.147171–188. 10.1016/S0079-6123(04)47013-7
70
VelascoV.ShariatiS. A.EsfandyarpourR. (2020). Microtechnology-based methods for organoid models.Microsyst. Nanoeng.6:76. 10.1038/s41378-020-00185-3
71
VoogdE. J. H. F.FregaM.HofmeijerJ. (2023). Neuronal responses to ischemia: Scoping review of insights from human-derived in vitro models.Cell. Mol. Neurobiol.433137–3160. 10.1007/s10571-023-01368-y
72
WalshK.MegyesiJ.HammondR. (2005). Human central nervous system tissue culture: A historical review and examination of recent advances.Neurobiol. Dis.182–18. 10.1016/j.nbd.2004.09.002
73
WangQ.PercM.DuanZ.ChenG. (2010). Impact of delays and rewiring on the dynamics of small-world neuronal networks with two types of coupling.Phys. Statist. Mechan. Appl.3893299–3306. 10.1016/J.PHYSA.2010.03.031
74
WattsD. J.StrogatzS. H. (1998). Collective dynamics of ‘small-world’ networks.Nature393440–442. 10.1038/30918
75
WattsM. E.PocockR.ClaudianosC. (2018). Brain energy and oxygen metabolism: Emerging role in normal function and disease.Front. Mol. Neurosci.11:216. 10.3389/fnmol.2018.00216
76
WeiY.UllahG.IngramJ.SchiffS. J. (2014). Oxygen and seizure dynamics: II. Computational modeling.J. Neurophysiol.112213–223. 10.1152/jn.00541.2013
77
WenJ.PeitzM.BrüstleO. (2022). A defined human-specific platform for modeling neuronal network stimulation in vitro and in silico.J. Neurosci. Methods373:109562. 10.1016/j.jneumeth.2022.109562
78
WilsonS. B.EmersonR. (2002). Spike detection: A review and comparison of algorithms.Clin. Neurophysiol.1131873–1881. 10.1016/S1388-2457(02)00297-3
79
ZaitsevA. V.PovyshevaN. V.Gonzalez-BurgosG.LewisD. A. (2012). Electrophysiological classes of layer 2/3 pyramidal cells in monkey prefrontal cortex.J. Neurophysiol.108595–609. 10.1152/jn.00859.2011
80
ZanelliS. A.RajasekaranK.GrosenbaughD. K.KapurJ. (2015). Increased excitability and excitatory synaptic transmission during in vitro ischemia in the neonatal mouse hippocampus.Neuroscience310279–289. 10.1016/J.NEUROSCIENCE.2015.09.046
81
ZhuJ.AjaS.KimE.-K.ParkM. J.RamamurthyS.JiaJ.et al (2012). Physiological oxygen level is critical for modeling neuronal metabolism in vitro.J. Neurosci. Res.90422–434. 10.1002/jnr.22765
Summary
Keywords
in silico modeling, neuron firing, oxygen metabolism, in vitro neuronal network, digitalized neuronal network
Citation
Fabbri R, Botte E, Ahluwalia A and Magliaro C (2025) Digitoids: a novel computational platform for mimicking oxygen-dependent firing of neurons in vitro. Front. Neuroinform. 19:1549916. doi: 10.3389/fninf.2025.1549916
Received
22 December 2024
Accepted
09 June 2025
Published
01 July 2025
Volume
19 - 2025
Edited by
Alberto Marchisio, New York University Abu Dhabi, United Arab Emirates
Reviewed by
Leonardo Dalla Porta, August Pi i Sunyer Biomedical Research Institute (IDIBAPS), Spain
Nina Doorn, University of Twente, Netherlands
Updates
Copyright
© 2025 Fabbri, Botte, Ahluwalia and Magliaro.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Chiara Magliaro, chiara.magliaro@unipi.it
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.