Mean-field analysis of synaptic alterations underlying deficient cortical gamma oscillations in schizophrenia

Deficient gamma oscillations in the prefrontal cortex (PFC) of individuals with schizophrenia (SZ) are proposed to arise from alterations in the excitatory drive to fast-spiking interneurons (E→I) and in the inhibitory drive from these interneurons to excitatory neurons (I→E). Consistent with this idea, prior postmortem studies showed lower levels of molecular and structural markers for the strength of E→I and I→E synapses and also greater variability in E→I synaptic strength in PFC of SZ. Moreover, simulating these alterations in a network of quadratic integrate-and-fire (QIF) neurons revealed a synergistic effect of their interactions on reducing gamma power. In this study, we aimed to investigate the dynamical nature of this synergistic interaction at macroscopic level by deriving a mean-field description of the QIF model network that consists of all-to-all connected excitatory neurons and fast-spiking interneurons. Through a series of numerical simulations and bifurcation analyses, findings from our mean-field model showed that the macroscopic dynamics of gamma oscillations are synergistically disrupted by the interactions among lower strength of E→I and I→E synapses and greater variability in E→I synaptic strength. Furthermore, the two-dimensional bifurcation analyses showed that this synergistic interaction is primarily driven by the shift in Hopf bifurcation due to lower E→I synaptic strength. Together, these simulations predict the nature of dynamical mechanisms by which multiple synaptic alterations interact to robustly reduce PFC gamma power in SZ, and highlight the utility of mean-field model to study macroscopic neural dynamics and their alterations in the illness.


Introduction
Deficits in certain cognitive processes, such as working memory, are among the core clinical features of schizophrenia (SZ) (Kahn and Keefe, 2013;McCutcheon et al., 2020).Working memory is associated with neural population activity oscillating at gamma frequency (30-100 hertz) in the prefrontal cortex (PFC) (Fries, 2009;Miller et al., 2018), and individuals with SZ exhibit lower PFC gamma power during working memory tasks (Chen et al., 2014;Cho et al., 2006;Minzenberg et al., 2010).PFC gamma oscillations are generated, in part, by a local circuit formed by regular-spiking excitatory neurons (RSEs) and fast-spiking GABAergic interneurons (FSIs) (Gonzalez-Burgos et al., 2015).In this circuit, recurrently connected local RSEs provide excitatory drive to FSIs (E→I), and FSIs in turn provide inhibitory drive back to RSEs (I→E).Following decay of phasic inhibition from FSIs, local RSEs fire synchronously to recruit FSIs, and the repeated cycles of this reciprocal synaptic activity are thought to produce gamma oscillations in the PFC (Whittington et al., 2011(Whittington et al., , 2000)).Thus, identifying neural substrates underlying lower PFC gamma power in SZ requires understanding the mechanisms by which E→I and I→E synapses regulate the generation of gamma oscillations.
Computational modeling of neural circuits has proven to be an invaluable tool for understanding dynamical mechanisms by which changes in synaptic parameters regulate the generation of gamma oscillations (Ermentrout and Chow, 2002;Jadi et al., 2016;Pittman-Polletta et al., 2015).For example, computational studies that utilized model networks of quadratic integrate-and-fire (QIF) neurons showed that key determinants of optimal gamma power include the strength of E→I and I→E synapses (Chung et al., 2022;Gonzalez-Burgos et al., 2015;Jadi and Sejnowski, 2014;Kirli et al., 2014;Rotaru et al., 2011;Spencer, 2009).In line with these predictions, prior postmortem human brain studies showed lower levels of molecular and structural markers for the strength of E→I (Chung et al., 2016a(Chung et al., ,b, 2018) ) and I→E (Chung et al., 2023 in press;Glausier and Lewis, 2011;Fish et al., 2021;Hashimoto et al., 2003) synapses in PFC of SZ, suggesting lower strength of these synapses in the illness.In addition, a recent study reported higher variance in the expression levels of protein markers for E→I synaptic strength across individual FSIs in PFC of SZ, and that increasing cell-to-cell variability in E→I synaptic strength disrupts gamma oscillations in a QIF model network (Chung et al., 2022).Furthermore, this study showed that SZassociated synaptic alterations (i.e., lower E→I and I→E synaptic strength and greater variability in E→I synaptic strength) non-linearly interact to synergistically (i.e., greater than the sum of individual effect) reduce gamma power.Together, these findings provide computational insights into how small changes in multiple synaptic parameters may result in robust deficits in PFC gamma power in SZ.However, dynamical mechanisms underlying the synergistic interactions of these synaptic alterations in PFC of SZ still remain elusive.
Recent studies showed that Ott-Antonsen ansatz (Ott and Antonsen, 2008) allows a derivation of mean-field description of QIF neurons in terms of their mean firing rate and mean membrane potential (Byrne et al., 2022;Devalle et al., 2017;Gast et al., 2020;Montbrió et al., 2015;Schmidt et al., 2018).These measures can then be mathematically examined by computing bifurcation diagrams, providing an analytical tool to systematically investigate how multiple synaptic parameters interact to regulate macroscopic dynamics of neural network (Bick et al., 2020).In this study, we utilized a mean-field model of QIF neurons to investigate the dynamical mechanisms by which SZ-associated synaptic alterations interact to synergistically reduce PFC gamma power.First, we derived a mean-field description of QIF model network that consists of all-to-all connected RSEs and FSIs via AMPA, NMDA and GABA synapses.We then validated our mean-field model by comparing its gamma oscillation dynamics to that of the QIF model network across a broad spectrum of synaptic parameter choices.Finally, we performed a series of numerical simulations and two-dimensional (2D) bifurcation analyses to investigate how multiple synaptic alterations between RSEs and FSIs found in SZ may interact to non-linearly affect the generation of gamma oscillations.

Description of a model network of spiking neurons
Our prior studies utilized a pyramidal interneuron gamma (PING) model network that generates robust oscillatory activity at gamma frequency to simulate how gamma power is regulated by various properties of excitatory and inhibitory synapses onto regular-spiking excitatory (RSEs) and fastspiking inhibitory (FSIs) neurons (Rotaru et al., 2011;Gonzalez-Burgos et al., 2015;Chung et al., 2022Chung et al., , 2023 in press) in press).In this study, we sought to derive a set of exact equations that describe the macroscopic dynamics of a PING model network.
First, we consider a PING model network that consists of RSE and FSI quadratic integrate-andfire (QIF) point neurons.The RSEs and FSIs are reciprocally connected in an all-to-all manner via AMPA, NMDA and GABA synapses.
The membrane potential, V (t), of each neuron is governed by the following QIF model: where C is the membrane capacitance, g l is the leak conductance, V l is the leak equilibrium potential, and V T is the threshold potential.In this model, neurons fire when V (t) reaches the spike voltage (V spike ).Once neurons spike, their V (t) returns to the reset voltage (V R ).Network activity is initiated by external tonic current (I appl ) that is applied to RSEs only.The synaptic currents (I syn ) are the sum of AMPA (g e ), NMDA (g n ) and GABA (g i ) synaptic conductance as described by the following equation: where V ex is the AMPA and NMDA equilibrium potential and V in is the GABA equilibrium potential.The dynamics of synaptic gating (the fraction of open channels) for AMPA (s e ), NMDA (s n ) and GABA (s i ) synapses are described by the following equations: where t k j is the kth spike of the jth presynaptic neuron and N is the total number of neurons.τ e , τ n and τ i are the decay time constants for AMPA, NMDA, and GABA synaptic conductance, respectively.We do not include the postsynaptic voltage effects on the NMDA channels.

Macroscopic description for the model network of spiking neurons
Following (Montbrió et al., 2015), we derive a set of equations that exactly describe the macroscopic dynamics of a network of all-to-all connected QIF neurons.First, we introduce a nondimensionalized variable u, where Hence, equation (1) becomes To rewrite equation (2), we introduce the following parameters: Substituting these parameters into equation (2), we have: Then, we rewrite equation (3) in a form of the following equation: where Here, we assume that g e , g n , g i and I appl conform to a Lorentzian distribution with the mode of ḡe , ḡn , ḡi and Ī, and the half-width of γ e , γ n , γ i and ∆, respectively.Subsequently, g e , g n , g i and I appl can be expressed by the following equations: where ξ e , ξ n , ξ i and η have the Lorentzian distribution, p(ξ) = 1/[π(1 + ξ 2 )].
We now let the size of the network go to infinity and set V spike = +∞ and V R = −∞, so that our QIF model is equivalent to a theta model and we can employ the methods of (Montbrió et al., 2015).Given the conservation of number of neurons in the network, we have the following continuity equation for the probability density function ρ(u|ξ e , ξ n , ξ i , η, t): where ρ(Au 2 + Bu + C) is the flux of membrane potentials of QIF neurons.Based on (Montbrió et al., 2015), we assume that the solution of equation ( 5) converges to a Lorentzian-shape function independent of the initial conditions.This ansatz is expressed by the following equation: (6) where β(ξ e , ξ n , ξ i , η, t)and α(ξ e , ξ n , ξ i , η, t) are the mode and the half-width, respectively, of the Lorentzian distribution.
We now abbreviate α(ξ e , ξ n , ξ i , η, t) as α and β(ξ e , ξ n , ξ i , η, t) as β, and let α ′ denote ∂α ∂t and β ′ denote ∂β ∂t .Substituting the ansatz into the continuity equation and equating powers of u, we find that: The firing rate of a neuron with parameters Ξ = (ξ i , ξ n , ξ i , η) ∈ R 4 is the flux at infinity: and the mean firing rate of the population of neurons is where P (Ξ) is the joint probability density for the four Lorentzian distributed parameters, Ξ.
The dynamics of synaptic gating for AMPA synapses are rewritten as: Similarly, the dynamics of synaptic gating for GABA synapses are rewritten as: As for the dynamics of synaptic gating for NMDA synapses, In order to obtain a closed system for the mean firing rate, mean potential, and synaptic gating variables, we need to evaluate the integral in equation ( 8).Since ξ e , ξ n , ξ i and η are Lorentzian distributed, Cauchy's residue theorem is used to perform integration as previously described (Montbrió et al., 2015), yielding: where σ e , σ n , σ i , σ ∈ {−1, +1} are to be determined by the signs of the terms in the ODEs, so that the firing rate does not fall below zero.
Similarly, the mean membrane potential of QIF neurons, β, with ξ e , ξ n , ξ i , and η at time t is expressed by the following equation: As ξ e , ξ n , ξ i and η are Lorentzian distributed, Cauchy's residue theorem is used to obtain the population membrane potential, u, so that where σ e , σ n , σ i and σ take the same values as those in equation ( 9).Substituting g e = ḡe +γ e σ e i, g n = ḡn + γ n σ n i, g i = ḡi + γ i σ i i, and I = Ī + ∆σi into equation ( 7), we have We must guarantee that the mean firing rate, r = α/π, is never negative and that dα dt > 0 when α = 0, which implies that: Therefore, ODEs for the population firing rate and membrane potential become, respectively, and Since the PING network consists of all-toall connected RSEs and FSIs, the ODEs for the population firing rate and membrane potentials can be expressed as below to construct a meanfield model of PING network; the parameters for RSEs and FSIs are denoted with (e) and (i), respectively: 3 Results To further compare gamma oscillation dynamics between the QIF network and mean-field model, we quantified changes in the gamma power and the peak gamma frequency in response to increasing external current (Figure .1B).Increasing the strength of external current to RSEs augmented gamma power to a similar extent between the QIF network and the mean-field model.Furthermore, increasing the external current led to identical increments in peak gamma frequency in both models.Together, these findings demonstrate the robustness of mean-field model in extending the gamma oscillation dynamics observed in the QIF model network to the macroscopic level.

Characterizing the effect of synaptic strength on gamma power in mean-field model
Multiple studies using in vivo models have shown that E→I synaptic strength is a critical determinant of cortical gamma power (Sohal et al., 2009;Carlen et al., 2012;Pelkey et al., 2015).To investigate how E→I synaptic strength regulates gamma oscillation dynamics in the mean-field model, we characterized the effect of mean AMPA-mediated (ḡ E→I ) and mean NMDA-mediated (ḡ N →I ) conductance in FSIs on gamma power.First, we numerically simulated the effect of ḡE→I on gamma power (Figure .2A).Gamma power sharply increased as ḡE→I increased from zero and reached a peak at ḡE→I =0.7.As ḡE→I continued to increase beyond this point, gamma power began to decline, eventually reaching zero at ḡE→I =5.These changes resulted in an inverted-U relationship between ḡE→I and gamma power, which mirrors the findings from a QIF model network in our prior study (Chung et al., 2022).The inverted-U relationship was also observed in the bifurcation diagram, which showed a Hopf bifurcation point at ḡE→I =0.26 marking the onset of gamma oscillations (Figure .2D).
Next, we assessed the effect of ḡN→I on gamma power (Figure .2B).Numerical simulation showed maximal gamma power at ḡN→I =0, which then gradually decreased as ḡN→I was increased to 0.12.Further increasing ḡN→I resulted in a sharp decrease in gamma power, which reached zero at ḡN→I =0.15.Consistent with these changes, a Hopf bifurcation point was observed at ḡN→I =0.15 corresponding to the offset of gamma oscillations (Figure .2E).These findings align with a prior finding in a QIF model network in which increasing NMDA-mediated conductance in FSIs reduced gamma power (Rotaru et al., 2011).
In addition to E→I synapses, I→E synaptic strength is proposed to be another critical determinant of cortical gamma power (Cho et al., 2015;Berryer et al., 2016).To investigate the relationship between I→E synaptic strength and gamma oscillation dynamics in the mean-field model, we characterized the effect of mean GABA-mediated conductance in RSEs (ḡ I→E ) on gamma power.
Numerical simulation showed that gamma power sharply increased as ḡI→E increased from zero and reached a peak at ḡI→E =0.68 (Figure .2C).Further increasing ḡI→E resulted in a sharp decrease in gamma power which reached zero at ḡI→E =1.52.In line with these findings, the bifurcation diagram showed a Hopf bifurcation point at ḡI→E = 1.52 corresponding to the offset of gamma oscillations (Figure .2F).Thus, similar to ḡE→I , ḡI→E also formed an inverted-U relationship with gamma power, consistent with its effect in a QIF model network in our prior study (Chung et al., 2022).

Characterizing the effect of synaptic variability on gamma power in mean-field model
Excitatory synaptic inputs to cortical neurons intrinsically vary in their strength (Ko et al., 2011;Hofer et al., 2011) and a recent study showed that increasing the variability in E→I synaptic strength monotonically reduces gamma power in a QIF model network (Chung et al., 2022).To characterize the relationship between the variability in E→I synaptic strength and gamma oscillation dynamics in the mean-field model, we assessed how gamma power changes in response to increasing the half-width of AMPA-mediated (γ E→I ) and NMDA-mediated (γ N →I ) conductance across FSIs.
Increasing γ E→I progressively reduced gamma power, which reached zero at γ E→I =1.55 (

Synergistic interactions of SZ-associated synaptic alterations in mean-field model
Above findings show that our mean-field model reliably captures changes in gamma oscillation dynamics in response to shifts in mean and variability of synaptic strengths at the macroscopic level.Prior postmortem human brain studies suggest that these synaptic parameters are altered in SZ (Table.1).For example, we previously reported fewer E→I synapses in PFC of SZ (Chung et al., 2016a).Also, multiple studies showed lower levels of the GABA-synthesizing enzyme glutamic acid decarboxylase 67 (GAD67), a marker for I→E synaptic strength, in PFC of SZ (Glausier and Lewis, 2011;Fish et al., 2021;Hashimoto et al., 2003).Moreover, we recently reported greater coefficient of variation in the levels of vesicular glutamate transporter 1 (VGlut1) and postsynaptic density 95 (PSD95), markers of E→I synaptic strength, within E→I synapses across FSIs in PFC of SZ (Chung et al., 2022).Finally, prior simulations in a QIF model network showed that these SZ-associated alterations (i.e., lower E→I strength, lower I→E strength and greater variability in E→I strength) interact to synergistically (greater than the sum of individual impact) reduce gamma power (Chung et al., 2022).In light of these findings, we explored if the synergistic interactions of these synaptic alterations can be captured by our mean-field model.First, we sought to identify a set of values to represent the parameter state in healthy individuals.To establish this parameter set, we simulated the effect of interaction between ḡE→I , ḡN→I and ḡI→E on gamma power and searched for a set of parameter values that produced maximal gamma power (Figure .4A).In this simulation, we modeled changes in the total E→I synaptic strength by adjusting ḡE→I and ḡN→I simultaneously while maintaining the fixed ratio between ḡN→I and ḡE→I at 1:10 to reflect the NMDA-to-AMPA ratio empirically observed in FSIs (Rotaru et al., 2011).Also, we set γ E→I and γ N →I at 0.01, corresponding to the lower end of the range for these parameters that results in a linear decrease in gamma power (Figure .3).The parameter values that produced maximum gamma power occurred at: where 'h' denotes the parameter value that reflects the state in healthy individuals (Table .2).
We then defined a parameter set that represents the state in SZ by referencing empirical data from prior postmortem studies, which showed percent differences in these parameters ranging from 10% to 30% relative to unaffected comparison subjects (Table .1).To prevent any single alteration from dominating the interaction effect, we uniformly applied a 20% difference to the parameter values for healthy state in order to represent a generalized model for the typical pathological state in SZ.The adjusted parameters for SZ state at 20% difference, as denoted by 'sz', were (Table .2):ḡE→I,sz = 1.04, ḡN→I,sz = 0.104, ḡI→E,sz = 0.48, γ E→I,sz = 0.012, γ N →E,sz = 0.012 (13) Next, we defined a scaling factor λ to systematically adjust the degree of alteration from the healthy to SZ states for each parameter.λ is defined as: Here, x h and x sz denote parameter values in healthy and SZ states, respectively.By varying λ from 0 to 1, the parameter x can be adjusted from its original value in the healthy state (0% difference where λ=0) to its altered value in the SZ state (20% difference where λ=1).This approach allowed us to methodologically study the effect of synaptic alterations on gamma power across a wide range of parameter values between healthy and SZ states.
Using this approach, we first investigated how varying individual synaptic parameters from healthy to SZ states affects gamma power (Figure .4B).Varying λ from 0 to 1 for ḡE→I and ḡN→I simultaneously to model lower E→I synaptic strength resulted in only a 5% decrease in gamma power.Similarly, varying λ for γ E→I and γ N →I simultaneously to model greater variability in E→I synaptic strength resulted in a 3% decrease in gamma power.Finally, varying λ for ḡI→E to model lower I→E synaptic strength resulted in an 18% deficit in gamma power.Consequently, the additive effect of all these parameter changes was expected to be a 26% reduction in gamma power.However, lower E→I synaptic strength, lower I→E synaptic strength and greater variability in E→I synaptic strength in concert reduced gamma power by 45% (Figure .4B),demonstrating a synergistic interaction among these synaptic parameters in the mean-field model.
To evaluate the robustness of this synergistic effect, we compared gamma power deficits predicted by additive changes with those from the simulated effect across a wide range of percent differences in the synaptic parameters.Our Numerical simulations of the interactions among the synaptic alterations found in PFC of SZ on gamma power.(A) Simulating the effect of interaction between E→I and I→E synaptic strengths reveals optimal gamma power (red square) at ḡE→I =1.3, ḡN→I =0.13 and ḡI→E =0.6 when γ E→I , γ N →I and γ I→E were set at 0.01; these parameter values represent the healthy state.The xaxis represents the parameter range for ḡE→I , while ḡN→I was set at 0.1 of ḡE→I across each trial to maintain the fixed ratio of these two parameters.(B) Scaling λ from the healthy state (0% difference where λ=0) to SZ state (20% difference where λ=1; see Table 2 for the parameter values for these two states) for ḡE→I and ḡN→I (open square) or γ E→I and γ N →I (open triangle) resulted in 5% or 3% deficit in gamma power, respectively, while scaling λ for ḡI→E (open circle) led to an 18% deficits in gamma power.Simulating these changes in concert (filled circle) reduced gamma power by 45%, significantly greater than the expected 26% deficits from the additive effect (filled square), demonstrating a non-linear interaction of these synaptic parameters.(C) Comparing gamma power deficits between the expected additive effect (filled square) and the simulated effect (filled circle) of the parameter interactions across a wide range of magnitude of difference shows that the non-linear interaction occurs when the parameter difference ranges from 7.5% to 32.5%.(D) The fold ratio of gamma power deficits between the simulated effect and the additive effect shows that the non-linear interaction (marked by fold ratio > 1) is greatest when the parameter difference ranges from 12.5% to 22.5% (represented by the peak of the graph).
analysis showed that simulating the combined alterations consistently produced greater deficits in gamma power than the additive predictions when the percent difference uniformly applied to all synaptic parameters was in the range of 7.5% to 32.5% (Figure .4C).Further analysis measuring the fold ratio of gamma power deficits between the simulated effect and the additive prediction Fig. 5 2D bifurcation analyses of the interactions among the synaptic alterations found in PFC of SZ on gamma power.
Each bifurcation diagram shows an interaction of two synaptic parameters on the x and y axes while varying the third parameter from its healthy to SZ states as shown on the top of each graph.The colored lines and areas represent the Hopf bifurcation points and the parameter space that support gamma oscillations, respectively, as the third parameter is varied from the healthy (blue) to SZ (red) states.The circles and the crosses indicate the values for the healthy and SZ states, respectively, for the parameters on the x and y axes.The reduction in parameter space that supports gamma oscillations is greatest when the varied parameter is ḡE→I and ḡN→I (A), intermediate with ḡI→E (B), and least with γ E→I and γ N →I (C).Thus, these findings suggest that lower E→I strength is the main synaptic pathology that drives the predicted synergistic effect on PFC gamma power in SZ.
showed that the synergy (marked by the fold ratio>1) is the greatest when the parameters differed by 12.5% to 22.5% from their healthy states (Figure .4D).Together, these findings demonstrate that the mean-field model robustly captures the synergistic property of interactions among multiple synaptic alterations found in PFC of SZ.

Two-dimensional bifurcation analyses of synergistic interactions in mean-field model
Finally, we characterized dynamical mechanisms underlying the synergistic interactions of the three SZ-associated alterations via a series of 2D bifurcation analyses of our mean-field model.Specifically, we assessed Hopf bifurcation curve of the two parameters while varying the third parameter from healthy to SZ states as shown in Table 2. 2D bifurcation diagrams of the two synaptic parameters all showed a reduction in the parameter space that supports gamma oscillations when the third parameter was varied from healthy to SZ states (Figure .5).Also, the points that marked the SZ states of the parameters on the x and y axes were placed significantly closer to the bifurcation curve than those that marked the healthy states, further highlighting the reduction in the parameter space that allows the generation of gamma oscillation.Finally, the reduction of this parameter space was the greatest when the varied parameter was E→I strength (Figure .5A),intermediate with I→E strength (Figure .5B),and least with variability in E→I strength (Figure .5C).Together, these observations suggest that although each synaptic alteration can influence the interaction between the other two, lower E→I strength is the main synaptic pathology driving the synergistic interaction and its non-linear effect on gamma power.

Discussion
In this study, we investigated the macroscopic dynamics of gamma oscillations and their alterations in SZ by deriving a mean-field description of the QIF model network that consisted of all-toall connected RSEs and FSIs coupled by AMPA, NMDA and GABA synapses.We scaled this QIF model network to the macroscopic level by transforming the membrane potential into a nondimensional variable and employing Lorentzian distributions for the synaptic strengths and external currents.This approach allowed us to derive a set of firing-rate equations for a network of heterogeneous RSE and FSI populations that are exact in the thermodynamic limit.The firing rate of the mean-field model matched the spiking activity of the finite-size QIF model network for both RSEs and FSIs over a wide range of parameter choices even though the reduction was not exact for the NMDA synapses (see below for a discussion).This close match allowed us to use the mean-field model to reliably explore changes in gamma oscillation dynamics at the macroscopic level.
Prior studies have shown that E→I and I→E synaptic strengths are the critical determinants of cortical gamma power (Sohal, 2022).To investigate how macroscopic dynamics of gamma oscillations are affected by E→I and I→E synaptic strengths in our model, we assessed the effect of ḡE→I , ḡN→I and ḡI→E on gamma power.Both the numerical simulations and the bifurcation analyses showed that ḡE→I and ḡI→E form an inverted-U relationship with gamma power, whereas increasing ḡN→I reduces gamma power.These findings replicate the effect of each synaptic parameter previously observed in the QIF model network (Rotaru et al., 2011;Chung et al., 2022), validating the utility of our mean-field model in studying the impact of synaptic parameters on the macroscopic dynamics of gamma oscillations.
Mean-field models constructed with Lorentzian distribution assume quenched variability due to the heavy-tailed nature of this distribution.In line with this idea, our mean-field model showed that increasing γ E→I and γ N →I each reduces gamma power monotonically.These findings are also consistent with a prior study that captured the effect of variability in E→I synaptic strength modeled by a Gaussian distribution in a QIF model network (Chung et al., 2022).Thus, the sensitivity of mean-field model to variability provides a computational tool for exploring how increases in synaptic variability might influence gamma oscillation dynamics at the macroscopic level.
A prior study showed that multiple synaptic alterations found in PFC of SZ (i.e., lower E→I strength, lower I→E strength and greater variability in E→I strength) can synergistically reduce gamma power due to the nature of non-linear dynamics of QIF neurons (Chung et al., 2022).This synergistic interaction was also observed in our mean-field model, demonstrating that the non-linearity of synaptic interactions is robustly captured by our model.In addition, the strength of this synergistic interaction was greatest when the magnitude of difference applied to each parameter was comparable to those reported in the postmortem human brain studies of SZ (Table.1).These findings highlight the potential utility of our mean-field model to investigate the dynamic interactions of synaptic alterations modeled by empirically findings in SZ.
The dynamic mechanism by which lower E→I strength, lower I→E strength and greater variability in E→I strength interact to synergistically reduce gamma power in SZ is not well understood.The 2D bifurcation analyses in this study revealed a reduction in the Hopf bifurcation-defined parameter space between two synaptic parameters as the third parameter shifted from healthy to SZ state.Furthermore, the reduction in the parameter space that supports gamma oscillations was most pronounced when the varied parameter was E→I strength.Thus, the shift in Hopf bifurcation due to lower E→I strength could be the primary driver of the synergistic interactions that allow small changes in multiple synaptic parameters to robustly reduce gamma power in PFC of SZ.
Our mean-field model included the dynamics of NMDA synapses, which we modeled as a slow excitatory synaptic current gated by the activity of the AMPA synapses.The dynamics of NMDA synapses typically depend on the preand post-synaptic activity, so that an exact meanfield is not technically possible.In low Mg ++ , the dependence on the postsynaptic potential is weak, so that the mean-field approximation is reasonable.Alternatively, recent work (Gast et al., 2021;Chen and Campbell, 2022;Gast et al., 2023) on short-term synaptic plasticity and spike frequency adaptation, all of which are neuron specific, has shown that a mean-field approximation can be justified if there is a dramatic separation of time scales.Since the dynamics of NMDA synapses are markedly slower than the other synaptic processes, the use of the mean-field approximation is reasonable for NMDA synapses.Consistent with this idea, our model showed a close alignment in gamma oscillation dynamics with the QIF model network that included NMDA synapses, validating our approach in modeling the dynamics of these synapses.
In conclusion, we derived a mean-field description of QIF model network to systematically explore the effect of synaptic alterations between RSEs and FSIs and their collective impact on gamma oscillation dynamics in SZ.Our mean-field model robustly captures the macroscopic dynamics of gamma oscillations and their alterations across a wide range of synaptic parameters that are thought to be the key determinants of cortical gamma power.Furthermore, our analyses reveal that the shift in Hopf bifurcation by lower E→I strength could be a critical mechanism by which multiple synaptic alterations reported in PFC of SZ non-linearly interact to synergistically reduce gamma power.As alterations in gamma oscillation dynamics are observed not only in SZ but across a wide spectrum of illnesses including bipolar disorder (Lu et al., 2022), autism (Simon and Wallace, 2016) and Alzheimer's dementia (Mably and Colgin, 2018), our mean-field model offers a powerful tool to systematically investigate the disease processes of altered neural dynamics in these diverse neuropsychiatric disorders at the macroscopic level.

Fig. 1
Fig. 1 Comparing gamma oscillation dynamics between QIF model network and mean-field model.(A) Representative average spiking activity (red) and raster plots (gray ) of QIF model network and the mean firing rate (blue) of mean-field model for RSEs (up) and FSIs (bottom) at I appl (QIF) or Ī (mean-field) = 2 .QIF network consisted of 800 RSEs and 200 FSIs.Superimposing the mean firing rate of mean-field model to the average spiking activity of QIF network shows an overlapping gamma oscillation dynamics for both RSEs and FSIs.(B) Increasing I appl and Ī in QIF network and mean-field model, respectively, results in congruent increases in gamma power and peak gamma frequency, highlighting a close alignment between the two models.Results for the QIF model network are the average of 100 trials.

Fig. 2
Fig. 2 Effect of E→I and I→E synaptic strengths on gamma power in mean-field model.(A-C) Numerical simulations showing the effect of ḡE→I (A), ḡN→I (B) and ḡI→E (C) on gamma power.ḡE→I and ḡI→E form an inverted-U relationship with gamma power, whereas increasing ḡN→I progressively reduces gamma power.(D-F) Bifurcation diagrams reveal the Hopf bifurcation points (blue dots) at ḡE→I =0.26 (D), ḡN→I =0.15 (E) and ḡI→E =1.52 (F).

Fig. 3
Fig. 3 Effect of E→I synaptic variability on gamma power in mean-field model.(A-B) Numerical simulations show a progressive decline in gamma power with greater γ E→I (A) and γ N →I (B).(C-D) Bifurcation diagrams reveal the Hopf bifurcation points (blue dots) at γ E→I =1.55 (C) and γ N →I =0.053 (D).
Fig.4Numerical simulations of the interactions among the synaptic alterations found in PFC of SZ on gamma power.(A) Simulating the effect of interaction between E→I and I→E synaptic strengths reveals optimal gamma power (red square) at ḡE→I =1.3, ḡN→I =0.13 and ḡI→E =0.6 when γ E→I , γ N →I and γ I→E were set at 0.01; these parameter values represent the healthy state.The xaxis represents the parameter range for ḡE→I , while ḡN→I was set at 0.1 of ḡE→I across each trial to maintain the fixed ratio of these two parameters.(B) Scaling λ from the healthy state (0% difference where λ=0) to SZ state (20% difference where λ=1; see Table2for the parameter values for these two states) for ḡE→I and ḡN→I (open square) or γ E→I and γ N →I (open triangle) resulted in 5% or 3% deficit in gamma power, respectively, while scaling λ for ḡI→E (open circle) led to an 18% deficits in gamma power.Simulating these changes in concert (filled circle) reduced gamma power by 45%, significantly greater than the expected 26% deficits from the additive effect (filled square), demonstrating a non-linear interaction of these synaptic parameters.(C) Comparing gamma power deficits between the expected additive effect (filled square) and the simulated effect (filled circle) of the parameter interactions across a wide range of magnitude of difference shows that the non-linear interaction occurs when the parameter difference ranges from 7.5% to 32.5%.(D) The fold ratio of gamma power deficits between the simulated effect and the additive effect shows that the non-linear interaction (marked by fold ratio > 1) is greatest when the parameter difference ranges from 12.5% to 22.5% (represented by the peak of the graph).

Table 1
Prior postmortem findings in PFC of SZ, their expected effects, and the corresponding parameters in the mean-field model.

Table 2
Comparison between synaptic parameters in mean-field model for healthy and SZ states