Multitarget Multiscale Simulation for Pharmacological Treatment of Dystonia in Motor Cortex

A large number of physiomic pathologies can produce hyperexcitability in cortex. Depending on severity, cortical hyperexcitability may manifest clinically as a hyperkinetic movement disorder or as epilpesy. We focus here on dystonia, a movement disorder that produces involuntary muscle contractions and involves pathology in multiple brain areas including basal ganglia, thalamus, cerebellum, and sensory and motor cortices. Most research in dystonia has focused on basal ganglia, while much pharmacological treatment is provided directly at muscles to prevent contraction. Motor cortex is another potential target for therapy that exhibits pathological dynamics in dystonia, including heightened activity and altered beta oscillations. We developed a multiscale model of primary motor cortex, ranging from molecular, up to cellular, and network levels, containing 1715 compartmental model neurons with multiple ion channels and intracellular molecular dynamics. We wired the model based on electrophysiological data obtained from mouse motor cortex circuit mapping experiments. We used the model to reproduce patterns of heightened activity seen in dystonia by applying independent random variations in parameters to identify pathological parameter sets. These models demonstrated degeneracy, meaning that there were many ways of obtaining the pathological syndrome. There was no single parameter alteration which would consistently distinguish pathological from physiological dynamics. At higher dimensions in parameter space, we were able to use support vector machines to distinguish the two patterns in different regions of space and thereby trace multitarget routes from dystonic to physiological dynamics. These results suggest the use of in silico models for discovery of multitarget drug cocktails.

A large number of physiomic pathologies can produce hyperexcitability in cortex. Depending on severity, cortical hyperexcitability may manifest clinically as a hyperkinetic movement disorder or as epilpesy. We focus here on dystonia, a movement disorder that produces involuntary muscle contractions and involves pathology in multiple brain areas including basal ganglia, thalamus, cerebellum, and sensory and motor cortices. Most research in dystonia has focused on basal ganglia, while much pharmacological treatment is provided directly at muscles to prevent contraction. Motor cortex is another potential target for therapy that exhibits pathological dynamics in dystonia, including heightened activity and altered beta oscillations. We developed a multiscale model of primary motor cortex, ranging from molecular, up to cellular, and network levels, containing 1715 compartmental model neurons with multiple ion channels and intracellular molecular dynamics. We wired the model based on electrophysiological data obtained from mouse motor cortex circuit mapping experiments. We used the model to reproduce patterns of heightened activity seen in dystonia by applying independent random variations in parameters to identify pathological parameter sets. These models demonstrated degeneracy, meaning that there were many ways of obtaining the pathological syndrome. There was no single parameter alteration which would consistently distinguish pathological from physiological dynamics. At higher dimensions in parameter space, we were able to use support vector machines to distinguish the two patterns in different regions of space and thereby trace multitarget routes from dystonic to physiological dynamics. These results suggest the use of in silico models for discovery of multitarget drug cocktails.

INTRODUCTION
A large number of physiomic pathologies can produce hyperexcitability in cortex. In motor cortex, this hyperexcitability will manifest as alterations in movement and muscle tone. At the most extreme, hyperexcitability leads to a seizure with uncontrolled movement, as seen in epilepsia partialis continuans. Lesser hyperexcitability produces a variety of hyperactive movement disorders, including tics, chorea, tremor, etc, whose pathophysiology is not restricted to cortex, but involves multiple brain areas including basal ganglia, thalamus, cerebellum, and others. We focus here on dystonia, a movement disorder that produces prolonged involuntary muscle contractions (Neychev et al., 2008;Crowell et al., 2012).
The large variety of dystonias of different etiologies may present with involvement of one or several parts of the body. Pediatric causes of dystonia include cerebral palsy and are generally distinct from adult-onset cases. Common adult dystonias are torticollis, causing involuntary head turning, and movement-overuse dystonias such as writers cramp. Despite these differences, dystonias in different patient populations are primarily treated with the same therapies. While most research in dystonia has focused on basal ganglia, much pharmacological treatment is provided directly at muscles. Similarly, we propose that treatment could be targeted elsewhere in the motor pathway, here focusing on motor cortex as a potential target for therapy.
As with many other movement disorders, the dystonias generally lack a reliable biomarker and are diagnosed by semiology, the assessment of signs and symptoms. However, all dystonias feature excessive muscle activation that is associated with hyperactivity in multiple motor areas associated with movement activation. Electrophysiological studies of dystonia patients confirms a pattern of hyperactivation in cortex. Healthy individuals show low beta oscillations (∼15-20 Hz) in motor cortical local field potential (LFP). This beta is reduced in amplitude and synchrony during movement (Jasper and Penfield, 1949;Pfurtscheller and Aranibar, 1979;Crone et al., 1998;Miller et al., 2007). In dystonia patients, motor cortex shows increases in neuronal activity levels (Nobrega et al., 1995;Pratt et al., 1995), with relatively high beta amplitude and high functional connectivity at the beta frequency (Schnitzler and Gross, 2005;Jin et al., 2011). There is also excessive neural synchrony both at rest and in certain phases of movement (Toro et al., 1994;Kristeva et al., 2005;Mallet et al., 2008;Crowell et al., 2012). Some dystonias, in common with several other movement disorders, are thought to have their origin in the basal ganglia. Other dystonias, such as those associated with cerebral palsy and with movement overuse, probably have a strong cortical component. In all cases, however, the interconnections of brain motor systems makes it clear that multiple brain areas will be "in the loop" of abnormal activity. Following some primary insult or insults to a brain area, a secondarily-involved brain area will contribute further to the disorder by reacting to the alterations in input activity through its own homeostatic responses. In some cases these homeostatic changes may be compensatory so as to reduce the severity of the symptoms. However, in other cases, plasticity may actually exacerbate the abnormal movements Neychev et al., 2008;Casellato et al., 2014).
There are at least two, and perhaps more, cerebello-thalamostriato-cortical loops that play a role in movement disorders. There may also be additional contributions from still longer loops involving recurrent connections from spinal cord or from muscle spindles. One or more of these sites may have associated pathology. Regardless of the locus of primary pathology, multiple sites are potential targets where therapy could interrupt pathophysiological dynamics. Currently, brain pharmacotherapy often fails and patients are treated with botulinum toxin to partially paralyze muscles by blocking nicotinic cholinergic transmission at the affected muscle. Another treatment is deep brain stimulation using implanted electrodes. In this paper, we take two or three steps back from the level of muscle treatment by proposing interventions at the level of motor cortex.
Complex multifocal diseases may require complex multitarget treatments (Viayna et al., 2013). In the context of brain disease, multitarget therapy can hit multiple brain regions or multiple receptors in a region or both. High-level models that include many brain areas can assist in understanding how different brain areas contribute to a disorder (Sanger and Merzenich, 2000;Sanger, 2003;Hendrix and Vitek, 2012;Kerr et al., 2013). However, these models typically lack biological detail, making them unsuitable for assessing the impact of specific pharmacological manipulations. Detailed models are not yet elaborated to the point of handling multiple brain areas but do provide the details needed to assess pharmacological intervention more directly.
Single agent treatments for disease are traditionally tested in vitro or in vivo. As noted above, single agent treatments for dystonia have not had much success. There is, however, the potential for success with multitarget drug cocktails that could target multiple locations in the brain, or multiple drug receptor targets at a single location, or both (Delnooz and van de Warrenburg, 2012). Due to combinatorial explosion, evaluating combinations of drugs in different dosages in this way can not be readily done in tissue and is most feasible in silico (Viceconti et al., 2008;Kohl and Noble, 2009;Lytton et al., 2014;Action, 2016;Viceconti et al., 2016). In this study, we use our detailed multiscale model of primary motor cortex to assess potential multitarget pharmacological therapies for treatment of dystonia. The model contains 6 cortical layers with multiple classes of excitatory and inhibitory neurons, using wiring based on mouse microconnectomic data (Shipp, 2005;Weiler et al., 2008;Kiritani et al., 2012;Hooks et al., 2013). Excitatory neurons contain intracellular molecular mechanisms that contribute to persistent activity and hyperexcitability . These mechanisms include endoplasmic reticulum associated calcium stores released by activation of IP 3 Rs, and ryanodine receptors, both with affinity for caffeine, an agent that can exacerbate dystonia symptoms (Richter and Hamann, 2001). Plasma membrane calcium, sodium, and potassium channels also contribute to cellular excitability.
Since our model does not include spinal cord and muscle, we defined dystonia pathology as a state of cortical hyperactivation characterized by increased beta oscillations with excessive and hypersynchronous firing in layer 5 corticospinal neurons. These layer 5 neurons project downward to brainstem and spinal cord, and their sustained firing would lead to the increased muscle contractions of dystonia. We distinguished the hyperexcitability of dystonia from the still greater hyperexcitability of a seizure by excluding simulations that showed higher levels of activity with higher frequency oscillation and a strong tendency to "latch-up" through multicell depolarization blockade (Lytton and Omurtag, 2007). Classification in 11-dimensional space demonstrated that we could identify different regions in parameter space for these different states-baseline normal, dystonia, epileptiform-and predict pharmacological combinations that would lead from pathology back to the physiological activity state. As in our previous investigations of epilepsy (Lytton and Omurtag, 2007), we found multiple parameter combinations that were consistent with the pathological state, as well as multiple parameter combinations to produce our baseline physiological state. Such parameter degeneracy is typical of complex neural systems (Edelman and Gally, 2001;Golowasch et al., 2002).

MATERIALS AND METHODS
Network simulations consisted of 1715 reduced compartmental cell models with single compartments for inhibitory cells and five compartments for pyramidal cells, arrayed by layer with connectivity taken from experimental results on motor cortex (Weiler et al., 2008; Figures 1A,B). Parallel-conductance electrophysiological simulation in the pyramidal cells was complemented by chemophysiological simulation focused on Ca 2+ handling, based on our prior models (Neymotin et al., 2015; Figure 1C).
We briefly describe the scales of the multiscale model from smaller to larger in the following sections ( Table 1). For more details, readers are referred to our previous papers (Neymotin et al., 2015.

Intracellular Molecular Scale
Our Ca 2+ dynamics ( Figure 1C), are based on . We modeled a one-dimensional RxD system of intracellular neuronal Ca 2+ signaling in all compartments of neocortical pyramidal (PYR) neurons. Within each compartment, we modeled cytosolic and endoplasmic reticulum (ER) sub-compartments by using a fractional volume for each.
IP 3 was produced through a reaction sequence initiated by glutamate binding to the metabotropic glutamate receptor (mGluR), based on a reaction scheme developed by Ashhad and Narayanan (2013) (ModelDB #150551). IP 3 diffused outward from the synapse location and decayed following first-order kinetics (τ IP 3 of 1 s). Baseline mGluR synaptic weight was normalized to represent the increase in the amount of glutamate bound to mGluR. Extracellular glutamate did not diffuse but was represented by a local Glu value that was incremented in response to an event delivered due to a presynaptic spike. Glu showed bind/unbind kinetics on mGluR and was eliminated by first-order degradation (lower left of Figure 1C).
The ER Ca 2+ model involves IP 3 receptors (IP 3 Rs), ryanodine receptors (RYR) (Sneyd et al., 2003), SERCA pumps, and a Ca 2+ leak. IP 3 R dynamics involved a slow Ca 2+ inactivation binding site state (De Young and Keizer, 1992;Li and Rinzel, 1994). The SERCA pump is a pump rather than a channel and so is modeled with Hill-type dynamics. Calcium buffering followed CaB where B is diffusible buffer with diffusion coefficients D = 0.043 µm 2 /ms for both B and CaB, about half the rate of Ca 2+ diffusion (Anwar et al., 2012). Calcium extrusion across the plasma membrane was modeled by first-order decay with τ ex = 5 ms.

Synapses
AMPA/NMDA synapses were modeled by standard NEURON double-exponential mechanisms ( Table 2). All excitatory projections were mixed AMPA (rise,decay τ : 0.05, 5.3 ms) and NMDA (rise,decay τ : 15, 150 ms). NMDARs were scaled by 1/(1 + 0.28 · Mg · exp(−0.062 · V)); Mg = 1mM (Jahr and Stevens, 1990). 13% of I NMDA was carried by Ca 2+ (Spruston et al., 1995). AMPA and NMDA receptors had reversal potentials of 0 mV. Inhibitory synapses were mediated by GABA A and GABA B receptors. GABA A synapses were modeled with a doubleexponential mechanism. The GABA B synapse had second messenger connectivity to a G protein-coupled inwardlyrectifying potassium channel (GIRK). LTS cells connected to apical dendrites of PYR cells using GABA A receptors (GABA A R; rise,decay τ : 0.2, 20 ms) and GABA B receptors (GABA B R) and onto somata of FS and other LTS cells with GABA A Rs (rise,decay τ : 20, 40 ms). GABA A Rs had reversal potentials of −80 mV, and GABA B Rs −95 mV. GABA B Rs provide longer-lasting activation compared to GABA A Rs.
Voltage-gated ionic current models were based on prior models of our own and others (McCormick and Huguenard, and second-(and third-etc) messenger modulation (red back-beginning arrows). We distinguish membrane-associated ionotropic and metabotropic receptors and ion channels involved in reaction schemes in red (in reality, it is likely that almost every membrane-bound protein is modulated). External events are represented by yellow lightning bolts-there is no extracellular diffusion; the only extracellular reaction is glutamate binding, unbinding, and degradation on mGluR1 after an event. Ca 2+ is shown redundantly in blue-note that there is only one Ca 2+ pool for extracellular, 1 pool for cytoplasmic, and 1 pool for ER (PLC, phospholipase C; DAG, diacyl-glycerol; cAMP, cyclic adenosine monophosphate; PIP 2 , phosphatidylinositol 4,5-bisphosphate). Adapted from Figure 1  1992; Migliore et al., 2004;Stacey et al., 2009;Neymotin et al., 2011bNeymotin et al., ,a, 2013. Voltage sensitive channels generally followed the Hodgkin-Huxley parameterization, wherebyẋ = (x ∞ − x)/τ x (x = m for activation particle and h for inactivation particle). Steady-state x ∞ and time constant τ x are either related to channel opening α(V) and closing kinetics β(V) as x ∞ = α/(α + β), τ x = 1/(α + β), or are directly parameterized: x ∞ (V), τ x (V). Kinetics for channels were scaled by Q 10 from an experimental temperature (where available) to simulation temperature of 37 • C. Q 10 = 3 was used when no experimental value was available. All cells contained leak current, transient sodium current I Na , and delayed rectifier current I K−DR , to allow for action potential generation. Additionally, PYR cells contained in all compartments I K−A , I K−M providing firing-rate adaptation (McCormick et al., 1993). Pyramidal cells also had I h , voltagegated calcium channels (VGCCs) in all compartments: I L , I T , I N (Kay and Wong, 1987;McCormick and Huguenard, 1992;Safiulina et al., 2010;Neymotin et al., 2015), and SK and BK    (Hagiwara and Irisawa, 1989;Schwindt et al., 1992;Chen et al., 2001;Wang et al., 2002;Robinson and Siegelbaum, 2003). The hyperpolarizationactivated HCN current I h used in pyramidal cells was modeled with second messenger and calcium dependence taken from Winograd et al. (2008) (ModelDB #113997), and modified to provide the faster voltage-sensitivity time constants found in cortex (Harnett et al., 2015), and provides PYR cells longerlasting firing activity via augmentation of the HCN channel conductance. The mechanism for Ca 2+ regulation of HCN channels in PYR cells in Winograd et al. (2008) is modeled empirically in order to produce the relationship between cytosolic Ca 2+ levels and I h activation without simulating the details of Ca 2+ effects on adenyl cyclase (see schematic for HCN chan in Figure 1C).
g h was 0.0025 S/cm 2 in PYR soma, basal dendrites and exponentially-increasing in apical dendrites with distance from soma with 325 µm space constant, hence e-fold augmented at 325 microns as described by Kole et al. (2006). Apical dendrite I K−DR , I K−A , I K−M density also increased with the same length constant, based on data showing HCN and Kv channel colocalization (Harnett et al., 2015(Harnett et al., , 2013.

Network Scale
The network consisted of 1715 cells ( Table 4). The network contained 157,507 synapses for an overall connection density of ∼5% (see Table 6). PYR cells synapsed onto each-other's dendrites. PYR-to-PYR synaptic locations on the dendrite were randomized between basal and apical compartments (Markram et al., 1997).

PYR cells synapsed onto somata of FS and LTS cells (single-compartment models).
Neuronal populations were arranged by cortical layer based on our prior models (Neymotin et al., 2011a,c;Chadderdon et al., 2014;Neymotin et al., 2016), with additional data from direct measurements from mouse motor cortex (Shipp, 2005;Weiler et al., 2008;Kiritani et al., 2012;Hooks et al., 2013), and recent experiments which demonstrate a thin layer 4 in mouse motor cortex (Yamawaki et al., 2014). The network consisted of 13 populations of 3 excitatory cell types, intratelencephalic (IT), pyramidal-tract (PT), and corticothalamic (CT), and 2 inhibitory cell types, low-threshold spiking (LTS) and fastspiking (FS). These were distributed across cortical layers 2/3, 4, 5a, 5b, and 6 (Harris and Shepherd, 2015), with cell numbers for each population based on estimated cell densities and volume ( Table 4). The volume of each population was calculated assuming a horizontal area (x and z dimensions) of 120 × 120 µm, and a layer-dependent cortical depth range (y dimension).
Connection probabilities p ij (Tables 5, 6) of presynaptic excitatory populations were dependent on pre-and pothstsynaptic type and layer. For presynaptic inhibitory populations, connection probabilities inversely scaled based on distance p ij =p ij · exp(− (dx 2 + dz 2 )/15), in x, z plane perpendicular to the y-direction of layering. Connection probabilities and weights are based on data from rodent motor cortex mapping (Weiler et al., 2008;Lefort et al., 2009;Anderson et al., 2010;Fino and Yuste, 2011;Apicella et al., 2012;Kiritani et al., 2012). Individual neurons were placed randomly with uniform distribution. Weights from E cells displayed in Table 6 are for the AMPA synapses, with colocalized NMDA weights at 10% of AMPA weights. Synaptic delays were randomized between 1.8 and 5 ms with additional delay based on distance.  Property Description E to E p ij , w ij dependent on pre-/post-synaptic cell type/layer E to I p ij , w ij dependent on pre-synaptic cell layer, and post-synaptic cell type/layer I to E/I p ij decreases exponentially with x,z plane distance between pre-/post-synaptic neurons; fixed w ij All delays Randomized 1.8-5 ms with additional delay based on distance p ij denotes probability of connection between type i and j; w ij denotes weight. Parameters by pre-and post-synaptic type in Table 6.
Background activity was simulated by excitatory and inhibitory synaptic inputs following a Poisson process, sent to all cells, representing ongoing drive from other cortical areas and other inputs. These inputs were selected to maintain low-frequency firing of neurons within the model, which would not fire otherwise, due to small network size and the requirement for multiple synaptic inputs to trigger a postsynaptic spike (Neymotin et al., 2011a). The strength of these background inputs was not based on the full source of inputs from all upstream brain areas, since these inputs are not completely mapped.

Simulation Variations
We ran over 5800 simulations, randomly varying each of the following parameters using an independent normal distribution: 1. E neuron mGluR density (mGluR); 2. E neuron ER RYR density (RYR); 3. E and I neuron HCN channel density; 4. E and I neuron fast Na + channel density (Na f ); 5. E neuron K dr channel density; 6. E neuron K a channel density; 7. E neuron K D channel density; 8. E neuron K M channel density; 9. E neuron SK calcium-activated potassium channel density; 10. E neuron BK calcium-activated potassium channel density; 11. E and LTS neuron voltage-gated calcium channel (VGCC) density.
Means and standard deviations differed for the different parameters and were selected to allow each simulation to maintain activity. A subset of the simulations was used for the analyses described (Table 7).
We ran simulations with initial calcium concentration in the ER set to 1.25 mM (Bygrave and Benedetti, 1996), to mimic the effects of ER calcium priming via prior excitatory synaptic stimulation (Ross et al., 2005;Hong and Ross, 2007;Fitzpatrick et al., 2009;Neymotin et al., 2016).
We categorized the simulations into distinct groups by noting major differences in activity across parameter sets (Table 8).
From the full set of 5867 simulations, 1505 did not display any firing due to random variations in ion channel densities which led to low neuronal activity ( Table 7). The remaining 4341 simulations were Active due to higher neuronal activity, e.g., partially caused by the higher average Na f density in these simulations. Of these 4341 Active simulations, 1077 exhibited epileptic latch-up dynamics-periods of intense activity which led to depolarization blockade (Na + channel inactivation; Lytton and Omurtag, 2007). These periods where neurons did not fire lasted 200-300 ms (gaps between E5P spikes: E5P gap in Table 8). We categorized the top and bottom 2nd percentile of the Active/non-latch-up simulations ranked by E5P firing rate into dystonia (n = 65) and physiological (n = 65) sub-sets. We used E5P firing rate as a criterion for dystonia classification because E5P neurons project downward to brainstem spinal cord, and sustained overactive E5P firing can lead to the tonic muscle contractions symptomatic of dystonia.

Data Analysis
We formed multiunit activity (MUA) time-series, which count the number of spikes in each bin (10 ms) for a given population. To calculate neuronal population rhythms, we took the power spectral density (PSD) of the mean-subtracted MUA time-series; we then calculated the peak frequencies and amplitudes in the PSD. We used the average Kendall's τ nonparametric rank-correlation coefficient (Kendall, 1938;Knight,

1966) between pairs of neuron binned spike train time-series
for calculating the synchronization of populations of neurons (denoted population-synchrony). Kendall's τ non-parametric rank correlation, defined as: , is used with these data. Kendall's τ is a normalized difference between concordant (n c ) and discordant pairs (n d ); ties are taken into account by the normalizing term, 1 2 n(n − 1) , which represents the total number of ordered pairs in the time-series. We used the Python scikit-learn library (Pedregosa et al., 2011) for performing principal component analysis (PCA) and support-vector machine (SVM) classification (Cortes and Vapnik, 1995;Orrù et al., 2012). Dystonia and physiological simulation classes were characterized on the basis of layer 5 corticospinal pyramidal neuron (E5P) firing rates. The clearest examples of both classes (bottom/top 2nd percentiles as physiological/dystonia classes) were used for the majority of the analyses described in the Results (Figures 3-8). The NuSVC variant of SVMs was used to classify physiological and dystonia simulations and to find which simulation parameters contributed the most to classification accuracy. SVM inputs TABLE 7 | Parameter ranges (average ± standard deviation) for all simulations (n = 5867), active simulations (n = 4341), latch-up simulations (n = 1077), active/non-Latch-up simulations (n = 3264), physiological simulations (n = 65), and dystonia simulations (n = 65).    were vectors consisting of normalized parameter values. Each of these input vectors was labeled into either of two distinct binary classes: physiological (0) or dystonia (1). SVM parameters, including kernel type (linear, polynomial, radial-basis function), γ , tolerance, ν, and polynomial degree were selected using a grid search with N-fold cross validation run 10 times for each combination of parameters. SVM classification accuracy surpassed the accuracy of other machine learning methods, including logistic regression (not shown). Figures were drawn with Matplotlib (Hunter, 2007).

Simulation Overview
We ran over 5800 network simulations, randomizing 11 ion channel/receptor densities independently. A typical 2 s simulation took ∼3 min using 24 cores on Linux with parallel NEURON. After running simulations, we calculated neuronal population firing rates, synchronization, and power spectra.

Characterization of Dystonia Pathophysiology
Simulations were grouped into physiological and pathological based on differences in firing patterns (Table 8, Figure 2). 1505 of 5867 simulations produced no activity. The remaining simulations were characterized as physiological or pathological. Pathological simulations showed increased activity. High activity alternating with latch-up condition was defined as an epileptiform simulation with periods of >200 ms of depolarization blockade across multiple cells (1077 simulations Figure 2A, but displayed either high or low FV similarity which overlapped with the range of values displayed by the physiological simulations. Epileptiform simulations had intermediate average E5P rates due to high activity alternating with periods of quiescence caused by depolarization blockade. Across simulation types, higher E5P firing increased the excitatory drive to I5 neurons, causing increased I5 neuron firing ( Table 8). Higher I5 and E5P neuron firing then caused higher E5P synchronization via recurrent E5P excitation and feedback inhibition ( Figure 2B). Stronger E5P and I5 interactions then increased beta rhythm amplitude (Figure 2B), however with substantial variability. Peak oscillatory frequency was held relatively stable across simulations (Table 8). Physiological and epileptiform simulations had lower overall E5P synchrony and beta power compared to the dystonia simulations, which occupied the upper-right quadrant of Figure 2B. E5P FV similarity showed temporal recurrences which further distinguished the three simulation types (Figure 2C). The physiological simulation showed intermediate self-similarity (0.17) due to sparse firing of different subsets of pyramidal cells at different times. In contrast, the dystonia simulation firing patterns showed strong self-similarity (0.56) and recurrence over time (recurring orange/red blobs in Figure 2C), indicating stereotyped dynamics. The example epileptiform simulation showed relatively weak self-similarity (0.16) due to its two distinct firing patterns: high E5P synchrony alternating with E5P silence produced by depolarization blockade. Epileptiform and dystonia simulations showed a brief period of high similarity when the epileptiform simulation showed strong oscillations during the initial period. There was weak similarity between epileptiform and physiological (0.12) and dystonia and physiological (0.22) simulations, indicating that both pathological dynamics were distinct from the physiological.
E5P neurons in a representative physiological model fired sparsely with low synchrony (population-synchrony = 0.01; Figures 3A,D; Supplementary Figure 1 has all physiological rasters), with multiple downstream effects. Low excitatory drive from E5P to I5 and I5L neurons caused them to fire slowly. This low L5 inhibition allowed E5a neurons to fire quickly. The weak E5P and L5 interneuron interactions produced only weak beta rhythms which were confined to layer 5 (Figure 4A).
In a representative dystonia simulation, E5P neurons had sustained, synchronous, rapid firing (Figures 3B,D; Supplementary Figure 2 shows all dystonia simulation rasters). This promoted strong, continuous layer 5 interneuron activation. The L5 interneurons then suppressed E5a intratelencephalic  neurons, which fired at reduced rates. In contrast, E5b firing increased with the faster E5P firing, due to excitation spreading in the network. The relatively high recurrent connectivity (18% density) and strong synaptic weights between E5P neurons allowed the E5P neurons to remain activated despite strong feedback inhibition. The strong feedback inhibition also activated the E5P HCN channels, which produced rebound excitation. The strong E5P activation coupled with the feedback inhibition also enabled E5P neurons to synchronize (population-synchrony = 0.83; vertical stripes in Figure 3B) at a strong beta rhythm (∼20 Hz; Figure 4B). These synchronous beta rhythms also spread to other populations and layers (E2, I5, I5L, E5b, and E6).
Epileptiform simulation also displayed strong intermittent beta oscillations and strong synchrony (population-synchrony = 0.05; Figures 3C, 4C), but this activity alternated with lengthy periods (200-300 ms) where E neurons were not firing due to depolarization blockade. Even with these periods of depolarization blockade, most E neurons fired at higher average rate than in the physiological simulations ( Figure 3D). Such increased synchrony with high excitatory cell activity is seen in epilepsy patients (Meisel et al., 2015). In contrast to the dystonia simulations, the synchronous periods of epileptiform oscillations were largely confined to layer 5 and did not spread to other layers.

Need for Multitarget Approach
No individual parameter determined physiological vs. dystoniadynamical-condition in the network (Figure 5). Therefore, no single parameter adjustment would routinely provide an effective "treatment" that would reliably restore physiological activity in most pathological models. We therefore went on to explore whether multitarget manipulation would be able to find such treatment routes.
Although no single parameter could predict physiological vs. pathological dynamics, the outliers of certain individual parameters were predictive. At the pathological margin, simulations had parameters which are expected to produce high activity: high Na + or Ca 2+ channels promoting inward currents, high HCN channel densities providing high resting membrane potential (RMP), and low K + channel densities again producing depolarization and reduced repolarization with spiking.
Further evidence for lack of predictability of dynamics based on parameters, comes from viewing the parameters in all 11 dimensions organized into 2 classes by dynamics. The parameter space showed substantial heterogeneity in the patterns producing pathology (Figure 6A), with weak intra-class clustering ( Figure 6B). Correlation between parameter vectors of each simulation averaged 0.06 for physiological simulations, 0.07 for pathological simulations, with weak -0.05 anticorrelation between pathological and physiological simulations. The low correlations in both the physiological simulations (lower-left corner of Figure 6B) and the pathological simulations (upperright corner of Figure 6B) demonstrate that there is widespread degeneracy in the parameter sets that produce either the physiological or pathological states. Some of this degeneracy is unsurprising: for example K + channels with similar time courses of activation can substitute for one another to some extent. Other degeneracy is more complex and involves higher-order dynamic compensation.

High Dimensional Separation of Physiological and Pathological Parameters
Because of the difficulty of separating pathological from physiological with these high dimensional parameter sets, we used a SVM classification to create a separation (termed a maximum margin hyperplane) separating parameter sets producing physiological dynamics from parameter sets producing pathological dynamics. We started by training SVMs using only two parameters in combination (Figure 7). In order to test the efficiency of this separation, we separated out our target sets (physiological vs. pathological) into two subsets of each to serve as training and testing sets to evaluate the adequacy of the separation. By trying various random training and testing sets we got a mean and standard error for each case. Many two-parameter predictions were below chance (0.5) indicating that the SVM could not separate physiological from pathological based on that parameter pair. Two-parameter SVMs could accurately classify when the parameter pair included Na f density-the strongest predictor of excitability. The best prediction came with high Na f and low K dr . Logistic regression methods were also tried to perform this two-dimensional separation but did not perform as well as SVM.
Going beyond 2 parameters, SVM classification accuracy increased regularly with the number of parameters used  (Figure 8), suggesting that a multi-target drug approach beyond two targets might produce greater effect. Moving to higher and higher dimensional spaces, we checked all possible parameter combinations at each dimensionality. In Figure 8, we report the parameter combination that was most predictive-e.g., at 6 dimensions we report just one of the 462 combinations of six from 11 parameters. Looking at the red blocks below, we can identify that the six dimensions that provide best prediction are Na f , four of the K + channels, and VGCC. Predictability increases up through six parameters, then plateaus, and then falls off due to the extreme sparseness of data. This sparseness was due to the so-called curse of dimensionality: given a constant number of data points n, the density falls off #bin-fold with each increase in dimension, where #bin is the binning of the space in one dimension. Because of this, any high-dimensional method will tend to underestimate predictive strength given a limited amount of data (Bishop, 2006;Noble, 2006). This multi-target SVM approach revealed the parameters that had the highest contribution to producing or preventing dystonia. Na f density was the most predictive parameter across all numbers of parameters used (horizontal red stripe at top of Figure 8B), as had been also shown using 2 dimensions alone (Figure 7). Again confirming the 2-dimensional result, the next most predictive parameters was K dr . Following these came K a , K d , BK, SK, and VGCC densities which also significantly contributed to accurate predictions, due to their strong influence on E neuron excitability. mGLUR, RYR, and K m densities showed lesser contributions. FIGURE 7 | Support vector machine classification accuracy of pathological vs. physiological simulations using two parameter values has high levels for certain parameter combinations (e.g., including Na f channel density) but overall accuracy is often below chance (0.5). (A) Accuracy as a function of specific parameter combinations [indicated at same horizontal location in (B) (Red indicates parameter (param) was used for classification; blue indicates the parameter was not used)] (solid line: mean cross-validation accuracy (n = 10); dotted line: standard error of cross-validation accuracies).
Increasing the percentile cutoffs for categorizing physiological from pathological simulations from the 2nd percentile to 7th percentile decreased prediction accuracy but still demonstrated the value of multitarget changes (Figure 9). The left column shows the same result as Figure 8: accuracy increased (colormap) as one goes from fewer to more parameters (bottom to top). By including more exemplars on both the physiological and pathological sides, we moved away from the best exemplars and obtained less distinction between the two parameter sets. However, at all percentiles, there was an initial increase in classification accuracy with continued increase up to or beyond 3 parameters. This increase then declined as the number of parameters increased further due to the aforementioned sparseness at high dimensionality.

DISCUSSION
We developed a multiscale model of primary motor cortex to explore multitarget pharmacological therapies for dystonia. We searched parameter space-channel and receptor densitiesto create a set of models to contrast dystonia dynamics with physiological dynamics. Dystonia simulations displayed high excitability and synchrony in layer 5 corticospinal neurons (E5P), and strong beta oscillations which spread between cortical layers (Figures 3B, 4B). Dystonia simulations could be distinguished from epileptiform simulations primarily by the presence of periods of latch-up with depolarization blockade in the epileptiform simulations. Physiological simulations had low excitability, asynchronous firing, and weak beta oscillations (Figures 3C, 4C). Attempts to use high-dimensional visualization techniques to find potential therapeutic directions in the parameter space were limited by the solution degeneracy in the 11-dimensional parameter space with scattered parameter vectors with low correlation (Figure 6). We therefore turned to a SVM classification to identify hyperplanes in high-dimensional space that would separate the two populations. As expected, the major spike generating channels, Na f and K dr were the primary determinants of excitability, followed by additional potassium channels and calcium channels. We did not assess pharmacological effects on synapses, which would be useful to do in the future.

Biological Degeneracy and Personalized Therapy
Degeneracy of mechanism is a major theme in biology (Edelman and Gally, 2001), meaning that there are many different ways that a biological system can produce a particular shape in the case of an immunoglobulin, or a particular dynamics in the case of a neural system. Such degeneracy has been shown directly in the stomatogastric ganglion of lobster, where a particular cell type produces its stereotyped dynamics using many different combinations of ion channel densities (Golowasch et al., 2002). Associated with this degeneracy is failure of averagingaveraging across parameter sets that produce the dynamics gives a set of parameter values that do not produce the same dynamics.
In the context of brain physiology, this means that we can expect that individuals differ in the details of how their FIGURE 9 | SVM classification accuracy increases with more parameters then decreases due to "curse of dimensionality"-sparseness of parameter vectors relative to dimension. Best classification accuracy from all combinations of y parameters (params) using bottom/top SPI firing rate percentiles on x-axis. motor cortex produces oscillations and contributes to movement. Similarly, we can expect that individuals differ in the details of their pathology. From a pharmacological perspective this argues that we may see greater benefit from personalized medicineidentifying the high-dimensional complex of pathological parameters in a particular patient in order to treat them with their own individualized cocktail of multitarget drugs to produce a dynamics that falls somewhere in the physiological regime. To this might also be added complementary individualized, perhaps multi-locus, brain stimulation (Kerr et al., 2012;Song et al., 2013;Chadderdon et al., 2014;Hiscott, 2014;Nelson and Tepe, 2014;Dura-Bernal et al., 2016). Such a personalized approach would require much more intensive, and more costly, diagnostic procedures of a type that is currently only used by epilepsy surgery centers, which typically require invasive methods to perform their diagnostic tests.
Due to the degeneracy, parameter averaging failed in our dataset-using the average of all parameters sets that produce pathological simulations does not give a pathological simulation. However, the ability of the SVM method to separate pathological from physiological populations in high dimensional parameter space does suggest that there may be some benefit to pushing all patients in that direction through a multitarget pharamacological cocktail. In future studies, we plan to test this explicitly in the simulations in order to determine what percentage improve, what percentage worsen and what percentage remain essentially unchanged with such an average treatment. This assessment will require a larger number of simulated patients than we have thus far accumulated.

Multilocus, Multitarget, Multiscale Approaches for Treating Dystonia
In general, single target pharmacology has not been effective in dystonia (Fahn, 1987). As in other complex diseases, many of the treatments for dystonia have highly variable effectiveness and must be used at high doses that produce side-effects (Jankovic, 2006). For these reasons, botulinum toxin, targeting the final endpoint -the muscle movement-is commonly used as a treatment (Jankovic, 2006;Sanger et al., 2007;Bragg and Sharma, 2014). Deep-brain stimulation, an invasive procedure, is also used to partially restore normal brain dynamics (Tarsy, 2007;Johnson et al., 2008;Air et al., 2011;Bhanpuri et al., 2014).
Multilocus, multitarget approaches may be particularly useful in movement disorders because movement produces coordination by utilizing coordination among multiple brain areas including the basal ganglia, thalamus, cerebellum, sensory, and motor cortices (Neychev et al., 2008;Crowell et al., 2012;Delnooz and van de Warrenburg, 2012). Pathology within any one region, or disturbances in communication between any of the regions can potentially lead to disorders. To begin to address these multiple challenges, we focused our computer modeling here on a multiscale model of motor cortex and multitarget pharmacology based in this one area. In the future, this model will be expanded to encompass more areas and will include synaptic receptor targets in each area.

AUTHOR CONTRIBUTIONS
All authors listed, have made substantial, direct and intellectual contribution to the work, and approved it for publication.

ACKNOWLEDGMENTS
The authors would like to thank Ben Suter and Gordon MG Shepherd (Northwestern University) for help with the model; Tom Morse (Yale) for ModelDB support; R.A.N. for helpful suggestions. The authors declare no competing financial interests. Research supported by NIH grant R01 MH086638, NIH grant U01 EB017695, NIH grant R01 NS064046, NIH grant R01 DC012947. The NIH had no role in study design; in the collection, analysis and interpretation of data; in the writing of the report; and in the decision to submit the article for publication.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphar. 2016.00157