Neuronal and Astrocytic Regulations in Schizophrenia: A Computational Modelling Study

According to the tripartite synapse model, astrocytes have a modulatory effect on neuronal signal transmission. More recently, astrocyte malfunction has been associated with psychiatric diseases such as schizophrenia. Several hypotheses have been proposed on the pathological mechanisms of astrocytes in schizophrenia. For example, post-mortem examinations have revealed a reduced astrocytic density in patients with schizophrenia. Another hypothesis suggests that disease symptoms are linked to an abnormality of glutamate transmission, which is also regulated by astrocytes (glutamate hypothesis of schizophrenia). Electrophysiological findings indicate a dispute over whether the disorder causes an increase or a decrease in neuronal and astrocytic activity. Moreover, there is no consensus as to which molecular pathways and network mechanisms are altered in schizophrenia. Computational models can aid the process in finding the underlying pathological malfunctions. The effect of astrocytes on the activity of neuron-astrocyte networks has been analysed with computational models. These can reproduce experimentally observed phenomena, such as astrocytic modulation of spike and burst signalling in neuron-astrocyte networks. Using an established computational neuron-astrocyte network model, we simulate experimental data of healthy and pathological networks by using different neuronal and astrocytic parameter configurations. In our simulations, the reduction of neuronal or astrocytic cell densities yields decreased glutamate levels and a statistically significant reduction in the network activity. Amplifications of the astrocytic ATP release toward postsynaptic terminals also reduced the network activity and resulted in temporarily increased glutamate levels. In contrast, reducing either the glutamate release or re-uptake in astrocytes resulted in higher network activities. Similarly, an increase in synaptic weights of excitatory or inhibitory neurons raises the excitability of individual cells and elevates the activation level of the network. To conclude, our simulations suggest that the impairment of both neurons and astrocytes disturbs the neuronal network activity in schizophrenia.

According to the tripartite synapse model, astrocytes have a modulatory effect on neuronal signal transmission. More recently, astrocyte malfunction has been associated with psychiatric diseases such as schizophrenia. Several hypotheses have been proposed on the pathological mechanisms of astrocytes in schizophrenia. For example, post-mortem examinations have revealed a reduced astrocytic density in patients with schizophrenia. Another hypothesis suggests that disease symptoms are linked to an abnormality of glutamate transmission, which is also regulated by astrocytes (glutamate hypothesis of schizophrenia). Electrophysiological findings indicate a dispute over whether the disorder causes an increase or a decrease in neuronal and astrocytic activity. Moreover, there is no consensus as to which molecular pathways and network mechanisms are altered in schizophrenia. Computational models can aid the process in finding the underlying pathological malfunctions. The effect of astrocytes on the activity of neuron-astrocyte networks has been analysed with computational models. These can reproduce experimentally observed phenomena, such as astrocytic modulation of spike and burst signalling in neuron-astrocyte networks. Using an established computational neuron-astrocyte network model, we simulate experimental data of healthy and pathological networks by using different neuronal and astrocytic parameter configurations. In our simulations, the reduction of neuronal or astrocytic cell densities yields decreased glutamate levels and a statistically significant reduction in the network activity. Amplifications of the astrocytic ATP release toward postsynaptic terminals also reduced the network activity and resulted in temporarily increased glutamate levels. In contrast, reducing either the glutamate release or re-uptake in astrocytes resulted in higher network activities. Similarly, an increase in synaptic weights of excitatory or inhibitory neurons raises the excitability of individual cells and elevates the activation level of the network. To conclude, our simulations suggest that the impairment of both neurons and astrocytes disturbs the neuronal network activity in schizophrenia.

INTRODUCTION
Astrocytes contribute to the complex cognitive function in humans (Santello et al., 2019), whereby they have structural, homeostatic, and metabolic roles. Astrocytes are connected to each other via gap junctions and form non-overlapping domains in complex networks. This glial cell type composes, together with a pre-and a postsynaptic neuron, the so-called tripartite synapse (Araque et al., 1999). They take up glutamate rapidly from the synaptic cleft to ensure a short glutamate exposition to the postsynaptic neuron for precise synaptic transmission. Subsequently, astrocytes respond to neuronal activity with intracellular calcium (Ca 2+ ) elevations and release ions and transmitter molecules in return.
Seemingly neurons and astrocytes play a role in psychiatric diseases such as schizophrenia (SCZ) (Uhlhaas and Singer, 2010;Takahashi and Sakurai, 2013;Moraga-Amaro et al., 2014;Mei et al., 2018). SCZ is associated with positive (such as delusions and hallucinations), negative (such as social withdrawal and lack of motivation), and cognitive symptoms (difficulties in memory and attention).
Due to the success of first generation antipsychotics, the dopamine hypothesis of SCZ was formed. The underlying assumption of this hypothesis is that excessive transmission of subcortical dopamine causes the positive symptoms of SCZ (Grnder and Cumming, 2016). While research has found significant support of this hypothesis, it is not pathognomonic and does not sufficiently account for the negative symptoms of SCZ.
Post-mortem examinations, human induced pluripotent stem cell experiments, and animal studies have revealed a reduced astrocytic and neuronal synaptic density in SCZ patients (Cotter et al., 2001;Brennand et al., 2011;Dietz et al., 2020;Naujock et al., 2020). Cell death is associated with excitotoxicity due to excessive glutamate amounts and glutamate receptor activation (Drouin-Ouellet et al., 2011). However, the mechanisms and pathways of this developmental brain disorder are not well understood yet. Different studies reach conflicting conclusions. Some papers report an upregulation in glial fibrillary acidic protein (GFAP) in SCZ astrocytes (review by Moraga-Amaro et al., 2014), whereas in other studies, a reduced GFAP expression has been described (review by Takahashi and Sakurai, 2013). Koskuvi et al. (2020) found that human stem cell-derived astrocytes from SCZ patients show highly sex-specific alterations in gene expression. The Gene Ontology pathways related to neuronal wiring and inflammation are altered in SCZ astrocytes. Healthy neurons co-cultured with SCZ astrocytes of either sex exhibited a significantly increased response to glutamate. Transplanted hiPSC-astrocyte progenitors from SCZ patients developed into matured astrocytes in a mouse model and led to subtle behavioural changes. Furthermore, they induced demyelination, synaptic dysfunction and affected the inflammation pathways (Koskuvi et al., 2020). Several studies have shown that the expression of the metabotropic glutamate receptor (mGluR) type 5 is increased in astrocytes in the prefrontal cortex (Ohnuma et al., 1998;Kosten et al., 2016).
For neurons, the glutamate hypothesis of SCZ has been established (de Bartolomeis et al., 2005;Moghaddam and Javitt, 2012;Takahashi and Sakurai, 2013;Cohen et al., 2015;Mei et al., 2018): a hypofunction of NMDA receptors (NMDARs) leads to excessive glutamate release, mainly prevalent in the hippocampus and prefrontal cortex. However, the underlying mechanisms are not yet clarified. A decreased level of D-serine, which is mainly released by astrocytes, affects the NMDAR function in neurons. Both glutamate and either D-serine or glycine are coagonists to the NMDAR binding sites. The astrocyte-derived NDMAR antagonist kynurenic acid (KYNA) shows increased levels in SCZ patients (Takahashi and Sakurai, 2013;Cohen et al., 2015). Interneurons might react to changes in NMDAR function with inhibition, leading to disinhibition and a hyperexcitation of glutamate levels in SCZ (Moghaddam and Javitt, 2012;Mei et al., 2018). Alternatively, Schobel et al. (2013) suggested that excess synaptic glutamate release results in a reduction of the interneuronal function and hippocampal disinhibition. The abnormality of glutamate transmission is also related to astrocytes since they regulate the glutamate uptake and release from and to the extracellular space (Danbolt, 2001;Kondziella et al., 2007). While the glutamate hypothesis is more recent and less well studied than the dopamine hypothesis, it has the advantage of being able to account for a wider range of SCZ symptoms (Hu et al., 2015;Mei et al., 2018).
A third hypothesis is that, instead, the excessive glutamate release leads to the NMDAR hypofunction. Higher glutamate levels may result in a rise of astrocytic Ca 2+ levels. Subsequently, astrocytes may release higher concentrations of the gliotransmitter adenosine triphosphate (ATP). Lalo et al. (2016) suggested that the activation of postsynaptic purinergic P2X receptors (P2XRs) by ATP downregulates NMDARs.
These findings indicate that different neuronal and astrocytic malfunctions could lead to the same disease, despite different underlying mechanisms. This poses a major challenge to effective treatments for psychiatric diseases (Herman et al., 2012;Stephan et al., 2016). Computational models can complement or supersede experimental studies (Manninen et al., 2018;Oschmann et al., 2018), where the latter are not feasible or too expensive. They can guide experimental researchers to the relevant pathways that are altered in disease and thus predict pharmaceutical responses (Moran et al., 2016). Siekmeier and Hoffman (2002) proposed a computational model in which they reduced the number of synapses in a neuronal network. This synaptic pruning leads to semantic priming and difficulties in accessing memories in SCZ patients. Mahdavi et al. (2014) developed a tripartite synapse model focusing on the dysregulation of D-serine released by the astrocyte. To our best knowledge, no biophysical neuron-astrocyte network models have been published that explore the role of astrocytes in SCZ.
In this study, we investigate how changes in mechanisms and topology in neuron-astrocyte networks alter the functional state, i.e., spike and burst rate of neurons as well as glutamate levels in astrocytes, in health and SCZ. For this, we use the previously published computational model called INEXA  to test the following four SCZ-related hypotheses: we alter (1) the number of neurons or astrocytes in the network, FIGURE 1 | Schematic of the tripartite synapse, which consists of a pre-and a postsynapse and an astrocyte, modelled with INEXA . The presynapse releases glutamate with release rate f . The glutamate is taken up by the astrocyte with uptake rate g .
(2) the effect of astrocytic ATP on the postsynaptic activity, (3) the release of glutamate from the presynapse and the uptake of glutamate by the astrocyte, and (4) the excitatory and/or inhibitory synaptic strength.

MATERIALS AND METHODS
To study the three aspects mentioned above in SCZ, we used the published model INEXA  with the main components as shown in Figure 1. In section 2.1, we describe the governing equations of the INEXA model. In section 2.2, we explain how the model is used to simulate neuron-astrocyte networks in healthy and SCZ states. In sections 2.3, 2.4, the statistical analysis of the resulting neuronal and astrocytic activity is described.

Governing Equations of the INEXA Model
The phenomenological model INEXA  includes both excitatory and inhibitory neurons and was initially developed to study the astrocytes' effect on neuronal firing rates in the cortex. The governing equation for simulating neurons is the one for the firing rate λ i of a postsynaptic neuron i for each time step t k of 5 ms: where c i is the stochastic noise of neuron i. The term y ij denotes the synaptic weight between the presynaptic neuron j and a postsynaptic neuron i, which can vary at each time step. The synapse can be either excitatory (y ij ∈ [0, 1], representing glutamatergic neurons) or inhibitory (y ij ∈ [−1, 0], representing GABAergic interneurons). The parameter s j is a binary term indicating whether neuron j transmitted a spike in the previous time step t k−1 . The second term in the equation denotes the depressing effect caused by astrocytes. The variable A ija is binary, indicating whether a synapse ij is enclosed by astrocyte a and if astrocyte a was active at the previous time step. If this is the case, the astrocyte applies a depressing effect, y Astro , on the synapse. Otherwise, A ija is zero, removing the astrocytic effect on neuronal activity.
Due to limited experimental data on the interaction between GABAergic neurons and astrocytes, the INEXA model currently only connects astrocytes with glutamatergic synapses. These excitatory synapses release glutamate into the synaptic cleft with rate f . The released glutamate binds to astrocytic mGluRs with rate g causing the inositol 1,4,5-trisphosphate (IP 3 )-mediated release of Ca 2+ from the endoplasmic reticulum (ER) into the astrocytic cytoplasm (Figure 1).
The governing equation for astrocyte dynamics describes the intracellular calcium dynamics [Ca 2+ ] ija in an astrocyte a that encloses the synapse ij: The calcium concentration is composed of the calcium concentration left from the last time step (Ca 2+ ija (t k−1 )), the IP 3mediated Ca 2+ -induced Ca 2+ -release from the ER stores, and the Ca 2+ uptake back to the ER by the SERCA (sarco/endoplasmic reticulum Ca 2+ -ATPase) pumps. To reflect the slow dynamics of the calcium release (which can take up to multiple seconds), we

Simulated Experimental Setups
The topology and connectivity of the neuron-astrocyte network models were generated as described in Lenk et al. (2020) (see parameters in Table 1). Hence, 107 astrocytes and 250 neurons were placed uniformly at random on a planar area of 750 × 750 × 10 µm 3 , resembling the electrode field of an in vitro multielectrode array (MEA). Neurons and astrocytes within a distance of 30 µm were reallocated to account for the cell sizes. Astrocyte-astrocyte connections were formed whenever the distance was smaller than 100 µm. Neuron-neuron and astrocyte-neuron connections were created according to a truncated Gaussian distribution based on their distance.
To study the effects of various spatial network topologies and network connectivities, we generated three spatially different network topologies A, B, and C. For each network A, B, C, we computed three independent connectivity configurations (CC 1, 2, 3) according to the rules described above, resulting in a total of nine different networks. We considered those networks to represent the healthy, non-pathological baseline. Different topological properties of these initial networks are described in Table 2. The stochastic noise c i was set to 0.02 in all simulations.
To simulate healthy, non-pathological neuron-astrocyte network activity, we tuned the parameters of the INEXA model so that they resemble in vitro MEA recordings. Thereby, we used a reference spike rate of 30-200 spikes per minute as measured from mouse neurons (Gramowski et al., 2006;Jenkinson et al., 2017) and human-derived neurons (Tukker et al., 2018;Kizner et al., 2019).
Based on hypotheses and experimental evidence about SCZ described in the Introduction, we modified four different sets of INEXA parameters and examined their effect on neuronal and astrocytic activity. More specifically, we simulated the following four pathological changes in the neuron-astrocyte interaction: 1. We stochastically removed 25% of astrocytes or 25% of neurons from premade network configurations to account for the reduced astrocytic and neuronal synaptic density in SCZ patients. 2. We amplified the depressing signal y Astro of astrocytes to study the effect of astrocytic ATP on the postsynaptic activity. 3. We scaled the two parameters f and g by the factors 0.5 and 0.25, which control glutamate release from presynapses and glutamate uptake by astrocytes, respectively (Figure 1). Simulations were performed for each combination of the distorted and unmodified parameter values f and g . 4. We increased the maximal synaptic weight Y + max and Y − max for excitatory and inhibitory neurons, which both result in a higher excitability of the cells.
The latter two manipulations allowed us to simulate the excessive release of glutamate in excitatory neurons and astrocytes as well as the dysfunctional interneurons assumed in the glutamate hypothesis of SCZ. A summary of the modified parameters can be found in Table 1. All other model parameters not listed in this table were used as described in the original paper by Lenk et al. (2020). Each simulation consisted of T = 5 min simulated time and was repeated 10 times to account for the statistical variability of our results.

Quantifying Neuronal and Astrocytic Activity Features
To analyze the neuronal activity, we calculated spike and burst describing features. Bursts are cascades of spikes (i.e., action potentials). The features were determined using a modified version of the cumulative moving average (CMA) algorithm (Kapucu et al., 2012;Välkki et al., 2017). Therewith, we calculated the mean spike rate in spikes per minute and the mean burst rate in bursts per minute across all neurons, respectively. Additionally, we quantified the astrocytic activity by calculating how often an astrocyte changes its state to 'active' during the simulated time T and denoted the count as mean number of activations per astrocyte. We will refer to these three features as response variables below. The results were depicted as box plots.

Testing the Effects of Different Neuronal and Astrocytic Dysfunctions
The objective of this section is twofold: we quantify the effects of the four simulated pathological changes in SCZ (see section 2.2) on the neuronal and astrocytic response variables and estimate how much the effects vary for different network configurations.
In the previously described simulated experimental setups, there were two sources of stochasticity. One came from the stochastic nature of the INEXA model, e.g., when creating the excitatory and inhibitory synaptic weights and in the update equations. The other was due to the fact that neuron-neuron and astrocyteneuron connections were formed stochastically in our setup.
To estimate the average random effect from both sources of stochasticity, we performed all simulations on the nine initial networks, as described in section 2.2. We analysed the data with linear mixed-effects (LME) models using the lme4-package (Bates et al., 2015) in R, version 3.6.3 (2020-02-29). In order to satisfy the LME model requirements, it was sometimes necessary to transform the response variable Z by some function ϕ : R → R. In the following subsections, we decomposed the variance of ϕ(Z) into intra-network terms , which capture the stochasticity from different network connectivity configurations, and the residual error variance σ 2 . We note that these are the variance terms for the transformed response ϕ(Z). If one is interested in prediction intervals for Z, these must first be computed for ϕ(Z) and then back-transformed.

Hypothesis 1: Astrocyte or Neuron Removal
In the first set of simulated experiments, we quantified the average effect of stochastically removing 25% of the astrocytes or neurons, respectively. In order to estimate the random effect of different network configurations, we made three identical copies of each of the nine networks described in section 2.2 and Table 2.
Then, we either randomly removed 25% of astrocytes or 25% of neurons including their connections to other cells, which resulted in 27 cell-reduced networks per cell type. For all of these networks, we ran 10 simulations and fit the resulting data using an LME model with a fixed effect for cell removal and a random effect for different network configurations. Using contrasts to distinguish the full network from reduced networks (Oehlert, 2010) (see Supplementary Material), we then determined the average reduction δ rem of the response variables resulting from cell removal in percent. Moreover, we decomposed the variance of the transformed response ϕ(Z) into an inter-network term σ 2 rem and the residual error variance σ 2 . For each response variable, we report the restricted maximum likelihood (REML) estimates, 95% confidence intervals (CIs) of δ rem as well as the p-value for the null-hypothesis that the mean response of reduced and default networks A, B and C coincides. We also present the REML estimates and the 95% CIs for the standard deviations σ rem /μ and σ/μ, which have been normalised by the meanμ of the transformed response ϕ(Z), when no cells were removed. A more detailed description of the reported quantities mentioned in this section and the LME model can be found in the Supplementary Material.

Hypothesis 2: Effect of Astrocytic ATP
In this part, we altered the impact of astrocytic ATP on the postsynaptic P2XR and NMDAR activity by modifying the parameter y Astro . Therefore, we first performed 10 simulations with y Astro ∈ {0.01, 0.025, 0.05, 0.075, 0.1} for each of the nine networks described in section 2.2. To analyze our simulated data, we used the LME model where the index i enumerates the nine network configurations and j ∈ {1, . . . , 10} indices the 10 independent repetitions. The variable µ is the average value of log(Z) if y Astro was equal to zero. The coefficients β a , β a 2 correspond to the average increase of log(Z) in dependence of y Astro and y 2 Astro , respectively. Random effects are captured in the vectors u i : = (u i,µ , u i,a ) for each i, which are independent with multivariate Gaussian distribution u i ∼ N (0, ) with mean 0 and diagonal covariance matrix = diag(σ 2 µ , σ 2 a ) ∈ R 2×2 . This means that we assumed the individual components to be independent. The stochasticity of the INEXA model and the residual error is captured in ε ij , which we assumed to be normally distributed with mean 0, equal variance σ 2 and independent of u i for all indices i. We reported REML estimates and CIs for µ, β a , β a 2 , and σ 2 .

Hypothesis 3: Glutamate Dynamics Parameters
Following the glutamate hypothesis of SCZ, we tested the effect of changing the parameters regulating the glutamate dynamics in neurons and astrocytes. The parameters w f and w g for the glutamate dynamics are defined as follows: running the LME model with parameter values w f and w g corresponds to the scaling of the presynaptic glutamate release rate by a factor w f (equivalent to a multiplication of f by w f ) and scaling the astrocyte glutamate uptake rate by a factor w g (equivalent to multiplying g by w g ) compared to the default rates of the INEXA model. For every response variable Z = Z(w f , w f ), we used the LME model where the index i enumerates the nine network configurations described in section 2.2 and j ∈ {1, . . . , 10} denotes 10 independent repetitions. The variable µ stands for the average value of the response Z if w f and w g were equal to zero.
The coefficients β f , β g , β f 2 , β g 2 , β fg ∈ R correspond to the average increase of Z in dependence of the parameters they are multiplied by. Random effects were captured in vectors u i : = (u i,µ , u i,f , u i,g , u i,fg ), which are independent with multivariate Gaussian distribution u i ∼ N (0, ) with mean 0 and covariance matrix . We assumed = diag(σ 2 µ , σ 2 f , σ 2 g , σ 2 fg ) ∈ R 4×4 was a diagonal matrix with positive entries, which means that the individual components were assumed to be independent. The random effects u i are specific to each network configuration i, however we assumed them to be identically distributed. The estimates for σ 2 µ , σ 2 f , σ 2 g , σ 2 fg quantify the variance of Z as a stochastic function of parameter values w f and w g , considering the stochasticity from random network connectivities. The stochasticity of the INEXA model as well as the residual error was captured in ε ij ∼ N (0, σ 2 ), which we assumed to be normally distributed with mean 0, equal variance σ 2 and independent of u i for all indices i.
We performed 10 simulations with the INEXA model for each of the nine networks described in section 2.2 and every combination of the parameter values w f , w g ∈ {1, 1/2, 1/4}. The corresponding values for f , g are listed in Table 1. Fitting this data in R, we reported the REML estimates and CIs for µ, β f , β g , β f 2 , β g 2 , β fg , , and σ 2 . We created contour plots for the response Z(w f , w g ) as a function of w f and w g using the FIGURE 2 | (A) Mean spike rate (top), mean burst rate (middle) and mean number of activations per astrocyte (bottom) resulting from the networks with reduced astrocyte or neuron density compared to the default networks A-C with connectivity configuration CC 1 (bold). (B) Average glutamate amount in astrocytes when randomly removing 25% of either astrocytes (blue lines) or neurons (brown lines) from a premade network configuration B, CC 1 (black line). estimates provided by the LME model fit and setting all stochastic terms equal to 0.

Hypothesis 4: Changes in the Excitatory and Inhibitory Synaptic Weights
In this set of simulated experiments, we tested the dependency of neuronal and astrocytic activity on the inhibitory and excitatory synaptic weight boundaries Y − max and Y + max (see section 2.1). In our statistical LME model, we considered first and second order fixed effects of synaptic weight boundaries on the response variables as well as a random effect term to account for the variation caused by different network configurations. For ease of notation, let y in and y ex denote the parameter values for Y − max and Y + max , respectively. For every response variable Z = Z(y in , y ex ), we assumed the statistical model where the index i enumerates the nine network configurations described in section 2.2 and j ∈ {1, . . . , 10} denotes 10 independent repetitions. The variable µ refers to the estimated mean of Z if y in and y ex were equal to zero. The coefficients β in , β ex , β in 2 , β ex 2 , β in,ex ∈ R correspond to the average increase of Z in dependence of the parameters they are multiplied by. Random effects were captured in the random variables γ i , which we assumed to be independent with zero-mean Gaussian distribution and equal variance σ 2 net . The variance term σ 2 net quantifies how much the response varies between different stochastically chosen network connectivities if all other parameters are held constant. The residual error was captured in ε ij , which we assumed to be normally distributed with mean 0, equal variance σ 2 and independent of γ i for all indices i.
We performed 10 simulations with the INEXA model for each of the nine networks described in section 2.2 and every combination of the parameter values for Y − max , Y + max listed in Table 1. Fitting this data in R, we reported the REML estimates and 95% CIs for µ, β in , β ex , β in 2 , β ex 2 , β in,ex , σ net and σ . Moreover, we created contour plots for the response Z(y in , y ex ) as a function of y in and y ex using the estimates provided by the LME model fit while setting all stochastic terms equal to 0.

Hypothesis 1: Astrocyte or Neuron Removal Leads to Reduction of the Network Activity
In our first set of simulations, we reduced the cell density of either astrocytes or neurons by 25%. The cell removal resulted in a general decrease of mean spike rate, mean burst rate and mean number of astrocyte activations (Figure 2A). Overall, the decline of the response variables was larger when the neurons were removed from the network. In Figure 2B, the average amount of glutamate prevalent in astrocytes is shown at each time step. While both forms of reduction caused a decrease in astrocytic glutamate, the effect was more pronounced in the case of neuron removal. The p-values for the null-hypothesis that the mean response of reduced and default networks coincide were smaller than 10 −10 for all response variables and both simulated experiments.
The results of the LME model fit are listed in Table 3 for the astrocyte removal and in Table 3 for the neuron removal. They show that cell removal leads to statistically significant activity reduction. The 95% CIs of δ rem are tighter for the neuron removal than for astrocyte removal, indicating that the neuron removal data was fit more accurately by the LME model. The results highlight stark differences between the inter-network standard deviations σ rem /μ of the transformed responses ϕ(Z). While in the astrocyte removal experiments, the inter-network variance of the logarithm of the mean spike rate is relatively low, the standard deviation term for the logarithm of the mean number of astrocyte activations is estimated to be at around 50%. This shows that the robustness to changing network configurations varies significantly for the different response variables considered. Similarly, the standard deviation σ/μ of the residual error varies considerably between different response variables, but is consistently smaller than the inter-network standard deviation σ rem /μ.

Hypothesis 2: Effect of Astrocytic ATP
For our second hypothesis, we increased the depressing effect of the astrocytic gliotransmitter ATP on excitatory postsynapses.
As can be seen in Figure 3A, the depressing effect results in a decreased mean spike rate. In networks A and C, the activity reduction can also be observed for the other two response variables, mean burst rate and astrocyte activity. The increase of the depression, y Astro , from 0.01 to 0.025 has the proportionally largest effect. In network B, the behaviour of the mean burst rate and the astrocyte activity follows no clear tendency. Figure 3B depicts the average amount of astrocytic glutamate ready for release in presence with an ATP release from the astrocyte. Compared to the baseline level, we observed temporary higher levels of glutamate when amplifying the depressing effect.
In Table 4, we listed the REML estimates and confidence intervals of the parameters of the LME model described in section 2.4.2. We observed a positive second order fixed effect (β a 2 ) and a negative first-order one (β a ) for all logarithmically transformed responses log(Z ij ). While we obtain tight confidence intervals for the baseline activityμ, the first and second order coefficientsβ a ,β a 2 have comparatively wider confidence intervals. The standard deviation termσ µ of log(Z ij ) associated with µ is considerably lower than the standard deviation term σ a associated withβ a , moreover the confidence intervals ofσ µ are narrower than those ofσ a , indicating that the first order term is difficult to fit by our model. The error standard deviation σ is comparatively low and has tight confidence intervals for all responses.

Hypothesis 3: Reducing the Glutamate
Release and Uptake Rate Results in Higher Network Activity Figure 4A depicts the behaviour of the response variables with respect to varying rates of synaptic facilitation, f , and varying recovery rates of the gliotransmitter receptors, g . In general, a reduction in both, g and f , resulted in an elevation of the response variables. However, the increase between two consecutive values for f was larger than between two values of g . An exception was observed for network B, for which the modification of g did not result in any considerable changes of the response variables.
We observed a similar behaviour for the average glutamate level when varying f and g . Figure 4B shows a larger effect of g compared to the effect of f . Decreasing either parameter resulted in elevated response variables compared to the default network. Also in this case, the effect of f was larger compared the effect of g .
The results of the LME model fit are reported in Table 5. We observed that the average response µ fitted for values w f = w g = 0 (baseline activity) was difficult to estimate, indicated by the wide confidence intervals. Also the estimated standard deviationσ µ of the baseline activity due to different network configurations was large relative to µ for all response variables, particularly for the mean number of astrocyte activations. For most response variables, the estimated standard deviationσ of the residual error was at around 10% of the baseline activity, significantly smaller thanσ µ .
The parameters β f , β f 2 had tight confidence intervals, indicating that the estimates were reliable. We saw a strong negative first order effect of the rate of presynaptic glutamate uptake w f on all responses and a positive second-order effect. For all responses, the standard deviation estimate σ f of the first order effect was small relative toβ f . Similarly, we found a negative first order and a positive second-order effect of the astrocyte glutamate uptake on all variables with small estimated standard deviationσ g . The corresponding 95% CIs for the parameters β g , β g 2 were slightly wider than those of β f and β f 2 . The firstorder interaction coefficient β fg of the parameters w f and w g was significant, but with a small effect for all responses. The fact thatσ fg and β fg lay in a similar range indicates that the interaction terms might vary considerably for different network configurations. However, these effects were all small compared to β f or β g . Contour plots of the response variables derived from the LME model fit are shown in Figure 5.

Hypothesis 4: Higher Network Activity When Increasing Synaptic Weight Parameters
In this section, we computationally investigated the neuronal and astrocytic response to synaptic perturbations caused by SCZ. Figure 6A depicts the behaviour of the three response variables while changing either the maximum inhibitory synaptic weight Y − max or in the maximum excitatory synaptic weight Y + max . Increases in Y − max and Y + max yielded in elevated response variables. The relationship between the parameter distortion and the effect on the observation variables were consistent among the three networks (A, B, C) and their corresponding connectivities (1, 2, 3). However, network B showed a comparatively reduced response to deviations in Y − max . The variance induced by alterations of Y + max was higher than the inter-network and the inter-connectivity variance. This was not the case for the variance induced by deviating Y − max , which was due the lower number of inhibitory cells in the networks. Moreover, compared to the unmodified parameter configuration (Y − max = −0.7 and Y + max = 0.7), the response variables from the distorted parameter configurations showed higher variances in most of the cases.
In Figure 6B, the average glutamate amount in all astrocytes is plotted against time for varying Y − max and Y + max . We noted that a change in the glutamate behaviour was achieved by all deviation levels of Y + max . For Y − max , a large deviation from the default values was required to see a change in the average glutamate levels. The smaller deviations of Y − max did not yield a notable elevation of neither the response variables ( Figure 6A) nor the time course of the average glutamate level (Figure 6B).
The parameter estimates of the corresponding statistical LME analysis are listed in Table 6. We refer to the average µ of the response variable when setting parameter values y in , y ex to zero as the baseline activity. Comparing the parameter estimates in Table 6 to the data in Figure 6A, we observe that the estimated baseline activityμ is significantly larger than the average response. This effect is compensated by large values of parameterŝ β in ,β ex ,β in 2 ,β ex 2 andβ in,ex . Thus, to analyse the magnitude of the standard deviation terms, we refer to Figure 6A. We observed that the estimated residual standard deviationσ was at around 1/10 to 1/2 of the average response value. The confidence intervals for σ net were wide, so the random effect due to different network connectivities was not clear. However, across all experiments,σ net was estimated to be smaller thanσ . Thus, it is likely that on average, responses vary more significantly for different repetitions than between networks with different connectivities. Further, we observed relatively wide confidence intervals for all fitted coefficients, indicating that the model fit and the provided estimates have limited expressiveness. Contour plots of the response variables derived from the LME model fit are shown in Figure 5.

DISCUSSION
Schizophrenia is characterised by diverse dysregulations and malfunctions in neurons and astrocytes (Takahashi and Sakurai, 2013;Moraga-Amaro et al., 2014;Mei et al., 2018). In astrocytes, the pathology is among others indicated by a higher GFAP expression, an altered morphology, loss of gap junctions between cells, and an elevated glutamate release (Mitterauer, 2009;Catts et al., 2014;Mei et al., 2018). In neurons, SCZ is indicated by morphological alterations, reductions in cell density, and excessive glutamate release, which is due to the hypofunction of NMDA receptors (Cotter et al., 2001;Mei et al., 2018). Spontaneous neuronal activity is a fine balance between excitation and inhibition (Lisman, 2012;Sohal and Rubenstein, 2019), which is between 30 and 200 spikes per minute in healthy human stem cell-derived neurons (Tukker et al., 2018;Kizner et al., 2019). However, it is disputed whether SCZ causes an increase or a decrease in neuronal and astrocytic activity (Mei et al., 2018;Naujock et al., 2020).
Computational models can support the understanding of biophysical pathways when experiments are not realisable or too expensive. So far, a multitude of biophysical (bottomup) and symptom focused (top-down) models have been developed. According to the reviews by Rolls et al. (2008) and Valton et al. (2017), top-down models mostly support the dopamine hypothesis. For example, Hoffman and McGlashan (2006) explored the effects of neuronal pruning and excessive dopamine. Using a recurrent backpropagation model of working memory, they showed how neuron reduction or hyperdopamigenic systems can produce spontaneous percepts simulating hallucinated speech. Furthermore, they were able provide an explanation for the onset of SCZ in late adolescence. For the glutamate hypothesis, only a few computational models are available (Valton et al., 2017). The models mostly focused on excitatory and inhibitory neurons, showing that deficits in glutamatergic or GABAergic activity lead to impairments in workingmemory and other cognitive processing (Valton et al., 2017).
In the last few years, attempts were made to simulate SCZ using the Bayesian inference hypothesis. This hypothesis assumes that the brain combines sensory evidence together with prior knowledge and expectation to interpret a stimulus.  Table 5, and as a function of the synaptic weight bound parameters y in , y ex (right column) using the parameter estimates from Table 6.
In SCZ, different delusions bias the expected value, leading to stronger perceptual biases in the future. So far, the hypothesis has mostly been tested for illusions (Valton et al., 2017). One of the first computational models in that area was made by Adams et al. (2013). Using a biological predictive coding scheme, they showed that psychosis in SCZ can be accounted by an increased prior precision and trait abnormalities (e.g., abnormalities of smooth pursuit eye movements, eventrelated brain potentials or anhedonia) by a decrease in prior precision.
The loss of gap junctions in astrocytes can be a hallmark in SCZ (Mitterauer, 2009). Using the INEXA model, Genocchi et al. (2020) studied the effect of gap junction uncoupling on the neuronal activity. Additionally, Lenk et al. (2021) showed that larger distances between astrocytes, and thus, less gap junctions, seem to centralise the information transmission in astrocytes. The spike rate and burst rate increased until a cell-cell distance of 100 µm and then decreased again. In this paper, we investigated the effect resulting from four different hypotheses of SCZ on neuronal and astrocytic activity.
Using the computational network model by Lenk et al. (2020), we separately examined the neuronal and astrocytic dynamics depending on (1) the neuron and astrocyte density in the network, (2) the depressing effect of astrocytes on excitatory neurons, (3) the glutamate release from the presynapse and the glutamate uptake by the adjacent astrocyte, and (4) the excitatory and/or inhibitory synaptic weight between neurons. Furthermore, we applied a set of statistical models to quantify the effect of these pathological changes.
Several studies report that the neuronal and astrocytic density decreases in several brain regions in SCZ (Cotter et al., 2001;Brennand et al., 2011;Williams et al., 2013), which may be caused by excessive glutamate levels (Drouin-Ouellet et al., 2011). Wagenaar et al. (2006) studied the effect of neuronal cell density on the spike and burst behaviour recorded with MEAs. They found that sparser cultures exhibit a lower spike rate as well as delayed network development. In our simulated experiments, we found that reducing the number of neurons or astrocytes by 25% (hypothesis (1)) causes a significant decrease in neuronal and astrocytic activity. In the original paper of our computational model , the number of astrocytes was increased, which led to an elevation of the neuronal activity. Conversely, this means that when astrocytes and their connections are eliminated, neuronal activity also decreases. Astrocytes are assumed to have a homeostatic effect on neurons, they can both enhance and dampen neuronal activity (Santello et al., 2019). Additionally, we found that the effect of reducing astrocytes or neurons can vary significantly for different network connectivity configurations. While the effect on the mean spike rate is robust to varying network connectivity configurations, the mean burst rate and the number of astrocyte activations are not. This might be useful for relating the behaviours of networks with the same number of cells, but different or unknown network configurations.
Furthermore, we have studied the effect of ATP release from the astrocyte on the postsynapse [hypothesis (2)]. Lalo et al. (2016) have demonstrated that ATP activates postsynaptic P2XRs resulting in a downregulation of neuronal NMDARs. This negative feedback mechanism might prevent cell damage due to excessive Ca 2+ and glutamate concentrations (Singh and Abraham, 2017). By the regulation of NMDARs, P2XRs might also be involved in synaptic plasticity and long-term potentiation (LTP) (Moraga-Amaro et al., 2014;Lalo et al., 2016;Singh and Abraham, 2017). This receptor type could even have a protective role by avoiding excitotoxicity and excessive LTP (Lalo et al., 2016).
An increase in neuronal and astrocytic activity was observed when altering the rate of presynaptic facilitation and the recovery rate of mGluRs [hypothesis (3)]. The decreased facilitation rate increases the amount of neurotransmitters ready for release to the synaptic cleft, thereby increasing the release probability and spike rate. Similarly, decreasing the recovery rate of gliotransmitter receptors results in a higher presynaptic potentiation (Helen et al., 1992;Agulhon et al., 2008). Subsequently, this leads to an increase of the neuronal activity.
According to the glutamate hypothesis of SCZ, NMDARs undergo a hypofunction leading to an excessive glutamate release (Mei et al., 2018). The changes in the NMDAR function may be caused by D-serine (Takahashi and Sakurai, 2013;Cohen et al., 2015;Mei et al., 2018), which we implicitly modelled by increasing the excitatory synaptic weights. In the future, we aim to add D-serine explicitly to the INEXA model for investigating its role in long-term potentiation and thus memory (Henneberger et al., 2010). Furthermore, GABAergic interneurons (Uhlhaas and Singer, 2010) might be impacted by SCZ due to reaction to the alterations in the NMDARs and an impaired synthesis and reuptake of GABA into the cell (Uhlhaas and Singer, 2010;Moghaddam and Javitt, 2012;Mei et al., 2018). When investigating hypothesis (4), the results show a higher neuronal and astrocytic activity, when increasing the inhibitory or excitatory synaptic weights (both yields in a neuronal excitation).
Further experimental observations are coherent with our modelling results. In SCZ, glutamine synthetase, which catalyzes glutamate, was shown to have reduced activity leading to higher levels of intracellular astrocytic glutamate, and thus, result in a neuronal hyperexcitability (Hu et al., 2015;Mei et al., 2018). This observation is in agreement with the results of our simulated experiments from hypotheses (3) and (4). The accumulation of glutamate was most pronounced in simulated experiments with a decreased mGluR recovery rate. In these experiments, glutamate accumulated quickly before stabilizing at an extremely high level. As to be expected, the neuronal activity increased with dysinhibition of the GABAergic neurons or the excitation of the glutamatergic neurons.
In the INEXA model, we mainly concentrated on the IP 3 -dependent Ca 2+ dynamics. Srinivasan et al. (2015) showed that somatic calcium is abolished in adult Ip3r2 −/− mice. In the future, we intend to extend the model by adding astrocytic NMDARs and α-amino-3-hydroxy-5methyl-4-isoxazolepropionic acid (AMPA) receptors as well as the glutamate transporter-dependent pathway to gain further insights into the possible dysfunctions in schizophrenia.
Dynamical causal modeling (DCM), a framework using Bayesian statistics to compare competing models (Friston et al., 2003), has given rise to differential diagnosis methods for psychiatric diseases (Stephan et al., 2017). In the future, we will apply Bayesian inference to obtain posterior distributions of certain astrocytic parameters from INEXA, building the basis for Bayesian model selection among competing disease mechanisms linked to astrocyte malfunction.

DATA AVAILABILITY STATEMENT
The datasets for this study can be found in Github: https://github. com/kerstinlenk/INEXA_SCZ/.

AUTHOR CONTRIBUTIONS
KL wrote the first draft of the manuscript. LF, JL, and FS collected simulation data. LF, FS, and KL visualized the results. FS designed and performed the statistical analysis. All authors designed and performed research and contributed to the manuscript writing and revision. In addition, they have read and approved the submitted version.