Mechanisms for dysregulation of excitatory-inhibitory balance underlying allodynia in dorsal horn neural subcircuits

Chronic pain is a wide-spread condition that is debilitating and expensive to manage, costing the United States alone around $600 billion in 2010. In a common type of chronic pain called allodynia, non-painful stimuli produce painful responses with highly variable presentations across individuals. While the specific mechanisms remain unclear, allodynia is hypothesized to be caused by the dysregulation of excitatory-inhibitory (E-I) balance in pain-processing neural circuitry in the dorsal horn of the spinal cord. In this work, we analyze biophysically-motivated subcircuit structures that represent common motifs in neural circuits in layers I-II of the dorsal horn. These circuits are hypothesized to be part of the neural pathways that mediate two different types of allodynia: static and dynamic. We use neural firing rate models to describe the activity of populations of excitatory and inhibitory interneurons within each subcircuit. By accounting for experimentally-observed responses under healthy conditions, we specify model parameters defining populations of subcircuits that yield typical behavior under normal conditions. Then, we implement a sensitivity analysis approach to identify the mechanisms most likely to cause allodynia-producing dysregulation of the subcircuit’s E-I signaling. We find that disruption of E-I balance generally occurs either due to downregulation of inhibitory signaling so that excitatory neurons are “released” from inhibitory control, or due to upregulation of excitatory neuron responses so that excitatory neurons “escape” their inhibitory control. Which of these mechanisms is most likely to occur, the subcircuit components involved in the mechanism, and the proportion of subcircuits exhibiting the mechanism can vary depending on the subcircuit structure. These results suggest specific hypotheses about diverse mechanisms that may be most likely responsible for allodynia, thus offering predictions for the high interindividual variability observed in allodynia and identifying targets for further experimental studies on the underlying mechanisms of this chronic pain condition.


Introduction
Understanding the neural circuitry in the spinal cord that processes pain signals is vital for understanding the mechanisms responsible for chronic pain [2,3], a wide-spread condition affecting ∼ 20% of adults in the US [1] and costing the US alone around $600 billion in 2010 [4].Indeed, the spinal cord is responsible for the initial processing of both tactile and pain-inducing stimuli at the periphery, and for relaying them to the brain [5][6][7].In particular, tactile and pain-inducing signals travel from the periphery to the spinal cord along different classes of afferent nerve fibers: Aβ-fibers which respond to innocuous stimuli such as gentle pressure or the brush of clothing on skin and C-fibers, and to some extent Aδ fibers, which respond to heat, noxious chemicals, or intense mechanical stimuli [5][6][7].Signals associated with painful stimuli, upon arriving at the spinal cord, are filtered by intermediate neural circuitry in the superficial laminae (primarily laminae I and II) of the spinal cord's dorsal horn [8].From there, the intensity of the painful signal is relayed to the brain through the firing activity of excitatory projection neurons in lamina I [8].
It is widely accepted that processing of afferent signals in the healthy pain-processing circuit of the dorsal horn relies on a balance between excitation and inhibition [2,8,9].
Melzack and Wall originally proposed this hypothesis in 1965 as the conceptual "gate control" model [10], and aspects of the theory remain influential and relevant [8,9,11].
Namely, the idea that inhibitory neurons "gate" the activity of excitatory neurons which relay pain signals towards the brain [2] remains a key hypothesis.In particular, while both excitatory and inhibitory neural populations receive input from Aβ fibers, inhibitory neurons suppress the firing of the excitatory neurons, thus blocking the response of projection neurons to pain-inducing signals from C-fibers.
Pathological changes to dorsal horn neural circuitry is frequently proposed as the culprit behind chronic pain [12,13], including allodynia, a type of chronic pain in which individuals feel pain in response to normally innocuous stimuli [14].There are several types of allodynia classified according to the type of innocuous stimuli that is painful [14].In this work, we focus on two well-studied types of allodynia: static and dynamic [15,16].In static allodynia, individuals feel pain in response to gentle pressure June 7, 2024 2/39 that normally would not be painful [16], while in dynamic allodynia individuals feel pain in response to brushing-type stimuli that normally would also not be painful [16].
Recent experimental results indicate that different types of excitatory and inhibitory dorsal horn neurons mediate static and dynamic allodynia, as discussed in Section 1.1.
Specifically, rodent experiments have shown that either static or dynamic allodynia can be induced by activating or inactivating (e.g.ablating) specific populations of excitatory or inhibitory interneurons in the dorsal horn [16][17][18][19], creating a disruption of excitatory-inhibitory (E-I) balance.In clinical conditions, however, allodynia likely occurs through more subtle circuit disruptions.Since E-I balance in a circuit can be achieved by diverse contributions of different excitatory and inhibitory neural populations, its disruption leading to allodynia can potentially occur through multiple pathways.
In this work, we use biophysically-based mathematical modeling to identify likely mechanisms by which dorsal horn neural subcircuits may be dysregulated to produce allodynia.While both Aβ fiber and C fiber activity play a role in gate control [20] and thus in allodynia, we mainly focus on the role of Aβ fiber activity in allodynia.In particular, we construct models of neural subcircuits implicated in mediating static and dynamic allodynia whose parameters are constrained to reproduce experimentally-observed behaviors.These constraints result in distributions of parameter sets representing populations of subcircuits that achieve E-I balance in different ways.We then identify the most sensitive mechanisms that disrupt E-I balance to result in allodynia in the subcircuit population.We find that the particular means of disruption varies across the subcircuit population and that the most sensitive mechanisms depend on subcircuit structure, thus predicting diverse, multiple mechanisms that may be most likely responsible for allodynia.

Proposed dorsal horn subcircuits mediating allodynia
Here, we describe subcircuit motifs that reflect recent experimental evidence for the structure of dorsal horn layer I-II networks mediating static allodynia (as evoked by a von Frey device) and dynamic allodynia (Fig 1A).Experimental studies in rodents have shown that static allodynia is reliant on activity in three putatively different types of excitatory interneurons: somatostatin-positive (SOM+) [16], calretinin-positive (CR+) [19], and protein kinase C γ-positive (PKCγ+) [17] cells.Additionally, inactivation of either dynorphin-positive (DYN+) [16] or parvalbumin-positive (PV+) [17] inhibitory interneurons is sufficient to produce static allodynia, suggesting that these inhibitory cells usually gate activation of SOM+, CR+ and PKCγ+ excitatory cells.Experiments have also identified direct synaptic connections from PV+ to PKCγ+ neurons [17].
Based on these experimental results, we consider a simplified subcircuit motif representing part of the pathway mediating static allodynia that consists of three neural populations (Fig 1B).One population, E, represents the collective activity of layer I-II SOM+, CR+ and PKCγ+ excitatory interneurons.The E population is inhibited by two distinct inhibitory interneuron populations, I1 representing DYN+ cells and I2 representing PV+ cells.All three populations receive input from Aβ fibers, (as suggested in [16,19,21]), while the E population is presumed to be additionally targeted by C fiber input (see e.g.[16,19]).
Rodent experimental studies probing the mechanisms for dynamic allodynia have identified that it relies on the activity of SOM+ excitatory interneurons that do not express the calbindin 2/calretinin gene (CR-).Specifically, ablating SOM+ neurons abolishes or greatly reduces dynamic allodynia in mice [18].On the other hand, ablating neurons expressing the calbindin 2/calretinin gene (CR+) had no effect on dynamic allodynia [18].Thus, the SOM+ neurons that are necessary for dynamic June 7, 2024 3/39 allodynia must also be negative for calretinin (SOM+/CR-).Excitatory interneurons expressing vesicular glutamate transporter 3 (VT3+) also contribute to dynamic allodynia since their ablation or silencing eliminates or attenuates dynamic allodynia induced by nerve injury or ablation in mice [18].Moreover, it has been proposed that the VT3+ neurons synapse onto the excitatory SOM+/CR-neurons [18].Further experimental results suggest that activity of the VT3+ excitatory interneurons are gated by DYN+ inhibitory neurons since VT3+ neurons fire action potentials in response to Aβ input when inhibitory signaling is blocked [18].Ablation of DYN+ cells is sufficient to induce dynamic allodynia [16] suggesting that DYN+ inhibitory cells may provide local inhibitory control of these two excitatory interneuron populations involved in dynamic allodynia.To account for these experimental results, we consider a simplified subcircuit representing part of the pathway mediating dynamic allodynia that consists of four neural populations (Fig 1C).An excitatory population E2 representing the collective activity of SOM+/CR-excitatory neurons is excited by a population E1 representing VT3+ neurons.E1 and E2 receive inhibitory input from populations I1 and I2, respectively, representing DYN+ cells.All populations receive Aβ input, (as suggested in [16,18]), and E2 is presumed to be additionally targeted by C fiber input (see [18]).The E population in the static allodynia subcircuit and the E2 population in the dynamic allodynia subcircuit are assumed to provide direct synaptic input to layer I projection cells that transmit painful signals to the brain [22].Hence, painful stimuli generating activity on C fibers would activate these excitatory cells and the downstream projection neurons.Under normal healthy conditions, responses of these excitatory populations to non-painful stimuli signaled by Aβ input is assumed to be gated by the inhibitory populations.Allodynia occurs when the E-I balance is disrupted such that non-painful Aβ input causes firing in these excitatory populations, leading to activation of the downstream projection neurons.While allodynia can also be induced by nociceptor-mediated sensitization of projection neurons [20], here we focus on allodynia caused by disruption of inhibitory gating of Aβ input.

Overview
The goal of this study is to identify and characterize "most likely" potential mechanisms that disrupt the processing of Aβ signals and E-I balance in these subcircuits leading to allodynia.We model the activity and interactions of the neural populations in these subcircuits using a population firing rate model formalism (see e.g.[23][24][25]) constrained by experimental data.We begin by identifying the full space of model parameters that generate responses to Aβ signals that is correlated to healthy conditions in each subcircuit, i.e. low firing of the E or E2 populations.Due to the relative mathematical simplicity of our firing rate model formalism, we are able to analytically define an "allowable parameter space" (APS) as a solution to a system of inequalities on model parameters.The APS for each subcircuit indicates the considerable range of parameter combinations that can reproduce healthy responses and demonstrates the myriad ways that E-I balance might be obtained.The population of subcircuits corresponding to the points in the APS represents the high variability in subcircuit structure that may occur physiologically.
Next, for each subcircuit, we define the analytic condition on model parameters that generates a response to Aβ input correlated with allodynia, i.e. high firing of the E or E2 populations.Solving this condition defines a hypersurface in the model parameter space, which we call the "allodynia surface", that separates parameter sets for which the subcircuit generates healthy, allodynia-free responses for all typical Aβ inputs, vs responses where allodynia may be present, (for at least some typical level of Aβ input).For model parameter sets (points) in the APS, we can then identify the minimal change in parameters that leads to allodynia by computing the shortest vector in parameter space from the point in the APS to the allodynia surface.The direction of these shortest vectors, corresponding to changes in specific model parameters, indicates the least variation of the subcircuit structure that results in the allodynia response.Since larger variations in parameter values can also generate an allodynia response, these shortest vectors represent the "most likely" potential mechanisms for allodynia to occur.
We then analyze the shortest vectors from the APS to the allodynia surface to identify distinct variations of parameters that dysregulate E-I balance leading to allodynia.Specifically, we find that the APS points form clusters based on parameters requiring the relatively smallest changes to reach the allodynia surface.These different clusters can be interpreted as representing different underlying mechanisms for generating allodynia.
We find that these proposed allodynia mechanisms generally involve the dysregulation of E-I balance occurring due to the release of excitatory cells from inhibitory gating or to the excitatory cells escaping the inhibitory gating.Significantly, the specific subcircuit components that are associated with the "release" or "escape" mechanisms differ in the subcircuits due to their differing network structures.As such, our results identify the diverse ways by which excitatory and inhibitory components within the subcircuits combine and interact to maintain E-I balance, and characterize the sensitivity of these subcircuits to disruptions that can lead to pathological responses.The paper is organized as follows: in Section 2.1, in order to illustrate our analysis method, we first analyze a simple subcircuit representing the canonical gate control model.In Sections 2.2 and 2.3, we then apply the analysis to the proposed subcircuits mediating static and dynamic allodynia, respectively.In Section 3, we summarize the results and put them into a greater physiological context.Details of the models and analysis methods are contained in Section 4. June 7, 2024 5/39 For our models of layer I-II dorsal horn neuronal subcircuits, we use a firing-rate model formalism [23][24][25] that describes the average membrane voltage V x (in mV) and average firing rate f x (x = E or I or x = Ej, Ij for j = 1, 2; in Hz) of the populations of excitatory and inhibitory interneurons (see Section 4).In this formalism, average voltages are governed by equations of the form: where V x,rest is the average resting voltage (in mV) and τ x is the time constant (in s) for the average voltage response of neurons in population x.Average population firing rates are computed by where the steady state firing rate activation function F ∞ x has a sigmoidal shape and its parameters are fit to experimental measurements of frequency-voltage relationships in dorsal horn neurons [26].Input to the subcircuits from Aβ fibers is modeled by a compound Poisson process with piecewise-constant rate, (see Section 4.3 for details).To fit subcircuit responses to normal healthy conditions and to identify mechanisms for allodynia, we vary the synaptic coupling strengths between populations and from the Aβ fibers, g yx (y = E, I, Ej, Ij or Aβ; in V-s).In particular, the APS for each subcircuit is a subset of the space of these parameters and the allodynia surface is a hypersurface defined within this space.

Analysis of a simple "gate control" subcircuit
We first illustrate our analysis methodology by applying it to a simple subcircuit representing the canonical "gate control" model (Fig 2A).This simple subcircuit consists of an inhibitory interneuron population (I ) that inhibits an excitatory interneuron population (E ).Both populations receive input from Aβ fibers.We assume that the E population directly synapses onto projection neurons, thus any sustained activity of the E population is a proxy for a painful response.Under normal healthy conditions, Aβ activity should not cause firing of the E population, i.e. a painful response should not occur.This is achieved through appropriate balancing of the inhibitory input from the I population and the excitatory Aβ input onto the E population.
For our analysis of this simple subcircuit, we first derive the conditions on the coupling strength parameters in the subcircuit, namely (g AβI , g IE , g AβE ), that generate normal healthy responses to Aβ inputs and define the APS.We then mathematically express the allodynia surface S that separates (g AβI , g IE , g AβE ) space into regions where the subcircuit displays healthy vs allodynia responses.Conveniently, we can visualize the APS and the allodynia surface S in the 3D parameter space, making this simple subcircuit especially helpful for illustrating our analysis method.Computing the shortest vectors in (g AβI , g IE , g AβE ) space from the parameter points in the APS to the allodynia surface S identifies the minimal parameter changes that can induce an allodynia response.By clustering the parameter points in the APS by their shortest path vectors, we characterize the different ways that the coupling strengths between subcircuit components maintain and disrupt E-I balance in the subcircuit, thus suggesting underlying mechanisms responsible for allodynia.

The APS for the simple subcircuit
To constrain the subcircuit coupling strength parameters, we impose conditions on the steady-state voltages of the E and I populations in response to constant Aβ input to June 7, 2024 6/39 account for experimental observations.To start, we require that steady-state voltages of the E and I populations, V ss E and V ss I , respectively, remain within biologically reasonable bounds, V min and V max , in response to a sustained, typical non-painful input on the Aβ fibers.Throughout this study, we assume that non-painful input corresponds to Aβ fiber firing rate f Aβ ∈ [10,20] Hz as observed in slowly adapting mechanoreceptors in response to ramp-and-hold mechanical stimulation [27].To incorporate the phenomenon of Aβ input gating responses to C fiber input, also known as pain inhibition, we require that the net signaling to the E population in response to typical Aβ activity is inhibitory, causing the voltage of the E population to decrease.
To enforce this gating response under normal healthy conditions, we require that the I population fires for all typical, non-painful f Aβ values, but the E population does not, thus representing a non-painful response.We also incorporate conditions accounting for experimental results showing that ablation of DYN+ inhibitory interneurons causes allodynia [16] by requiring that the E population fires in response to typical Aβ input in the absence of I population activity.
Each of these conditions leads to an inequality that steady-state voltages V ss E and V ss I must satisfy (Table 1).These inequalities are derived by requiring the steady state solution of the average population voltage equation (Eq. 1) in response to a constant f Aβ input remains below the maximum average voltage V x,max , above the minimum average voltage V x,min , above the voltage threshold for firing V x,thr , or below the resting voltage V x,rest for the pain inhibition condition (x = E, I; see Section 4.2 for values of these bounds).
These conditions can be re-written as a system of 6 inequalities, some nonlinear, on the coupling strengths (g AβI , g IE and g AβE ).The set of all coupling strength 3-tuples, (g AβI , g IE , g AβE ), that satisfy these inequalities constitutes the APS for this subcircuit.In particular, subcircuit instantiations in the APS do not produce allodynia for typical non-painful f Aβ signals (i.e. when f Aβ ∈ [10, 20] Hz).Because these inequalities must be satisfied for a range of input f Aβ signals, solving them explicitly requires finding the solution to a system of nonlinear optimization problems.As described in Section 4.4.1, this system for the simple subcircuit can be explicitly solved using Lambert functions.For the static and dynamic allodynia subcircuits, however, we do not find an explicit solution, although Lambert functions are used to greatly simplify the system of inequalities (Sections 4. To illustrate that the points within the APS correspond to subcircuit instantiations which yield the desired behaviors, we simulate instances of the simple subcircuit model for 20 randomly sampled APS points in response to a noisy Aβ signal with mean firing rate in [10,20]  From the distributions of the coupling strength values (g AβI , g IE and g AβE ) across the APS (Fig 2C ), we see that g AβI has the largest range of allowable values (from about 2.6 -7.1 mV/Hz), with g AβE having the second largest range (about 3.5 -6.9 mV/Hz), and g IE having the smallest range (about 0.9 -2.1 mV/Hz).Because g AβI June 7, 2024 7/39  exceeds the firing threshold for some typical non-painful, sustained f Aβ : We use this inequality to define the allodynia surface S in (g AβI , g IE , g AβE ) space, above which the corresponding subcircuit instantiation produces allodynia for at least one value of f Aβ ∈ [10,20] Hz.We use "above" in the sense that the g AβE component of a point in parameter space defines a height for that point.This surface, S, is then the following set of points: We compute S by solving this minimization problem in the unnormalized parameter space (see Section 4.7 for details).Plotting S in the normalized (ĝ AβI , ĝIE , ĝAβE ) space shows that it always lies above the APS (Fig 3A).
For each subcircuit instantiation associated with a point in the APS, we identify its vulnerability to allodynia by computing the shortest path from that point to the allodynia surface S using a customized global optimization scheme (see Section 4.8).
Because the shortest path indicates how to reach the allodynia surface by altering coupling strengths as little as possible, it represents the direction in parameter space in which the subcircuit instantiation is most vulnerable to allodynia.In this way, the June 7, 2024 9/39 components of the shortest path vector suggest which subcircuit components will need to change, along with the corresponding magnitudes of their relative change, in order to disrupt E -I balance and induce allodynia.To identify differences in vulnerabilities to allodynia (measured as changes in the parameters to reach the allodynia surface along the shortest path), we assign each point in the APS to a cluster based on the components of its shortest-path vector.The clustering algorithm, using density-based scanning [28], identifies two clusters in the APS for the simple subcircuit (Fig 3A).For points in Cluster 1 (red), the shortest direction to the allodynia surface primarily involves decreasing ĝAβI , while for points in Cluster 2 (cyan), the shortest path involves decreasing ĝIE and increasing ĝAβE .
Breaking down the points in each cluster, we find that points in Cluster Thus, the shortest paths to the allodynia surface from points in Cluster 2 are always in exactly the same direction, namely a ∼ 30% decrease in ĝIE and a ∼ 20% increase in ĝAβE .The most efficient way to induce allodynia for points in Cluster 2 involve a lowering of the inhibition onto the excitatory population through reduction of ĝIE coupled with an over-excitation of the excitatory population through an increase in ĝAβE .We characterize this mechanism for allodynia as the E population "escaping" inhibitory control.
The differences in these mechanisms for allodynia can be illustrated by considering how population firing-rate changes in response to Aβ input for coupling strengths sampled from the allodynia surface S. For example, Fig 3E (top panels) shows that a typical subcircuit instantiation in Cluster 1 with parameter values drawn from the APS (solid curves) exhibits I population firing rates around 50 Hz during non-painful Aβ input.However, when coupling strengths are set to the associated closest point on S, I population firing rates are much lower.In this way, the E population is "released" from inhibitory control and is able to fire.
On the other hand, for a typical subcircuit instantiation in Cluster 2 (bottom panels), the I population firing rate is saturated at its maximum value (80 Hz) during Aβ input for coupling strengths in the APS, as well as for coupling strengths set at their S values.Here, E population firing is promoted in the allodynia case because the response of the E population to the inhibitory input is altered through a decreased g IE (the inhibitory firing rate remains the same as in the non-painful case).Thus, the E population fires because it is able to "escape" inhibitory control.
For the simple subcircuit, the APS is basically evenly split into the 2 clusters with ∼ 52% of points in Cluster 1 and ∼ 48% in Cluster 2, suggesting that the release and escape mechanisms for allodynia are equally likely to occur.While these mechanisms for allodynia are intuitively clear and perhaps unsurprising, as shown below for the static and dynamic allodynia subcircuits, the most likely allodynia mechanism may be biased towards one of these mechanisms and can be generated by different subcircuit components when the subcircuit structure is more complex.

Analysis of the subcircuit mediating static allodynia
The static subcircuit consists of two inhibitory populations (I1 and I2 ), and one excitatory (E ) population with all three populations receiving Aβ input (Fig 4A).We assume that the E population directly relays signals to projection neurons and allodynia is defined as any firing of the E population in response to typical non-painful Aβ input.The APS is defined in the 5-dimensional space of coupling parameters (g AβI1 , g I1E , g AβE , g AβI2 , g I2E ) and the allodynia surface is a hypersurface in this space.June 7, 2024 11/39 While the I1 and I2 populations represent classes of inhibitory interneurons with different molecular markers (DYN+ and PV+, respectively), we assume they have similar response properties [26] and use the same model parameters for each population.In this way, the subcircuit is symmetric in the sense that there is nothing to distinguish the two inhibitory populations from one another.As shown below, this symmetry is reflected in the structure of the APS for this subcircuit, as well as in features of the predicted most likely mechanisms for allodynia.However, in contrast to the simple subcircuit, the most likely allodynia mechanisms are biased towards different modes of E population release from inhibitory control.

The APS for the static subcircuit
To define the APS, we impose conditions on the steady-state voltages of the E, I1, and I2 populations so that all instantiations of the static subcircuit display desired experimentally identified behaviors.In particular, we again require that all steady-state voltages remain within reasonable bounds.Since ablation of either inhibitory population can induce allodynia, we require both of the inhibitory populations to be active in order to maintain pain inhibition.Thus, under control conditions, typical Aβ stimuli induce both I1 and I2 firing, preventing E from firing, and reducing V E below its resting voltage.However, if the I1 population is ablated, then the excitatory signaling from Aβ input is strong enough to overcome the remaining inhibition from I2 to induce E firing.Likewise, if I2 is ablated, we expect E to fire.We do not incorporate ablation of both I1 and I2 in developing our conditions for the static subcircuit, as that would lead to extreme over-excitation of the E population.
Each of these conditions results in an inequality on steady-state population voltages, summarized in Table 2, that can be re-written as a set of inequalities for the coupling strengths that must be satisfied for all f Aβ input levels in the range [10,20] Hz (see Table 7).The APS for the static subcircuit is then defined as the sets of 5-tuples of coupling strengths (g AβI1 , g I1E , g AβE , g AβI2 , g I2E ) which satisfy the system of inequalities and optimization problems.As discussed in Section 4.4.2,we simplify the optimization problems using Lambert functions and then uniformly sample from the defined APS using our customized sampling algorithm (Section 4.6).
We proceed with the analysis just as we did for the simple subcircuit.We first illustrate the response of subcircuit instantiations in the APS to typical non-painful Aβ input by simulating instantiations of the subcircuit, defined by 20 points randomly sampled from the APS (Fig 4B).Subcircuit responses show that the firing rates of the I1 and I2 populations always rise to at least 25% of the maximum firing rate and often approach the 80 Hz maximum.The resulting inhibition from I1 and I2 prevents the E population from firing, causing the E voltage to drop, as required for pain inhibition.
Violin plots (Fig 4C ) of the distributions of coupling strengths in the APS highlight the symmetry of the static subcircuit.In particular, the distribution of g AβI1 is very similar to that of g AβI2 .Likewise, the distribution of g I1E is very similar to the distribution of g I2E .Also notable is that the range of g AβE is largest, followed by the ranges of g AβI1 and g AβI2 .The smallest ranges belong to g I1E and g I2E , suggesting that the subcircuit is least sensitive to changes in g AβI1 , g AβI2 , and g AβE , while most sensitive to changes in g I1E and g I2E .In addition, values of g AβE in the APS are far larger than in the simple subcircuit, reflecting the need for g AβE to balance out two sources of inhibition under control conditions, as well as to overcome the inhibition from one inhibitory population under ablation of either the I1 or I2 population.
be very strong, but because g I1E and g I2E are very small, the impact of inhibitory signaling on the E population is likely relatively weak.
We also see from the violin plots that values of g AβI1 , g AβI2 , g I1E , and g I2E are strongly bimodal.This bimodality is reflected in the two trends seen in the parallel plot (Fig 6D ) representing the sampled sets of coupling strengths in normalized parameter space.In one trend, ĝI1E is large (blacker lines), whereas the other trend (redder lines) involves larger ĝI2E values.From the parallel plot, we can also begin to see a number of correlations between coupling strengths.For instance, when ĝAβI1 is large, ĝI1E is small, and similarly, when ĝAβI2 is large, ĝI2E is small.Additionally, when ĝI1E is large, ĝI2E is small.These correlations are confirmed by the Pearson's correlation coefficients between each normalized coupling strength (Fig 6E).Generally these correlations reflect an E-I balance manifested as the excitatory (Aβ) inputs to the excitatory population balanced by a combination of inputs from the two inhibitory populations.Specifically, inputs from I1 and I2 compensate for each other where if the inhibition from I1 is large, then the inhibition from I2 is small.

Mechanisms of allodynia in the static subcircuit
To identify sensitivities to allodynia of the static subcircuit, we consider the allodynia surface S stat , which in normalized parameter space is the set of points in (ĝ AβI1 , ĝI1E , ĝAβE , ĝAβI2 , ĝI2E )-space above which subcircuit instantiations produce allodynia for at least some typical f Aβ (Section 4.7).We then compute the shortest vectors from each sampled point in the APS to S stat .Clustering of the APS based on the direction of the shortest vectors using density-based scanning (see Section 4.9) indicates that the APS divides into four clusters: Cluster 1 (green), Cluster 2 (blue), Cluster 3 (cyan), and Cluster 4 (gray) (Fig 5 Correlation Coefficients and ĝI2E values are larger than the APS mean.This reflects that subcircuit instantiations in Cluster 1 display high I1 population activity that weakly inhibits the E population, and a weakly active I2 population whose inhibitory effect on E is strong. For the subcircuit instantiations in Cluster 1, most efficiently reaching the allodynia surface primarily involves a decrease in ĝAβI2 values (Fig 5C).This parameter change reduces the weak I2 population activity and releases the E population from the inhibitory control provided by I2.The shortest vectors to S stat simultaneously involve small decreases in ĝI1E values and small increases in ĝAβE values.These parameter changes additionally promote E population firing through an escape mechanism by slightly reducing the inhibitory effect of I1 activity and slightly increasing E responses to excitatory Aβ input.The asymmetry in the inhibitory control of the E population under control conditions is maintained in the mechanism for allodynia in that the more weakly active I2 population becomes weaker and the more strongly active I1 population is not affected.This is evident in the numerical simulations of Cluster  Solid lines correspond to a subcircuit instantiation with coupling strengths given by the mean values for the particular cluster (except for Cluster 4, which is spatially disjoint, where we use a representative sampled set of clustering strengths), and dash-dotted lines correspond to the subcircuit instantiations with coupling strengths given by the corresponding closest point on the allodynia surface.Aβ input frequencies are chosen as the smallest value that induces allodynia for each cluster.

Analysis of the subcircuit mediating dynamic allodynia
In the dynamic subcircuit, allodynia is defined as firing of the E2 population in response to typically non-painful Aβ input.In this case, the APS is defined in the 7-dimensional space of coupling strength parameters (g AβI1 , g I1E1 , g AβE1 , g E1E2 , g AβI2 , g I2E2 , g AβE2 ) and the allodynia surface is a hypersurface in this space.The structure of the dynamic subcircuit consists of two simple subcircuits coupled together, one consisting of E1 and I1, and the other consisting of E2 and I2 (Fig 6A).However, our analysis shows that the most likely allodynia mechanisms do not exhibit the same symmetries as in the simple subcircuit; instead the escape from inhibition mechanism dominates.
The APS for the dynamic subcircuit Again, to define the APS, we impose conditions on the steady-state voltages of the E1, E2, I1, and I2 populations so the dynamic subcircuit displays experimentally-observed June 7, 2024 16/39 behaviors under normal healthy conditions (Table 3).Specifically, we again require that all steady-state voltages remain within reasonable bounds.In response to Aβ input in the typically non-painful range (i.e.f Aβ ∈ [10, 20] Hz), the I1 and I2 populations should fire and provide inhibitory control of the activity of the E1 and E2 populations, respectively.In particular, to account for the phenomenon of pain inhibition, we require that steady-state voltages of both E1 and E2 are hyperpolarized from resting voltage by this inhibitory input.However, if the E1 population is ablated, the E2 voltage remains within the reasonable bounds.We further require that if the I1 population is ablated, then the E1 population fires in response to typical Aβ stimuli, and in turn excites the E2 population to firing with both their voltages remaining within the reasonable bounds.Finally, if the I2 population is ablated, then we require that the E2 population fires in response to typical Aβ stimuli.
Each of these conditions leads to an inequality on steady state population voltages (Table 3) that is re-written as a set of inequalities for the coupling strengths that must be satisfied for all f Aβ input levels in the range [10,20] Hz (Table 9 in Section 4.4.3).
The APS for the dynamic subcircuit is then defined as the sets of 7-tuples of coupling strengths (g AβI1 , g I1E1 , g AβE1 , g E1E2 , g AβI2 , g I2E2 , g AβE2 ) which satisfy the system of inequalities and optimization problems.
Table 3. Conditions that the proposed subcircuit mediating dynamic allodynia must satisfy and the resulting inequalities on steady-state voltages in response to Aβ input.We ensure that the subcircuit exhibits these behaviors by imposing conditions (middle column) on the subcircuit.Each condition is exhibited in either control, E1-ablation, I1-ablation, or I2-ablation conditions (left-most column), and is realized as an inequality (right-most column) on the steady-state voltage of a population.

Condition Type Condition Steady state voltage inequality-for all f Aβ
Control conditions To illustrate the response of dynamic subcircuit instantiations in the APS to Aβ g AβE2

Discussion
While the basic tenets of Melzack and Wall's "gate control" theory [10] are still pertinent, updated conceptual models for spinal cord pain signaling are needed to account for recent results identifying diverse types of excitatory and inhibitory interneurons and their circuit structure in the dorsal horn.In this study, we analyzed biophysically-motivated subcircuit models that represent common motifs in dorsal horn layer I-II pain processing neural circuits to identify the diversity of mechanisms that maintain and disrupt E-I balance in subcircuit responses to Aβ inputs.In our model subcircuits in normal healthy conditions, E-I balance was tuned to suppress activity of excitatory interneurons, that are presumed to directly target layer I projection neurons, in response to non-painful Aβ inputs.Computation of the APS for each subcircuit defined all possible models that exhibit this healthy E-I balance and also replicate experimentally observed responses to neural ablation or silencing manipulations.To identify most likely mechanisms that disrupt E-I balance in individual subcircuits, we defined allodynia surfaces and computed shortest vectors to it from each point in the APS.The direction in parameter space of the shortest vector identified the minimal alterations to the subcircuit that resulted in a disruption of E-I balance leading to activation of the excitatory interneurons driving the pain projection neurons.
Our results show that in each subcircuit, E-I balance can be disrupted by the excitatory interneurons escaping their inhibitory control or by an attenuation of inhibitory signaling that releases them from inhibitory control, also referred to as disinhibition.Many reported changes in spinal circuits induced in physiological models of chronic pain and allodynia are examples of these mechanisms.In particular, in the escape mechanism, inhibitory signaling was unaffected but its effect on post-synaptic excitatory interneurons was diminished.Physiologically, this could occur through changes in intracellular chloride concentration in post-synaptic excitatory cells that attenuates GABA-receptor mediated inhibitory synaptic currents [29,30], down-regulation of GABA receptors or loss of synapses from inhibitory to excitatory populations [17].Additionally, escape from inhibition could occur due to increased excitability of excitatory cells [15]   ) for each population (columns) and for each cluster (rows).Solid lines correspond to a subcircuit instantiation with coupling strengths given by the mean values for the particular cluster, and dash-dotted lines correspond to the subcircuit instantiation with coupling strengths given by the corresponding closest point on the allodynia surface.Aβ input frequencies are chosen as the smallest value that induces allodynia for each cluster.
potentially arising from sprouting of synapses [31] or upregulation of neurotransmitter release from Aβ synapses [15].In the release from inhibition mechanism, inhibitory signaling is diminished which could occur physiologically through reduced pre-synaptic GABA or glycine levels, reduced GABA or glycine release, or lower firing rates in inhibitory interneurons [15,17,32].
In the simple subcircuit, model analysis predicts that E-I balance dysregulation is equally likely to occur through escape or release from inhibition, but slightly higher magnitude changes are necessary for the escape mechanism.While this result is intuitively clear due to the very simple subcircuit structure, it validates that our methodology does not have any inherent biases towards either mechanism.
In the static allodynia subcircuit, inhibitory gating of the excitatory population is distributed across two distinct interneuron populations.With this structure, our results predict that the disruption of E-I balance is biased towards reduced inhibitory signaling (i.e. release mechanisms) where the locus of decreased inhibitory signaling can occur in different parts of the subcircuit.The predicted most likely mechanism, represented by the largest cluster (Cluster 3), is when both inhibitory populations contribute equally to inhibitory control of the excitatory population and a decrease in firing or signaling June 7, 2024 21/39 occurs in both populations.The next most likely mechanisms, represented by Cluster 1 and Cluster 2, occur when one of the inhibitory populations is more weakly active than the other in healthy conditions, and loss of inhibitory control occurs when the activity of the weaker population decreases further.
In the dynamic allodynia subcircuit, Aβ input is amplified by two distinct excitatory interneuron populations, where the E2 population is downstream from the E1 population.With this structure, model results predict that disruption of E-I balance is equally likely to occur through escape (Cluster 2) or release (Cluster 1 and 3) mechanisms.The site of the escape mechanism is the E2 population, suggesting that increased activity in more downstream excitatory interneurons is more disruptive compared to similar changes in upstream excitatory interneurons.For the release mechanism, the site for loss of inhibitory control is not biased to either the upstream or downstream excitatory population but is equally likely to occur at either population.
While our analysis has centered around the minimal coupling strength changes that result in allodynia, it is important to remember that allodynia can be induced by other, higher magnitude disruptions in the subcircuit.For instance, in all three subcircuits, it is possible to induce an allodynia response by sufficiently increasing the impact of Aβ signaling on the most downstream excitatory population alone, so that it overcomes inhibitory control and fires in response to innocuous stimuli.Additionally, allodynia can be induced by sufficiently decreasing the impact of inhibitory signaling on the most downstream excitatory population.However, since larger magnitude changes to the subcircuit would be required in these scenarios, it may be presumed that they are less likely to occur than the minimal changes identified in our results.

Application to more complex dorsal horn circuits
These results can be applied to help understand how E-I balance is maintained and likely disrupted in more complex dorsal horn circuits [33,34].In models of more complex circuits, the space of unconstrained parameters is so high-dimensional that identifying the diversity of parameter sets that satisfy desired model behaviors is difficult, even when computational optimization algorithms are implemented.Our results can provide constraints on regions of parameter space where E-I balance in subcircuits of the network is achieved.For example, the dorsal horn layer I-III neural circuit modeled by Medlock et al. [33] consists of 5 excitatory interneuron populations and two inhibitory interneuron populations with each population consisting of a network of Hodgkin-Huxley-type model neurons.In the circuit, 3 of the excitatory populations act as upstream amplifiers of Aβ input to 2 downstream excitatory populations that make direct connections to projection neurons.The synaptic pathways of Aβ input through the upstream excitatory populations to one of the downstream excitatory populations are similar to the dynamic subcircuit modeled here.Our results predict that there should be parameter sets that limit responses to Aβ input in the downstream excitatory population through these synaptic pathways by mechanisms reflected by the clusters found for the dynamic subcircuit.Namely, there should be parameter sets in which the downstream excitatory neuron responses are primarily gated by direct inhibition from the inhibitory cells targeting them (corresponding to Clusters 1 and 2), and parameter sets in which their responses are limited by inhibitory control of the upstream excitatory populations (corresponding to Clusters 3 and 4).
Based on Medlock et al.'s [33] finding that, in order to maintain healthy E-I balance, small reductions in inhibitory control of upstream excitatory populations required larger increases in inhibitory input to the downstream excitatory population suggests that perhaps their model parameter set may be analogous to Cluster 3 or 4 parameter sets in the dynamic subcircuit.However, since the Medlock et al. [33] network contains additional components, there may be additional constraints on achieving E-I balance June 7, 2024 22/39 that are not contained in the smaller dynamic subcircuit.

Model limitations
A limitation of the firing rate model formalism we implemented, which only models average population voltages and firing-rates, is that specific excitability characteristics or response features that are evident on the single neuron and spike levels are not considered.Nevertheless, we did constrain the dependence of population average firing rates on membrane voltages by fitting the steady state firing rate activation functions to frequency-voltage relationships measured in dorsal horn neurons [26].However, these relationships do not take into account some spiking patterns observed in dorsal horn excitatory interneurons, such as delayed onset of firing or transient firing [16,20,35,36].
Recent development of next-generation firing rate models that can be directly reduced from networks of individual neuron models [37,38] provide the framework to include specific spiking patterns into a mean-field reduction, such as spike frequency adaptation [39,40].We expect that delayed firing could similarly be accounted for in a firing-rate model reduction of simplified neuron models that are fit to observed firing patterns of dorsal horn cells.Including such spiking properties in a firing rate population network would allow analysis of the interactions of these cellular firing patterns with network structure in maintaining or disrupting E-I balance in dorsal horn circuits.
Our models and analysis focused on the response to stimuli arriving on Aβ fibers only and did not include activity of layer I projection neurons.The model subcircuits could be extended to explicitly include the effects of painful signals arriving on C fibers that directly connect to projection neurons and interneuron populations.In previous work modeling spinal subcircuits using a firing rate model formalism, we included C fiber input that was mediated through NMDA receptors in post-synaptic populations [41].The NMDA-receptor mediated connection strength depended on post-synaptic average voltage to account for voltage-dependent removal of a Mg 2+ block.Such an extended model would be able to account for wind-up of projection neuron activity in response to repetitive brief C fiber input and also explicitly account for pain inhibition, namely the attenuation of C fiber response in the projection neurons in the presence of simultaneous Aβ input.We note, however, that our current results would not be qualitatively affected as Aβ and C fibers are parallel pathways and their inputs do not interact.

Application of methodology
In this study, we introduced an analysis methodology for neural population model circuits that can determine parameter values optimized to account for experimental observations and for identifying sensitivities of model behaviors to parameter variations.The methodology involves the following steps: 1. Translate normal and pathological experimental behaviors that a circuit should replicate into analytical conditions that model variables must satisfy.
2. Re-frame these conditions into systems of inequalities and optimization problems that the parameters of interest must satisfy.In our work, we focused on the parameters governing the coupling strengths between populations.These analytically determined conditions described distinct regions of parameter space in which the corresponding subcircuit instantiations displayed normal behaviors (the APS) and above which they displayed pathological behavior (above the allodynia surface).
June 7, 2024 23/39 3. Determine the most likely mechanisms that induce the pathological condition, in our case allodynia, by finding the shortest path from each parameter set in the APS to the surface at the boundary of the pathological region This methodology can be applied to study any circuit of neural populations, although the implementation of Step 2 is more tractable for feedforward circuits with a relatively small number of populations.Thus, we expect that the methodology may be useful for analyzing propagation of signaling and E-I balance in other sensory processing circuits with feedforward structure.

Conclusions
In all, successful treatment of chronic pain conditions, such as allodynia, require full understanding of the underlying physiological causes.Dysregulation of E-I balance in dorsal horn circuitry is a compelling cause supported by an array of pre-clinical studies, but our incomplete understanding of the circuitry and the building evidence for its complexity leave many questions for how the dysregulation occurs.The identification of multiple types of excitatory and inhibitory interneurons in dorsal horn circuits suggests that dysregulation may occur through multiple mechanisms that may be dependent on the nature of the injury or insult that induces the chronic pain condition [42].Our modeling results identifying diverse mechanisms underlying allodynia in dorsal horn subcircuits are a first step to systematically unravel the interactions among multiple interneuron populations for the maintenance of E-I balance in healthy conditions and its disruption in allodynia.Continued experimental work identifying dorsal horn circuit structure and the functional relationships among diverse interneuron types will provide constraints necessary to construct and analyze more complete circuit models that can participate in the development of therapies for this debilitating condition.

Materials and methods
This section is organized as follows: Sections 4.1 -4.3 contain descriptions of the subcircuit models and their parameter values; Section 4.4 defines the allowable parameter spaces for each subcircuit; Section 4.5 describes the algorithm used to normalize the APS; Section 4.6 describes the algorithm to uniformly sample the APS; mathematical descriptions of the allodynia surfaces for each subcircuit are contained in Section 4.7; and mathematical details for the computation of distance between APS points and the allodynia surface, and clustering of APS points based on that distance are in Sections 4.8 and 4.9, respectively.

Population firing rate model
For our models of layer I-II dorsal horn neuronal subcircuits, we implement a well-established firing rate model formalism that models the average membrane voltage and average firing rates of neuronal populations (see e.g.[23][24][25]).Average firing rates f x are computed from average voltages V x with a sigmoidal activation function of the form: where m x is the maximum firing rate of the population, β x is the half-activation voltage and α x governs the slope of the population's firing-rate response to voltage changes.We match these parameters to experimental measurements of frequency-voltage relationships for dorsal horn neurons as described in Section 4.2.We expect that in the June 7, 2024 24/39 absence of inputs, V x remains at a rest value V x,rest .In response to synaptic inputs, V x deviates from its rest value according to the following differential equation: where the time constant τ x describes how quickly V x changes in response to inputs.The inputs are computed as the sum of the firing rates of all pre-synaptic populations to population x, denoted here as y 1 , y 2 , . .., together with the firing rate of the Aβ fiber input, f Aβ , weighted by the corresponding coupling strengths: We thus expect that in the presence of steady inputs, V x approaches the steady-state value given by

Subcircuit Population
Governing Equation simple subcircuit

Parameters of population firing-rate models
We choose parameters for the activation functions of excitatory and inhibitory populations based on the experimental measurements of membrane properties and firing behavior in rat dorsal horn neurons reported in Ruscheweyh et al. [26].In our subcircuits, we assume all excitatory populations have the same parameters and assume the same for inhibitory populations.We assume that average resting voltages V I,rest and V E,rest are −60 mV, approximately the values reported in [26] Maximum firing rates of excitatory and inhibitory populations are set to 50 and 80 Hz respectively, based on [33].
We use the frequency-current relations and current-voltage relations reported in [26] to extract frequency-voltage relations.The firing-rate activation functions given in Eq (2) are fit to these frequency-voltage relations using the trust-region-reflective non-linear-least-squares algorithm via Matlab's "fit" function [43] to obtain the values of β x and α x (x = E, I) for excitatory and inhibitory populations.
Maximum and minimum voltages are set to 12α x mV above and below, respectively, the half-activation voltage value β x .The voltage thresholds for firing are defined as June 7, 2024 25/39 β x − α x .This ensures that both populations have similar behavior in response to proportional changes in input.Finally, we choose membrane time constants so they are roughly on the same time-scale as those of [44].Notably, since many of the results depend on steady-state voltages, small changes in the particular values of the time constants τ E and τ I have little effect on the qualitative behavior of the results.We choose τ E = 0.01 seconds and τ I = 0.02 seconds.All model parameters are listed in Table 5.

Aβ stimuli model
In simulations of subcircuit responses to time-varying input on Aβ fibers, we drive neural populations with input representing the average firing rate of a bundle of 300 Aβ fibers (which is on the order of magnitude of the number of Aβ fibers in afferent nerve fibers from rat skeletal muscle, as in e.g [45]).In particular, we describe the spiking activity on each Aβ-fiber in the bundle via a Poisson process with a variable rate depending on the presence or absence of peripheral stimuli.Namely, in the absence of stimuli, each fiber transmits action potentials at a background rate of 1 Hz, while in the presence of a stimulus the rate increases to f Aβ Hz.The input to subcircuit populations is the average spiking rate across all Aβ fibers (in Hz), namely a noisy time-varying input with average firing rate f Aβ Hz.In all subcircuit simulations, stimuli are applied from 0.2 − 0.7 seconds.
June 7, 2024 26/39 In the parameter sensitivity analysis of model subcircuits, f Aβ is taken as a constant value between [10,20] Hz.This range of Aβ fiber firing activity has been observed in slowly adapting mechanoreceptors in response to ramp-and-hold mechanical stimulation [27].

Defining the allowable parameter space (APS) for the subcircuits
Our parameter sensitivity analysis method consists of translating normal experimental behaviors that the subcircuits should replicate into analytical conditions on model variables, specifically on average voltages.These conditions are then re-written into systems of inequalities and optimization problems that coupling strength parameters must satisfy.The parameter sets that satisfy these systems constitute the allowable parameter space (APS).In this section, we derive these systems for our model subcircuits.Full details of the derivation are shown for the simple subcircuit only, as the derivation follows similarly for the static and dynamic subcircuits.

Simple subcircuit
For the simple subcircuit, it is possible to make considerable progress towards deriving an explicit definition of the allowable parameter space.As described in Table 1 in Section 2.1, the conditions on average voltages that ensure the subcircuit replicates experimentally appropriate behaviors are given as follows: The first two inequalities can be rewritten to yield bounds on g AβI : However, as these conditions must hold for all f Aβ ∈ [10, 20] Hz, we need that max which can be re-written as Likewise, the last two inequalities in Eq (4) yield analogous bounds on g AβE and the middle two inequalities yield analogous bounds also on g AβE .We summarize the resulting system of inequalities for coupling strength parameters in Table 6.
Since f I is a nonlinear (hyperbolic tangent) function of f Aβ and g AβI , the optimization problems in the last line of Table 6 are generally difficult to solve explicitly.Nevertheless, we can simplify the maximization problemmax f Aβ -and in fact explicitly solve the minimzation problem -June 7, 2024 27/39 f Aβ -making it much easier to numerically approximate the APS.To do so, note that both of those optimization problems can be rewritten as: Consequently, the solutions to the preceding optimization problems occur either at x = g AβI f Aβ,min /α I , g AβI f Aβ,max /α I , or at one of the critical points given by We show in S1.1 Appendix that when C = 0, the solutions of Eq 5 are given in terms of the 0 th , W 0 , and −1 st , W −1 , branches of the Lambert-W function: which yields In S1.1 Appendix,we address the case when C = 0 by showing how to take advantage of the structure of the optimization problem to solve it numerically.
As an alternative approach, it is possible to address the case where C = 0 by rewriting the middle two equations of Eq 4 so they bound g IE (see S1.2 Appendix).
The resulting optimization problems may then be solved explicitly in terms of the Lambert W functions.While we do not implement this alternative approach for the simple circuit, we do apply it in our analysis of the static subcircuit (see below).

Static subcircuit
For the static subcircuit, the inequalities that population voltages must satisfy to replicate experimentally-observed behaviors are listed in Table 2 in Section 2.2.
June 7, 2024 28/39 Table 7. Inequalities expressed as upper and lower bounds on coupling strengths define the allowable parameter space for the static subcircuit.These inequalities are obtained by algebraically manipulating the inequalities on the voltages of various populations from Table 2 that define the APS for the static subcircuit so that the inequalities are written explicitly in terms of coupling strengths.(UB = upper bound, LB = lower bound, abl = ablation, PI = pain inhibition).

I1
Fires + UB Following a similar derivation as for the simple subcircuit, we rewrite the conditions as the system of inequalities and optimization problems on coupling strength parameters given in Table 7.The inequalities in the last three rows of Table 7 involve maximizing or minimizing a quantity over the range of f Aβ values.To make the APS easier to compute, we find explicit solutions to the four optimization problems in the last two rows of Table 7 using the alternative approach described above for the simple subcircuit and outlined in S1.2 Appendix and S1.3 Appendix.The solutions to the four optimization problems (Table 8, 3rd column) can be written in terms of the Lambert W 0 function (4th column) with different constants A (5th column).
To more easily sample from the APS, it is helpful to compute bounds on E population coupling strength values that are independent of other coupling strength parameters.For example, to find an upper bound on g I1E , we can use the upper bound on E during I1 -ablation, (rewritten so the inequality is expressed as bounds on g AβE rather than g I1E ), along with the E lower bound given no ablations to obtain that which requires in particular that for all June 7, 2024 29/39 which in turn implies that To find a lower bound on g I1E , on the other hand, we can use the E lower bound under I1 ablation, (rewritten so the inequality is expressed as bounds on g AβE rather than g I1E ), along with the pain inhibition condition, i.e. that: which, via an argument analogous to the preceding argument used to find an upper bound on g I1E , implies that An analogous procedure produces bounds on g I2E .

Dynamic subcircuit
We rewrite the inequalities in Table 3 in Section 2.3 as the system of inequalities and optimization problems for the coupling strength parameters given in Table 9.
Table 9. Inequalities expressed as upper and lower bounds on coupling strength values define the allowable parameter space for the dynamic subcircuit.These inequalities are obtained by algebraically manipulating the inequalities on the voltages of various populations from (Table 3) that define the APS for the dynamic subcircuit so that the inequalities are written explicitly in terms of coupling strength parameters.
(UB = upper bound, LB = lower bound, abl = ablation, PI = pain inhibition) We further simplify the inequalities in the second line of Table 9 by explicitly solving the optimization problems in terms of Lambert-W functions analogously to our treatment of the optimization problems in the last row of (Table 6) for the simple subcircuit.Upper and lower bounds on g I1E1 are straightforward to find, because the I1 -E1 portion of the dynamic subcircuit has identical constraints to those of the simple subcircuit.To find upper and lower bounds on g I2E2 , we use a more computationally intensive approach.Namely, given the parameter values for the I1 -E1 portion of the subcircuit and given g E1E2 , we find the set of g I2E2 values such that the inequalities in the bottom-most three lines of Table 9 have a solution.That is, we choose g I2E2 so that the upper bounds on g AβE2 are indeed larger than the lower bounds appearing in the last three lines of Table 9.
June 7, 2024 30/39 4.5 Normalizing the allowable parameter space (APS) We normalize the APS so that changes in different coupling strength parameters can be compared.Since normalization requires upper and lower bounds on each coupling strength parameter, we construct a rectangular hypercube in parameter space that contains the APS.To do this, we leverage the hierarchical nature of the sets of inequalities on coupling strength parameters for each subcircuit in Tables 6, 7 and 9.
The hierarchy is formed by the dependencies of inequalities on the coupling strength parameters, with inequalities higher in the hierarchy depending on parameters defined by inequalities lower in the hierarchy.Table 10 lists coupling strength parameters for each subcircuit in the order of the hierarchy formed by their inequalities (from low to high).Here we describe the algorithm implemented to construct the rectangular hypercube in parameter space that contains the APS based on the hierarchy of inequalities on coupling strength parameter values.In general, consider a hierarchical set of inequalities that define bounds on elements of a parameter vector x = {x 1 , . . ., x n } in R n : for real numbers a 1 ≤ b 1 , . . .a n ≤ b n , and real functionals L 1 ≤ U 1 , . . ., L n ≤ U n .As written, the x n inequalities are the highest in the hierarchy because they depend on x 1 , . . ., x n−1 .To identify the interval of values [x i,min , x i,max ] that satisfies the inequality for each x i , we generate a large (at least 1000 elements), random (not necessarily uniform) sample of x values satisfying the inequality system as follows: • Given the value of x 1 , uniformly at random choose a value of x 2 ∈ [L 2 (x 1 ), U 2 (x 1 )] if such an interval exists.If such an interval doesn't exist, start over.
• Repeat to choose x 3 , . . ., x n sequentially.If, at any step, an interval for x i is not defined, start over with a new choice for x 1 .
June 7, 2024 31/39 We define the minimum (maximum) value of each coupling strength parameter in the APS as the minimum (maximum) value across all samples.This defines a hypercube in R n that contains the APS.To normalize, we map each element of x to [0, 1] as follows: xi = x i − x i,min x i,max − x i,min .
The normalized APS thus lies in the unit hypercube.

Sampling from the allowable parameter space (APS)
To perform analyses, we generate a uniform sampling of the APS.However, the APS is high-dimensional, can be non-convex, and is generally of an unknown shape, making it difficult to sample uniformly.In this section we briefly outline the algorithm we developed to sample the APS for each subcircuit; further details are contained in S2 Appendix.
Having found a rectangular subspace that contains the APS (Section 4.5), one approach would be to uniformly at random sample a point from the rectangular subspace, check to see if the point is in the APS, and keep it if it is.However, this naive sampling scheme is dependent on the volume of the APS in such a way that its computational complexity is exponential in n (see S2.2 Appendix).

Volume-independent sampling algorithm
Here we describe a spatially uniform sampling algorithm whose computational complexity is independent of the volume of the APS in R n .To do so, we sample from a cover of the normalized APS consisting of R n hyperrectangles that approximates the normalized APS.Our algorithm takes the following strategy: 1. Define a set of R n hyperrectangles that contains and approximates the allowable parameter space.
2. Uniformly at random select a hyperrectangle 3. Uniformly at random select a point from the hyperrectangle.
• If the point is not in the parameter space, discard it 4. Otherwise, keep the point with probability proportional to the volume of the hyperrectangle.
Step (4) ensures that the sample is indeed uniform-in-space.

Defining the allodynia surface
In this section, we derive the conditions describing the allodynia surface for each subcircuit.We present the derivation in terms of a generalized circuit and conditions for a general target state to occur.
In particular, we restrict our attention to circuits where the output signal is relayed by a single neural population, which we denote by y.The target state is represented by the voltage of this output population V y increasing above a specified threshold in June 7, 2024 32/39 response to an input signal.For our firing rate model, the steady-state average voltage V y is given by V y,steady = g in,y f in + n k=1 g x k y f x k ( g, f in ) + V y,rest .
where f in is the input signal and f x k are firing rates of the populations pre-synaptic to population y.The vector g contains all the coupling strength parameters in the circuit, g ji for the weight of the connection from pre-synaptic population j to postsynaptic population i.The circuit is in the target state when V y,thr ≤ V y,steady .
In our subcircuits, the output populations y are the excitatory interneuron populations that directly target the projection neurons that relay signals to the brain, namely the E, E, and E2 populations in the simple, static and dynamic subcircuits, respectively.The target state for allodynia is that the average voltage of these excitatory populations increases above the firing threshold in response to Aβ input.These allodynia conditions are summarized in Table 11.
Table 11.Allodynia conditions on each subcircuit.Allodynia conditions specifying precisely when the subcircuit is producing allodynia-when it is relaying pain-inducing stimuli in response to innocuous f Aβ signals.

Subcircuit
Conditions under which allodynia is produced Using the target state condition, we can identify the coupling strength values g for which the target state is attainable.To do so, we rewrite the condition defining the target state as an inequality on the coupling strength g in,y between the input signal and the population y: g in,y f in ≥ (V y,thr − V y,rest ) − n k=1 g x k y f x k ( g, f in ) ⇔ (7) g in,y ≥ (V y,thr − V y,rest ) − To illustrate this more concretely, the conditions on the coupling strengths for which allodynia is attainable are summarized in Table 12 for the simple, dynamic, and static subcircuits.
Table 12.Conditions for a subcircuit to produce allodynia.These conditions specify the sets of coupling strength values when the corresponding subcircuit will produce allodynia for some typical f Aβ input.

Subcircuit
Eq (9) thus defines a target state boundary surface S which divides the sets of coupling strength values for which the circuit is in the target state from those for which it is not.We can express this surface as S = g : g in,y = min fin∈[fin,min,fin,max] (V y,thr − V y,rest ) − n k=1 g x k y f x1 ( g, f in ) f in (10) For the simple, static, and dynamic subcircuits, this boundary (Table 13) separates regions of parameter space in which the subcircuit can produce allodynia from regions where it does not.In this section, we discuss the computational algorithm that computes the shortest path from points in the APS to the allodynia surface.Similarly as in Section 4.7, we describe the algorithm for a generalized circuit and a generalized target state boundary surface S.
June 7, 2024 34/39 The length of the shortest path between a point g in the space of coupling strength values to the target state boundary surface S indicates how easy it is to move the circuit into the target state.It also identifies which coupling strengths need to change to reach the target state.We illustrate the problem of finding the shortest path to S by considering an arbitrary set of coupling strength values-a point in the APS which we denote gp.To compute the shortest path from gp to S, we need to find the point gs nearest on S closest to gp.To do so, we need to solve the optimization problem: where || • || represents the Euclidean norm.
Solving the optimization problem in Eq (11) directly would require knowing the target state surface itself or knowing important properties of it such as its gradient.
However, the target state boundary surface in our work is defined via solving a minimization problem over the space of coupling strength values (excluding the coupling strength g in,y ) and is thus difficult to include in existing optimization algorithms.Thus, we seek to solve this problem without computing an explicit representation of the target state surface S.
where S * is the following set of coupling strength-input signal pairs defined by removing the minimization from Eq (10): S * := ( g, f in ) : g in,y = (V y,rest − V y,thr ) + n j=1 g xiy f xi ( g, f in ) f in and Eq (12) presents a constrained optimization problem that can be solved with high likelihood using global optimization algorithms based on stochastic gradient descent.As a result, Eq (12) is far more tractable than the original problem posed in Eq (11).
Moreover, if gs * nearest = gs nearest , the two minimization problems are equivalent, and by solving Eq (12), we will have solved the original minimization problem Eq (11).In S3 Appendix, we describe how we solve Eq (12) using Matlab's stochastic gradient descent-based algorithm fmincon in a multi-start global optimization scheme and show that if we solve Eq (12), we do indeed solve Eq (11).
4.9 Clustering the data based on shortest paths to the allodynia surface We use a density-based scanning clustering algorithm coupled with data visualization to identify clusters in the APS.To do so, we work with the uniformly sampled points in the APS and cluster them according to their shortest paths to the allodynia surface.
To apply density-based clustering to the data, we use Matlab's dbscan function.
Briefly, dbscan divides the data into equivalence classes, where two datapoints are equivalent if they are sufficiently close, and identifies equivalence classes as a cluster if they contain a point which exceeds a minimum number of sufficiently close neighbors.The function takes three arguments: 1.The data to be clustered: We use the set of shortest path vectors from sampled points in the APS to the allodynia surface.
June 7, 2024 35/39 2. A sufficiently close distance : we use the smallest under the euclidean metric that leads to no outliers in the data.
3. The minimum number of sufficiently close neighbors to identify an equivalence class as a cluster: we take this to be 5.

Fig 1 .
Fig 1. Proposed subcircuits in layer I-II of the dorsal horn mediating static and dynamic allodynia.(A) Schematic of neural circuitry in layer I-II of the dorsal horn consisting of populations of inhibitory (magenta shades) and excitatory (blue shades) interneurons that filter Aβ and C fiber inputs to projection neurons (green) that transmit signals to the brain.(B) A schematic of the proposed subcircuit mediating static allodynia.I1 and I2 represent the populations of dynorphin-positive (DYN+) and parvalbumin-positive (PV+) inhibitory interneurons, respectively, and E represents a collective population of somatostatin-positive (SOM+), calretinin-positive (CR+) and protein kinase C γ-positive (PKCγ+) excitatory interneurons.(C) A schematic of the proposed subcircuit mediating dynamic allodynia.I1 and I2 represent populations of DYN+ inhibitory neurons, and E1 and E2 represent the populations of VGLUT3-positive (VT3) and somatostatin-positive/calretinin-negative (SOM+/CR-) excitatory interneurons.In all panels, Aβ and C represent inputs relayed from the periphery along Aβ and C fibers, respectively.
4.2 and 4.4.3).Instead of explicitly solving this system of optimization problems, we compute a uniformly-distributed sample of the APS.Because the APS may be non-convex (see Fig 2D), we must develop a sampling algorithm (Section 4.6) to generate the uniform sample shown in Fig 2C -Fig 2D.
Hz.The subcircuit responses show that the firing rate of the I population rises to at least 80% of its maximum of 80 Hz (Fig 2B, left panel) and the inhibition from the I population causes voltages of the E population to drop (Fig 2B, middle panel), as we require for pain inhibition.Thus, as expected, E population firing rates remain low except for a brief increase at the onset of the Aβ stimulus (Fig 2B, right panel).On the other hand, if we simulate ablation of the I population (by setting g IE to zero), the voltage of the E population sufficiently increases so that the E population fires in response to non-painful Aβ stimuli, thus simulating allodynia, (see S1 Fig in the supplementary material).

Fig 2 .
Fig 2. Allowable parameter space (APS) for the simple "gate control" subcircuit.(A) A schematic of the simple "gate control" subcircuit, where E and I represent populations of excitatory and inhibitory interneurons, respectively.(B) Simulation results showing the mean (black lines) and range (shaded gray areas) for the firing-rate responses of the I population (left panel), and the average voltage (middle panel) and firing rate (right panel) of the E population, simulated for 20 sampled points in the APS each with a different random Aβ stimulus in range [10, 20] Hz (during t ∈ [0.2, 0.7] s).(C) Violin plots showing the distribution of each coupling strength parameter with mean (square marker) and range of values that lie within one standard deviation (vertical bar) indicated.(D) A scatter plot of 5000 uniformly sampled APS points in the normalized (ĝ AβI , ĝIE , ĝAβE ) space.(E) Normalized Pearson correlation coefficients between coupling strength parameters in the normalized APS sample.

Fig 3 .
Fig 3. Allodynia mechanisms represented as shortest paths from the APS to the allodynia surface in the simple "gate control" subcircuit.(A) A scatter plot of the sampled APS points (5000 points) in the normalized (ĝ AβI , ĝIE , ĝAβE ) space with the allodynia surface S overlaid.Based on the direction of their shortest paths, APS points separate into 2 clusters (Cluster 1 (red) and Cluster 2 (blue)).The nearest point on the allodynia surface from each sampled APS point is also shown (darker red and blue points on S). (B) Violin plots of the coupling strength distributions for the APS points in each cluster.Black * shows the mean APS values, colored square shows the mean values in each cluster.(C) The probability distribution of the shortest distances to the allodynia surface (overall profile).The shading represents the contributions from each cluster to the overall profile.For instance, a bar that is 70% light blue indicates that cluster 2 constitutes 70% of subcircuit instantiations with the corresponding distance to the allodynia surface.(D) Parallel plot representation of the components of the shortest path vectors from APS points to their corresponding nearest points on S, colored according to cluster membership.(E) Firing rate responses f E (left panels) and f I (right panels) to a noisy Aβ input signal with amplitude f Aβ ∈ [10, 20] Hz occurring during t ∈ [0.2, 0.7]s for each cluster, (top row: Cluster 1, bottom row: Cluster 2).Solid lines correspond to the subcircuit with coupling strength values set to the mean values for the cluster and dash-dotted lines correspond to the simple subcircuit instantiation with coupling strengths set to the corresponding closest point on the allodynia surface.Aβ input mean frequencies are chosen as the smallest value that induces allodynia for each cluster.
1 have larger values of ĝIE and ĝAβE and smaller values of ĝAβI relative to the mean of each coupling strength in the whole APS (Fig 3B).We also see that points in Cluster 1 are June 7, 2024 10/39 generally closer to the allodynia surface than Cluster 2 (Fig 3C), indicating that points in Cluster 1 are more sensitive to dysregulation than points in Cluster 2. In terms of the mechanism for allodynia, the shortest paths from Cluster 1 points to S consist of a decrease in ĝAβI and an increase in ĝAβE (Fig 3D).Thus, the most efficient means of producing allodynia for instances of the subcircuit in Cluster 1 involve disinhibition of the E population by lowering the response of the I population to Aβ input by reducing ĝAβI , and over-exciting the E population by increasing ĝAβE (Fig 3E).We characterize this mechanism of disrupting the E-I balance as the E population being "released" from inhibitory control.Points in Cluster 2, on the other hand, are characterized by larger ĝAβI values and smaller ĝIE and ĝAβE values (Fig 3B), relative to all points in the APS.Points in Cluster 2 are typically farther from the allodynia surface, indicating that these instances of the subcircuit may be more protected against E -I disruption (Fig 3C).The shortest paths from points in Cluster 2 to the allodynia surface involve decreasing ĝIE and increasing ĝAβE (Fig 3D) in a fixed ratio, in contrast to points in Cluster 1 (see S4 Fig).
Further, as seen in Fig 4C and Fig 4D, there are a few points with values of g AβI1 and g AβI2 much larger than is typical.This indicates that there are long and narrow regions of the APS with very low-volume, appearing for high values of g AβI1 and g AβI2 .In the corresponding subcircuit instantiations, firing rate responses of I1 and I2 would June 7, 2024 12/39

Fig 4 .
Fig 4. Allowable parameter space (APS) for the proposed subcircuit mediating static allodynia.(A) A schematic of the static subcircuit.I1 and I2 represent the populations of inhibitory neurons and E represents the population of excitatory neurons that synapse onto projection neurons.Aβ represents inputs to the subcircuit relayed from the periphery along Aβ fibers.(B) The mean (black lines) and range (shaded gray areas) for the firing rate (top left) and voltage (top right) of the E population, as well as for the firing rates of the I1 (bottom left) and I2 (bottom right) populations calculated from 20 sampled sets of coupling strengths each with a different random Aβ input stimulus (active during t ∈ [0.2, 0.7] s).(C) Violin plots showing the distribution of each coupling strength parameter with mean (square marker) and range of values that lie within one standard deviation (vertical bar) indicated.(D) Parallel plot representation of the sampled sets of normalized coupling strengths.A line gives the values of each coupling strength in a sampled set, colored on a gradient from light red to dark red according to g AβI1 so as to easier differentiate between individual lines.(E) Normalized Pearson correlation coefficients between sampled sets of normalized coupling strengths.
1 subcircuit instantiations shown in Fig 5D (top row).In particular, a typical Cluster 1 subcircuit instantiation exhibits a decrease in I2 activity in the allodynia condition compared to control, but no change in I1 activity.Subcircuit instantiations in Cluster 2 (blue) exhibit the opposite asymmetry in June 7, 2024 14/39 inhibitory control of the E population with higher ĝAβI2 values compared to ĝAβI1 and lower ĝI2E values compared to ĝI1E (Fig 5A, 2nd panel from the top).Thus, in these subcircuit instantiations, the I1 population is weakly active compared to I2 and the direction of the shortest vectors to S stat primarily involves decreases in ĝAβI1 .This results in lower I1 activity in the allodynia condition and release of the E population from I1 inhibitory control (Fig 5C and D).Clusters 1 and 2 combined make up about half of the APS points and they display approximately the same mean distances to the allodynia surface (Fig 5B).While subcircuit instantiations in Clusters 1 and 2 both exhibit asymmetry in the activity levels of the two inhibitory populations, in Cluster 3 (cyan) the inhibitory control of the E population is generally equally shared by I1 and I2.Indeed, the coupling strength distributions are very similar between ĝAβI1 and ĝAβI2 , as well as between ĝI1E and ĝI2E (Fig 5A, third panel).Cluster 3 is the largest cluster and its points lie closest to the allodynia surface (Fig 5B), suggesting that an equally distributed inhibitory gate on E activity is the most likely configuration under our conditions, but that it is also most easily disrupted.To induce allodynia in Cluster 3 subcircuit instantiations, it is most efficient to decrease both ĝAβI1 and ĝAβI2 by roughly equal amounts and to increase ĝAβE (Fig 5C), reflecting release from both I1 and I2 inhibitory control.As shown in Fig 5D, (third row), for typical subcircuit instantiations in Cluster 3, neither I1 or I2 firing rates are at their maximum values in control or allodynia conditions, and slight decreases in both their firing rates is sufficient to allow E firing in the allodynia condition.Clusters 1, 2, and 3 make up almost all of the APS, but ∼ 2% falls into Cluster 4 (gray).Points in Cluster 4 are also symmetric in the distribution of coupling strengths related to the inhibitory populations, and they are more bimodal than the APS as a whole.This is particularly reflected in the distributions of ĝAβI1 and ĝAβI2 (Fig 5A, bottom panel), suggesting that Cluster 4 has at least two groups of subcircuit instantiations that are far apart in parameter space, as is evident in the parallel plot displayed in S3 Fig.In fact, these two groups show similarities to parameter values in Cluster 1 or in Cluster 2, respectively.However, unlike Clusters 1 and 2, allodynia is most easily induced for all Cluster 4 subcircuit instantiations by decreasing both ĝI1E and ĝI2E and increasing ĝAβE , reflecting an escape of the E population from inhibitory control (Fig 5C).As shown in Fig 5D, (bottom panel), a typical subcircuit instantiation in Cluster 4 has inhibitorypopulations saturated near the maximum firing-rates, despite the asymmetry mentioned above.Thus, decreasing ĝI1E and ĝI2E reduces the effect of the two inhibitory populations.The fact that Cluster 4 is very small and its subcircuit instantiations are considerably further from the allodynia surface compared to the other clusters (Fig 5B)suggests that this mechanism for E-I balance and its disruption is less likely to occur in this subcircuit structure.In summary, the most likely mechanisms for allodynia in the static subcircuit involve release of the E population from inhibitory control where the locus of decreased inhibitory signaling can occur in different parts of the subcircuit.The largest cluster (Cluster 3), and thus, the predicted most likely mechanism, is when activity levels are similar in the I1 and I2 populations, and their simultaneous decrease releases the E population to activate in response to non-painful Aβ input.The other more likely mechanisms (Cluster 1 and Cluster 2) occur when one of the inhibitory populations activates more weakly than the other in control conditions, and allodynia occurs when the activity levels of the weaker population decrease further.These results suggest that when inhibitory gating of excitatory cell activity is distributed among separate neuronal populations, the disruption of E-I balance is biased towards reduced inhibitory signaling (i.e. release mechanisms) compared to enhanced excitatory cell responses (i.e.escape June 7

Fig 5 .
Fig 5. Mechanisms for static allodynia represented as shortest paths to the allodynia surface for the static allodynia subcircuit.(A) Violin plots for coupling strength distributions in the 4 clusters of shortest path lengths to the static allodynia surface: Cluster 1 (green, top panel), Cluster 2 (blue, 2nd panel), Cluster 3 (cyan, 3rd panel), and Cluster 4 (gray, bottom panel).Black * represents the mean APS coupling strength value, colored square is mean value for the cluster.(B) Probability distribution of the shortest distances to the allodynia surface (overall profile); shading represents the contributions from each cluster to the overall distribution.For instance, a bar that is 70% light blue indicates that cluster 3 constitutes 70% of subcircuit instantiations with the corresponding distance to the allodynia surface.(C) Parallel plot representation of the shortest paths to the allodynia surface from each sampled point from the APS.A line gives the components of the displacement vector corresponding to one such shortest path.Lines are colored based on clusters with a color gradient for better visualization.(D) Firing-rate responses to a noisy Aβ input signal (active during t ∈ [0.2, 0.7] s) for each population (columns) and for each cluster (rows).Solid lines correspond to a subcircuit instantiation with coupling strengths given by the mean values for the particular cluster (except for Cluster 4, which is spatially disjoint, where we use a representative sampled set of clustering strengths), and dash-dotted lines correspond to the subcircuit instantiations with coupling strengths given by the corresponding closest point on the allodynia surface.Aβ input frequencies are chosen as the smallest value that induces allodynia for each cluster.

input, we simulate 20
instantiations of the subcircuit with noisy Aβ input (Fig 6B;average f Aβ in[10,20] Hz).Firing rates of the I1 and I2 populations rise to nearly their maximum firing rates of 80 Hz preventing E1 and E2 from firing and decreasing E2 average voltage, as required for pain inhibition.June 7, 2024 17/39 Violin plots of the sampled APS points (Fig 6C) show that the coupling strengths governing the response of the populations to Aβ input, namely g AβI1 , g AβE1 , g AβI2 and g AβE2 , have the largest ranges, reflecting low sensitivities of these parameters for obtaining normal healthy subcircuit responses.On the other hand, the coupling strengths from the inhibitory populations to the excitatory populations, namely g I1E1 and g I2E2 , have the smallest ranges, suggesting that excitatory population responses to inhibitory signaling are more constrained to maintain healthy subcircuit responses.Correlations between coupling strength values indicate how E-I balance is maintained in the subcircuit.A parallel plot of normalized parameter sets (Fig 6D) indicates positive correlations by flat lines between two coupling strength values such as between ĝI1E1 and ĝAβE1 , and between ĝI2E2 and ĝAβE2 .These strong positive correlations are likewise apparent in computed correlation coefficients (Fig 6E) and reflect the balance of responses to inhibitory and excitatory inputs by the E1 and E2 populations, similar to what we found for the simple subcircuit.Also similar as in the simple subcircuit, negative correlations between g AβI1 and g I1E1 , as well as between g AβI2 and g I2E2 , indicate a preservation of inhibitory signaling to each excitatory population such that weak responses of the inhibitory populations are compensated by higher sensitivity of the excitatory populations to their activity.Further, coupling strengths pertaining specifically to the I2 -E2 component are only weakly correlated with coupling strengths pertaining to the I1 -E1 component, indicating that excitation and inhibition are being balanced separately within each subcircuit component.

Fig 6 .
Fig 6.Allowable parameter space (APS) for the proposed subcircuit mediating dynamic allodynia.(A) A schematic of the dynamic subcircuit.I1 and I2 represent the populations of inhibitory neurons, and E1 and E2 represent the populations of excitatory neurons.We assume the E2 population relays signals to projection neurons.Aβ represents inputs to the subcircuit, relayed from the periphery along Aβ fibers.(B) The mean (black lines) and range (shaded green areas) for the firing rate (top left) and voltage (top right) of the E2 population, as well as for the firing rates of the I2 (middle left), E1 (middle right), and I1 (bottom left) populations calculated across 20 sampled sets of coupling strengths each with a different random input Aβ stimulus (active during t ∈ [0.2, 0.7] s).(C) Violin plots of distributions of coupling strength values in the APS.The vertical bar represents the values within one standard deviation of the mean (square marker) for the corresponding coupling strength.(D) Parallel plot representation of the sampled sets of normalized coupling strengths.A line gives the values of each coupling strength in a sampled set, colored on a gradient from light red to dark red according to g AβI1 .(E) Normalized Pearson's correlation coefficients between coupling strength values across sampled sets in the APS.

Fig 7 .
Fig 7. Mechanisms for dynamic allodynia represented as shortest paths to the allodynia surface for the dynamic allodynia subcircuit.(A) Violin plots for coupling strength distributions in the 4 clusters of shortest path lengths to the dynamic allodynia surface: Cluster 1 (green, top panel), Cluster 2 (blue, 2nd panel), Cluster 3 (red, 3rd panel), and Cluster 4 (cyan, bottom panel).Black * represents the mean APS coupling strength value, colored square is mean value for the cluster.(B) Probability distribution of the shortest distances to the dynamic allodynia surface (overall profile); shading represents the contributions of each cluster to the overall profile.For instance, a bar that is 50% dark blue indicates that cluster 2 constitutes 50% of subcircuit instantiations with the corresponding distance to the allodynia surface.Only ∼2% of sampled subcircuit instantiations belong to cluster 4 with distances greater than 3.5, thus contributing minimally to distribution.(C) Parallel plot representation of the shortest paths to the allodynia surface from each sampled point in the APS.A line gives the components of the shortest displacement vector from an APS point to the allodynia surface.Lines are colored based on clusters with a color gradient for better visualization.(D) Firing-rate responses to a noisy Aβ input signal (active during t ∈ [0.2, 0.7]) for each population (columns) and for each cluster (rows).Solid lines correspond to a subcircuit instantiation with coupling strengths given by the mean values for the particular cluster, and dash-dotted lines correspond to the subcircuit instantiation with coupling strengths given by the corresponding closest point on the allodynia surface.Aβ input frequencies are chosen as the smallest value that induces allodynia for each cluster.
Fig 8 provides an illustration of the sampling algorithm.(See S2 Appendix for details on the implementation of this algorithm, a discussion of the volume-independence of the implementation, and proof that this produces a uniform in space sampling).

Fig 8 .
Fig 8. Schematic of the volume-independent sampling algorithm.

Figures
Figures 3D, 5C, and 7C show parallel plots of the shortest paths to the allodynia surface for the simple, static, and dynamic subcircuits, respectively; and are colored according to the clusters assigned by density-based scanning.The shortest paths in those figures clearly divide visually into spatially-separated clusters.

Table 1 .
[10,20]ons and resulting inequalities used to constrain the simple subcircuit and define its allowable parameter space (APS).Condition type and Condition (first 2 columns) describe the rationale behind each condition.The resulting inequality on population steady-state voltages V ss E and V ss I are given in the 3rd column.These inequalities must hold for typical non-painful f Aβ values, which we take to be f Aβ ∈[10,20]Hz.Vx,max and V x,min are the maximum and minimum limits on average voltage, respectively, Vx,rest is the resting membrane voltage and V x,thr is the voltage threshold for firing (x = E, I).

Table 2 .
Conditions that the proposed subcircuit mediating static allodynia must satisfy and the resulting inequalities on steady-state voltages.We ensure that the subcircuit exhibits these behaviors by imposing conditions (middle column) on the subcircuit exhibited in either control, I1-ablation, or I2-ablation conditions (left-most column), and is realized as an inequality (right-most column) on the steady-state voltage of a population.
To analyze sensitivity of this subcircuit to allodynia, we consider the allodynia surface S dyn , which in normalized parameter space is the set of points in (ĝ AβI1 , ĝI1E1 , ĝAβE1 , ĝE1E2 , ĝAβI2 , ĝI2E2 , ĝAβE2 ) space above which subcircuit instantiations produce E2 firing for at least some typical f Aβ (Section 4.7).Clustering points of the APS according to the directions of the shortest vectors to S , similar to Cluster 1 in the simple subcircuit.An additional similarity is that points in Cluster 1 are generally closer to S dyn compared to the other clusters (Fig7B).The shortest vectors to S dyn for Cluster 1 involve traveling only in directions associated with parameters of the I2 -E2 component of the subcircuit (Fig 7C).Specifically, to induce allodynia primarily involves decreases in ĝAβI2 leading to a dramatic reduction in the activity of I2 in response to normal Aβ input and the release of the inhibitory gate on E2.Simulations of Cluster 1 subcircuit instantiations with parameter values on S dyn show this effect (Fig 7D, top row, dash-dotted curves) where, in response to normal Aβ input, activity of the I2 population is sufficiently decreased so that E2 fires, thereby producing allodynia.The majority of points in the APS belong to Cluster 2 (blue, Fig 7B) and, as a result, the coupling strength values for points in Cluster 2 do not greatly differ from the mean or from the overall distributions of the APS values (Fig 7A).The direction of the shortest paths to S dyn involves decreasing ĝI2E2 and increasing ĝAβE2 .This corresponds to an escape of E2 from its inhibitory control, similar to the allodynia mechanism for Cluster 2 in the simple subcircuit.This is reflected in simulations of Cluster 2 subcircuit instantiations (Fig 7D) that show that on the allodynia surface (dash-dotted curves), the I2 population remains at highest activity levels in response to Aβ input in the normal range, but the E2 population is able to fire.Parameter sets in Cluster 3 differ from the mean APS values in coupling strengths associated with the I1 -E1 component of the subcircuit.Specifically, ĝAβI1 values are lower than mean values and ĝI1E1 and ĝAβE1 values are higher (Fig 7A, red).Cluster 3 points are generally quite close to the allodynia surface, almost as close as points from Cluster 1 (Fig 7B).The direction of the shortest path to the allodynia surface involves changes in all coupling parameters but the largest parameter changes reflect an allodynia mechanism involving the release of E1 from inhibitory control.In

Table 4
summarizes the model equations for each subcircuit we analyze.

Table 4 .
Equations for each subcircuit model.

Table 5 .
Model parameter values for inhibitory and excitatory populations.Parameters for each population type (first column) are divided into 3 categories (2 nd column): those that pertain to activation functions, those that give voltage cutoffs, and those that appear in the model differential equations.All these parameter values are fixed for each inhibitory or excitatory population.

Table 6 .
Inequalities on coupling strengths that define the allowable parameter space for the simple subcircuit.Inequalities on simple subcircuit coupling strength parameters that define the subcircuit APS and are obtained from the inequalities on population voltages in Table1.(UB = upper bound, LB = lower bound, abl = ablation, PI = pain inhibition)

Table 8 .
Explicit solutions of f Aβ for optimization problems in Table7.The value of f Aβ at extrema satisfying the four optimization problems in the last two rows of Table7(3rd column).All solutions (4th column) have the same form involving the Lambert W 0 function with different constants A (5th column) (UB = upper bound, abl = ablation).

Table 10 .
Hierarchical order formed by the inequalities on coupling strength parameters for each subcircuit.Order of the hierarchy formed by the inequalities for coupling strength parameters shown in Tables 6, 7,and 9 listed from lowest to highest in the hierarchy.AβI , g IE , g AβE Static g AβI1 , g AβI2 , g AβE , g I1E , g I2E Dynamic g AβI1 , g AβI2 , g I1E1 , g E1E2 , g I2E , g AβE2 For the target state to be attainable, we don't need to reach the target state for all values of f in .Instead, we can reach the target state for the value of f in that minimizes the right-hand side of the preceding equation.Thus, the circuit with the set of coupling strengths g can attain the target state if and only ifg in,y ≥ min fin∈[fin,min,fin,max] (V y,thr − V y,rest ) − n k=1 g x k y f x k ( g, f in ) f in .