Original Research ARTICLE
Multitarget Multiscale Simulation for Pharmacological Treatment of Dystonia in Motor Cortex
- 1Department Physiology and Pharmacology, SUNY Downstate Medical Center, State University of New York, Brooklyn, NY, USA
- 2Department Neuroscience, Yale University School of Medicine, New Haven, CT, USA
- 3Nathan S. Kline Institute for Psychiatric Research, Orangeburg, NY, USA
- 4Department Biomedical Engineering, University of Southern California, Los Angeles, CA, USA
- 5Division Neurology, Child Neurology and Movement Disorders, Children's Hospital Los Angeles, Los Angeles, CA, USA
- 6Department Neurology, SUNY Downstate Medical Center, Brooklyn, NY, USA
- 7Department Neurology, Kings County Hospital Center, Brooklyn, NY, USA
- 8The Robert F. Furchgott Center for Neural and Behavioral Science, Brooklyn, NY, US
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. 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 (Sanger et al., 2003; Neychev et al., 2008; Casellato et al., 2014).
There are at least two, and perhaps more, cerebello-thalamo-striato-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 (Neymotin et al., 2016). These mechanisms include endoplasmic reticulum associated calcium stores released by activation of IP3R s, 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).
2. 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 Ca2+ handling, based on our prior models (Neymotin et al., 2015, 2016; Figure 1C).
Figure 1. Model schematics. (A,B) Motor cortex architecture. Circles represent neuronal populations (red: excitatory cells; green: fast-spiking interneurons; blue: low-threshold firing interneurons). Circle size denotes number of cells in population. Lines (with arrows) indicate connections between the populations. Thickness of lines proportional to synaptic weights. E cells synapse with AMPAR/NMDARs; I cells synapse with GABAA R / GABAB Rs. Circles with self-connects denotes recurrent connectivity. (A) Excitatory connections. E5P projects to spinal cord (not modeled). (B) Inhibitory connections. (C) Chemical signaling in pyramidal cells showing fluxes (black arrows) 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. Ca2+ is shown redundantly in blue—note that there is only one Ca2+ pool for extracellular, 1 pool for cytoplasmic, and 1 pool for ER (PLC, phospholipase C; DAG, diacyl-glycerol; cAMP, cyclic adenosine monophosphate; PIP2, phosphatidylinositol 4,5-bisphosphate). Adapted from Figure 1 of Neymotin et al. (2016).
Simulations were run in the NEURON (version 7.4) simulation environment (Carnevale and Hines, 2006) utilizing the reaction-diffusion (RxD) Python module (McDougal et al., 2013a,b) and NMODL (Hines and Carnevale, 2000). Two seconds of simulation time took ~3 min using 24 nodes on a Linux cluster with parallel NEURON, run with a fixed time-step of 0.1 ms. The full model is available on ModelDB (https://senselab.med.yale.edu/ModelDB/ShowModel.cshtml?model=189154).
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, 2016).
2.1. Intracellular Molecular Scale
Our Ca2+ dynamics (Figure 1C), are based on (Neymotin et al., 2016). We modeled a one-dimensional RxD system of intracellular neuronal Ca2+ 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.
IP3 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). IP3 diffused outward from the synapse location and decayed following first-order kinetics (τIP3 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 Ca2+ model involves IP3 receptors (IP3R s), ryanodine receptors (RYR) (Sneyd et al., 2003), SERCA pumps, and a Ca2+ leak. IP3R dynamics involved a slow Ca2+ 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 where B is diffusible buffer with diffusion coefficients D = 0.043 μm2∕ms for both B and CaB, about half the rate of Ca2+ diffusion (Anwar et al., 2012). Calcium extrusion across the plasma membrane was modeled by first-order decay with τex = 5 ms.
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 INMDA was carried by Ca2+ (Spruston et al., 1995). AMPA and NMDA receptors had reversal potentials of 0 mV.
Inhibitory synapses were mediated by GABAA and GABAB receptors. GABAA synapses were modeled with a double-exponential mechanism. The GABAB synapse had second messenger connectivity to a G protein-coupled inwardly-rectifying potassium channel (GIRK). LTS cells connected to apical dendrites of PYR cells using GABAA receptors (GABAA R; rise,decay τ: 0.2, 20 ms) and GABAB receptors (GABAB R) and onto somata of FS and other LTS cells with GABAA Rs (rise,decay τ: 20, 40 ms). GABAA Rs had reversal potentials of −80 mV, and GABAB Rs −95 mV. GABAB Rs provide longer-lasting activation compared to GABAA Rs.
2.3. Cell Scale
The network consisted of pyramidal cells (PYR; 3 apical dendrite compartments, 1 basal dendrite compartment, 1 somatic compartment), fast spiking soma-targeting interneurons (FS; one compartment), and dendrite-targeting low-threshold spiking interneurons (LTS; one compartment; Wang and Buzsaki, 1996; Wang, 2002; Monyer and Markram, 2004; Bartos et al., 2007; Neymotin et al., 2011a,b; Tables 3, 4). Reaction-diffusion mechanisms (Ca2+,IP3,buffer) were restricted to the PYR cells in this network. Properties of pyramidal neurons in the different layers were identical except for apical dendrite length which is longer in deep pyramidal neurons than in superficial (Hay et al., 2011; Castro-Alamancos, 2013): 900 μm in L5-6; 450 μm in L2/3 and L4.
Table 4. Network Population, including normalized and nominal cortical depth range (ynormRange, yRange, neuron density, and number of cells).
Voltage-gated ionic current models were based on prior models of our own and others (McCormick and Huguenard, 1992; Migliore et al., 2004; Stacey et al., 2009; Neymotin et al., 2011b,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 Q10 from an experimental temperature (where available) to simulation temperature of 37°C. Q10 = 3 was used when no experimental value was available. All cells contained leak current, transient sodium current INa, and delayed rectifier current IK−DR, to allow for action potential generation. Additionally, PYR cells contained in all compartments IK−A, IK−M providing firing-rate adaptation (McCormick et al., 1993). Pyramidal cells also had Ih, voltage-gated calcium channels (VGCCs) in all compartments: IL, IT, IN (Kay and Wong, 1987; McCormick and Huguenard, 1992; Safiulina et al., 2010; Neymotin et al., 2015), and SK and BK calcium-activated potassium currents (IKCa). LTS cells contained IL, non-Ca2+-dependent Ih, SK, and Ca2+ decay.
HCN channels in different cell types have somewhat different voltage dependence and different kinetics (Hagiwara and Irisawa, 1989; Schwindt et al., 1992; Chen et al., 2001; Wang et al., 2002; Robinson and Siegelbaum, 2003). The hyperpolarization-activated HCN current Ih 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 longer-lasting firing activity via augmentation of the HCN channel conductance. The mechanism for Ca2+ regulation of HCN channels in PYR cells in Winograd et al. (2008) is modeled empirically in order to produce the relationship between cytosolic Ca2+ levels and Ih activation without simulating the details of Ca2+ effects on adenyl cyclase (see schematic for HCN chan in Figure 1C).
gh was 0.0025 S/cm2 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 IK−DR, IK−A, IK−M density also increased with the same length constant, based on data showing HCN and Kv channel colocalization (Harnett et al., 2015, 2013).
2.4. 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 fast-spiking (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 pij (Tables 5, 6) of presynaptic excitatory populations were dependent on pre- and pothst-synaptic type and layer. For presynaptic inhibitory populations, connection probabilities inversely scaled based on distance , 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.
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.
2.5. 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 (Naf); 5. E neuron Kdr channel density; 6. E neuron Ka channel density; 7. E neuron KD channel density; 8. E neuron KM 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).
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).
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 Naf 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.
Table 8. Dynamic measures (average ± standard deviation) for All simulations (n = 5867), Active simulations (n = 4341), Latch-up simulations (n = 1077), Active/Non-Latch-up (n = 3264), physiological simulations (n = 65), and dystonia simulations (n = 65).
2.6. 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 τ non-parametric 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 (nc) and discordant pairs (nd); ties are taken into account by the normalizing term, , 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 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).
3.1. 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.
3.2. 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). 1077 simulations were classified as epileptiform based on activity latch-up resulting in sustained periods. The different classes of simulations formed distinct clusters in multiple slices of excitatory corticospinal (ESP) activity feature-space (Figure 2). Physiological simulations showed E5P rates ≤ 2 Hz with low to intermediate E5P firing vector (FV) similarity. Dystonia simulations primarily occupied the upper-right quadrant of the scatterplot in 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.
Figure 2. Distinct dynamics in in physiological, dystonia, and latch-up simulations. (A) Average E5P firing rate vector (FV) similarity vs. average E5P firing rate for individual simulations. (B) E5P MUA Beta oscillation amplitude vs. E5P synchrony for individual simulations. (A,B) [light blue: physiological, purple: dystonia, orange: epileptiform, black: remaining Active simulations, large circles represent simulations shown in (C) and Figure 3]. (C) Pearson correlations between all pairs of E5P FVs. Solid black lines demarcate FVs from example physiological, dystonia, and epileptiform simulations. All FVs used 50 ms intervals, forming 40 FVs per 2 s of simulation.
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).
Figure 3. Distinct firing patterns in physiological, dystonia, and epileptiform (epileptic) simulations. (A) Physiological model has sparse, asynchronous E5P firing, relatively low I5 firing, and activated E5a/E5b populations. (B) Pathological model shows high-frequency, synchronous activity in E5P neurons, causing higher I5 firing, which suppresses E5a activity. (C) “Epileptiform” (epileptic) model shows high-frequency, synchronous activity with intermittent 200–300 ms gaps in firing of E neurons, caused by depolarization blockade (Na+ channel inactivation). (A–C) Left Dots represent individual neuron spike times (red: E cells, blue: LTS cells, green: FS cells). Cells arranged from layer 2/3 (top) to layer 6 (bottom). Scale-bar: 100 ms. (A–C) Right Population firing rates (25 ms bins) arranged vertically to roughly correspond to position on raster plot to the left. Scale-bar: 40 Hz (Same color code as raster; apparently flat lines indicate low variation in firing rate). (D) Population firing rates from simulations in (A–C) (Average ± standard error of the mean).
Figure 4. Motor cortex models produce different beta oscillations. Power spectrum of multiunit activity vectors of examples in Figure 3. Power (y-axis) in arbitrary units. (A) Physiological model shows weak beta (22 Hz) oscillations with power of < 0.1% of the pathological model. (B) Pathological model produces strong beta (20 Hz) oscillations with additional harmonic at 40 Hz. (C) Epileptiform model produces strong beta (19 Hz) oscillations with additional harmonic at 38 Hz.
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.
3.3. Need for Multitarget Approach
No individual parameter determined physiological vs. dystonia-dynamical-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.
Figure 5. Individual parameters do not distinguish physiological from dystonia activity. (A) Dystonia (purple) vs. physiological (light blue) simulations. of simulations sorted by E5P firing rate (N = 65 for each group). (B) Cumulative probability distributions for each parameter in the dystonia (purple) and physiological (light blue) simulations. Parameter values normalized to a distribution with zero mean and unit variance (zero mean does not indicate zero density of a given ion channel/receptor). Simulations shown are obtained from bottom and top 2nd percentile based on dynamic measures.
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 Ca2+ 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 (upper-right 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.
Figure 6. All parameters of pathological and physiological simulations reveals weak intra-class clustering. (A) 11-dimensional parameters for physiological and pathological simulations. Colorbar is normalized parameter values as in Figure 5. (B) Pearson correlations between all pairs of parameter vectors.
3.4. 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 Naf density—the strongest predictor of excitability. The best prediction came with high Naf and low Kdr. Logistic regression methods were also tried to perform this two-dimensional separation but did not perform as well as SVM.
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 Naf 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).
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 Naf, 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).
Figure 8. SVM classification accuracy generally increases when using 1–10 parameters, indicating utility of multitarget pharmacy approach to treating dystonia. (A) Best classification accuracy from all combinations of x parameters (solid line: mean cross-validation accuracy (n = 10); dotted line: standard error). (B) Best parameter (param) combinations (red: parameter used; blue: parameter not used). x-axis in (A,B) indicates number of parameters used.
This multi-target SVM approach revealed the parameters that had the highest contribution to producing or preventing dystonia. Naf 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 Kdr. Following these came Ka, Kd, BK, SK, and VGCC densities which also significantly contributed to accurate predictions, due to their strong influence on E neuron excitability. mGLUR, RYR, and Km densities showed lesser contributions.
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.
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.
We developed a multiscale model of primary motor cortex to explore multitarget pharmacological therapies for dystonia. We searched parameter space—channel and receptor densities—to 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, Naf and Kdr 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.
4.1. 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 averaging—averaging 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 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 medicine—identifying 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.
4.2. 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.
All authors listed, have made substantial, direct and intellectual contribution to the work, and approved it for publication.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
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.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fphar.2016.00157
Action, A. C. S. (2016). In silico Clinical Trials: How Computer Simulation will Transform the Biomedical Industry. Available online at: http://avicenna-isct.org/wp-content/uploads/2016/01/AvicennaRoadmapPDF-27-01-16.pdf (Accessed May 2, 2016).
Air, E., Ostrem, J., Sanger, T., and Starr, P. (2011). Deep brain stimulation in children: experience and technical pearls: clinical article. J. Neurosurg. Pediatr. 8, 566–574. doi: 10.3171/2011.8.PEDS11153
Anderson, C. T., Sheets, P. L., Kiritani, T., and Shepherd, G. M. G. (2010). Sublayer-specific microcircuits of corticospinal and corticostriatal neurons in motor cortex. Nat. Neurosci. 13, 739–744. doi: 10.1038/nn.2538
Apicella, A. J., Wickersham, I. R., Seung, H. S., and Shepherd, G. M. (2012). Laminarly orthogonal excitation of fast-spiking and low-threshold-spiking interneurons in mouse motor cortex. J. Neurosci. 32, 7021–7033. doi: 10.1523/JNEUROSCI.0011-12.2012
Ashhad, S., and Narayanan, R. (2013). Quantitative interactions between the A-type K+ current and inositol trisphosphate receptors regulate intraneuronal Ca2+ waves and synaptic plasticity. J. Physiol. 591, 1645–1669. doi: 10.1113/jphysiol.2012.245688
Bhanpuri, N. H., Bertucco, M., Ferman, D., Young, S. J., Liker, M. A., Krieger, M. D., et al. (2014). Deep brain stimulation evoked potentials may relate to clinical benefit in childhood dystonia. Brain Stimul. 7, 718–726. doi: 10.1016/j.brs.2014.06.003
Casellato, C., Maggioni, S., Lunardini, F., Bertucco, M., Pedrocchi, A., and Sanger, T. (2014). “Dystonia: altered sensorimotor control and vibro-tactile emg-based biofeedback effects,” in XIII Mediterranean Conference on Medical and Biological Engineering and Computing 2013 (Seville: Springer), 1742–1746.
Chadderdon, G. L., Mohan, A., Suter, B. A., Neymotin, S. A., Kerr, C. C., Francis, J. T., et al. (2014). Motor cortex microcircuit simulation based on brain activity mapping. Neural Comput. 26, 1239–1262. doi: 10.1162/NECO_a_00602
Chen, S., Wang, J., and Siegelbaum, S. (2001). Properties of hyperpolarization-activated pacemaker current defined by coassembly of HCN1 and HCN2 subunits and basal modulation by cyclic nucleotide. J. Gen. Physiol. 117, 491–504. doi: 10.1085/jgp.117.5.491
Crone, N. E., Miglioretti, D. L., Gordon, B., Sieracki, J. M., Wilson, M. T., Uematsu, S., et al. (1998). Functional mapping of human sensorimotor cortex with electrocorticographic spectral analysis. i. alpha and beta event-related desynchronization. Brain 121, 2271–2299. doi: 10.1093/brain/121.12.2271
Crowell, A. L., Ryapolova-Webb, E., Ostrem, J. L., Galifianakis, N. B., Shimamoto, S., Lim, D. A., et al. (2012). Oscillations in sensorimotor cortex in movement disorders: an electrocorticography study. Brain 135, 615–630. doi: 10.1093/brain/awr332
De Young, G. W., and Keizer, J. (1992). A single-pool inositol 1, 4, 5-trisphosphate-receptor-based model for agonist-stimulated oscillations in Ca2+ concentration. Proc. Natl. Acad. Sci. U.S.A. 89, 9895–9899.
Dura-Bernal, S., Li, K., Neymotin, S., Francis, J., Principe, J., and Lytton, W. (2016). Restoring behavior via inverse neurocontroller in a lesioned cortical spiking model driving a virtual arm. Front. Neurosci. 10:28. doi: 10.3389/fnins.2016.00028
Fitzpatrick, J. S., Hagenston, A. M., Hertle, D. N., Gipson, K. E., Bertetto-D'Angelo, L., and Yeckel, M. F. (2009). Inositol-1,4,5-trisphosphate receptor-mediated Ca2+ waves in pyramidal neuron dendrites propagate through hot spots and cold spots. J. Physiol. 587, 1439–1459. doi: 10.1113/jphysiol.2009.168930
Harnett, M. T., Magee, J. C., and Williams, S. R. (2015). Distribution and function of HCN channels in the apical dendritic tuft of neocortical pyramidal neurons. J. Neurosci. 35, 1024–1037. doi: 10.1523/JNEUROSCI.2813-14.2015
Harnett, M., Xu, N., Magee, J. C., and Williams, S. R. (2013). Potassium channels control the interaction between active dendritic integration compartments in layer 5 cortical pyramidal neurons. Neuron 79, 516–529. doi: 10.1016/j.neuron.2013.06.005
Hay, E., Hill, S., Schürmann, F., Markram, H., and Segev, I. (2011). Models of neocortical layer 5b pyramidal cells capturing a wide range of dendritic and perisomatic active properties. PLoS Comput. Biol. 7:e1002107. doi: 10.1371/journal.pcbi.1002107
Hooks, B. M., Mao, T., Gutnisky, D. A., Yamawaki, N., Svoboda, K., and Shepherd, G. M. (2013). Organization of cortical and thalamic input to pyramidal neurons in mouse motor cortex. J. Neurosci. 33, 748–760. doi: 10.1523/JNEUROSCI.4338-12.2013
Johnson, M. D., Miocinovic, S., McIntyre, C. C., and Vitek, J. L. (2008). Mechanisms and targets of deep brain stimulation in movement disorders. Neurotherapeutics 5, 294–308. doi: 10.1016/j.nurt.2008.01.010
Kerr, C. C., Neymotin, S. A., Chadderdon, G. L., Fietkiewicz, C. T., Francis, J. T., and Lytton, W. (2012). Electrostimulation as a prosthesis for repair of information flow in a computer model of neocortex. IEEE Trans. Neural Syst. Rehabil. Eng. 20, 153–160. doi: 10.1109/TNSRE.2011.2178614
Kerr, C. C., van Albada, S. J., Neymotin, S. A., Chadderdon, G. L., Robinson, P. A., and Lytton, W. W. (2013). Cortical information flow in Parkinson's disease: a composite network/field model. Front. Comput. Neurosci. 7:39. doi: 10.3389/fncom.2013.00039
Kiritani, T., Wickersham, I. R., Seung, H. S., and Shepherd, G. M. (2012). Hierarchical connectivity and connection-specific dynamics in the corticospinal-corticostriatal microcircuit in mouse motor cortex. J. Neurosci. 32, 4992–5001. doi: 10.1523/JNEUROSCI.4759-11.2012
Kole, M. H., Hallermann, S., and Stuart, G. J. (2006). Single Ih channels in pyramidal neuron dendrites: properties, distribution, and impact on action potential output. J. Neurosci. 26, 1677–1687. doi: 10.1523/JNEUROSCI.3664-05.2006
Kristeva, R., Chakarov, V., Losch, F., Hummel, S., Popa, T., and Schulte-Mönting, J. (2005). Electroencephalographic spectral power in writer's cramp patients: evidence for motor cortex malfunctioning during the cramp. Neuroimage 27, 706–714. doi: 10.1016/j.neuroimage.2005.05.004
Lefort, S., Tomm, C., Floyd-Sarria, J. C., and Petersen, C. C. (2009). The excitatory neuronal network of the C2 barrel column in mouse primary somatosensory cortex. Neuron 61, 301–316. doi: 10.1016/j.neuron.2008.12.020
Mallet, N., Pogosyan, A., Sharott, A., Csicsvari, J., Bolam, J. P., Brown, P., et al. (2008). Disrupted dopamine transmission and the emergence of exaggerated beta oscillations in subthalamic nucleus and cerebral cortex. J. Neurosci. 28, 4795–4806. doi: 10.1523/JNEUROSCI.0123-08.2008
Markram, H., Lübke, J., Frotscher, M., Roth, A., and Sakmann, B. (1997). Physiology and anatomy of synaptic connections between thick tufted pyramidal neurones in the developing rat neocortex. J. Physiol. 500, 409–440.
Meisel, C., Schulze-Bonhage, A., Freestone, D., Cook, M. J., Achermann, P., and Plenz, D. (2015). Intrinsic excitability measures track antiepileptic drug action and uncover increasing/decreasing excitability over the wake/sleep cycle. Proc. Natl. Acad. Sci. U.S.A. 112, 14694–14699. doi: 10.1073/pnas.1513716112
Migliore, M., Messineo, L., and Ferrante, M. (2004). Dendritic Ih selectively blocks temporal summation of unsynchronized distal inputs in CA1 pyramidal neurons. J. Comput. Neurosci. 16, 5–13. doi: 10.1023/B:JCNS.0000004837.81595.b0
Miller, K. J., Leuthardt, E. C., Schalk, G., Rao, R. P., Anderson, N. R., Moran, D. W., et al. (2007). Spectral changes in cortical surface potentials during motor movement. J. Neurosci. 27, 2424–2432. doi: 10.1523/JNEUROSCI.3886-06.2007
Monyer, H., and Markram, H. (2004). Interneuron diversity series: molecular and genetic tools to study gabaergic interneuron diversity and function. Trends Neurosci. 27, 90–97. doi: 10.1016/j.tins.2003.12.008
Neychev, V. K., Fan, X., Mitev, V. I., Hess, E. J., and Jinnah, H. A. (2008). The basal ganglia and cerebellum interact in the expression of dystonic movement. Brain 131, 2499–2509. doi: 10.1093/brain/awn168
Neymotin, S. A., Hilscher, M. M, Moulin, T., Skolnick, Y., Lazarewicz, M. T., and Lytton, W. W. (2013). Ih tunes theta/gamma oscillations and cross-frequency coupling in an in silico CA3 model. PLoS ONE 8:e76285. doi: 10.1371/journal.pone.0076285
Neymotin, S., Lazarewicz, M., Sherif, M., Contreras, D., Finkel, L., and Lytton, W. (2011b). Ketamine disrupts theta modulation of gamma in a computer model of hippocampus. J. Neurosci. 31, 11733–11743. doi: 10.1523/JNEUROSCI.0501-11.2011
Neymotin, S., Lee, H., Park, E., Fenton, A., and Lytton, W. (2011c). Emergence of physiological oscillation frequencies in a computer model of neocortex. Front. Comput. Neurosci. 5:19. doi: 10.3389/fncom.2011.00019
Neymotin, S. A., McDougal, R. A., Bulanova, A. S., Zeki, M., Lakatos, P., Terman, D., et al. (2016). Calcium regulation of HCN channels supports persistent activity in a multiscale model of neocortex. Neuroscience 316, 344–366. doi: 10.1016/j.neuroscience.2015.12.043
Neymotin, S. A, McDougal, R. A., Sherif, M. A., Fall, C. P., Hines, M. L., and Lytton, W. W. (2015). Neuronal calcium wave propagation varies with changes in endoplasmic reticulum parameters: a computer model. Neural Comput. 27, 898–924. doi: 10.1162/NECO_a_00712
Nobrega, J. N., Richter, A., Burnham, W. M., and Lôscher, W. (1995). Alterations in the brain GABAA/benzodiazepine receptor-chloride ionophore complex in a genetic model of paroxysmal dystonia: a quantitative autoradiographic analysis. Neuroscience 64, 229–239. doi: 10.1016/0306-4522(94)00334-2
Orrù, G., Pettersson-Yeo, W., Marquand, A. F., Sartori, G., and Mechelli, A. (2012). Using support vector machine to identify imaging biomarkers of neurological and psychiatric disease: a critical review. Neurosci. Biobehav. Rev. 36, 1140–1152. doi: 10.1016/j.neubiorev.2012.01.004
Pfurtscheller, G., and Aranibar, A. (1979). Evaluation of event-related desynchronization (ERD) preceding and following voluntary self-paced movement. Electroencephalogr. Clin. Neurophysiol. 46, 138–146.
Pratt, D., Richter, A., Möhler, H., and Löscher, W. (1995). Regionally selective and age-dependent alterations in benzodiazepine receptor binding in the genetically dystonic hamster. J. Neurochem. 64, 2153–2158.
Richter, A., and Hamann, M. (2001). Effects of adenosine receptor agonists and antagonists in a genetic animal model of primary paroxysmal dystonia. Br. J. Pharmacol. 134, 343–352. doi: 10.1038/sj.bjp.0704268
Robinson, R. B., and Siegelbaum, S. A. (2003). Hyperpolarization-activated cation currents: from molecules to physiological function. Annu. Rev. Physiol. 65, 453–480. doi: 10.1146/annurev.physiol.65.092101.142734
Ross, W. N., Nakamura, T., Watanabe, S., Larkum, M., and Lasser-Ross, N. (2005). Synaptically activated Ca2+ release from internal stores in CNS neurons. Cell Mol. Neurobiol. 25, 283–295. doi: 10.1007/s10571-005-3060-0
Safiulina, V. F., Caiati, M. D., Sivakumaran, S., Bisson, G., Migliore, M., and Cherubini, E. (2010). Control of GABA release at mossy fiber-CA3 connections in the developing hippocampus. Front. Synaptic Neurosci. 2:1. doi: 10.3389/neuro.19.001.2010
Sanger, T. (2003). Childhood onset generalised dystonia can be modelled by increased gain in the indirect basal ganglia pathway. J. Neurol. Neurosurg. Psychiatry 74, 1509–1515. doi: 10.1136/jnnp.74.11.1509
Sanger, T., Delgado, M., Gaebler-Spira, D., Hallett, M., Mink, J., and on Childhood Motor Disorders, T. F. (2003). Classification and definition of disorders causing hypertonia in childhood. Pediatrics 111, e89–e97. doi: 10.1542/peds.111.1.e89
Sanger, T. D., Kukke, S. N., and Sherman-Levine, S. (2007). Botulinum toxin type B improves the speed of reaching in children with cerebral palsy and arm dystonia: an open-label, dose-escalation pilot study. J. Child Neurol. 22, 116–122. doi: 10.1177/0883073807299975
Schwindt, P. C., Spain, W. J., and Crill, W. E. (1992). Effects of intracellular calcium chelation on voltage-dependent and calcium-dependent currents in cat neocortical neurons. Neuroscience 47, 571–578. doi: 10.1016/0306-4522(92)90166-Y
Sneyd, J., Tsaneva-Atanasova, K., Bruce, J. I., Straub, S. V., Giovannucci, D. R., and Yule, D. I. (2003). A model of calcium waves in pancreatic and parotid acinar cells. Biophys. J. 85, 1392–1405. doi: 10.1016/S0006-3495(03)74572-X
Song, W., Kerr, C. C., Lytton, W. W., and Francis, J. T. (2013). Cortical plasticity induced by spike-triggered microstimulation in primate somatosensory cortex. PLoS ONE 8:e57453. doi: 10.1371/journal.pone.0057453
Spruston, N., Schiller, Y., Stuart, G., and Sakmann, B. (1995). Activity-dependent action potential invasion and calcium influx into hippocampal CA1 dendrite. Science 8, 297–300. doi: 10.1126/science.7716524
Stacey, W., Lazarewicz, M., and Litt, B. (2009). Synaptic noise and physiological coupling generate high-frequency oscillations in a hippocampal computational model. J. Neurophysiol. 102, 2342–2357. doi: 10.1152/jn.00397.2009
Toro, C., Deuschl, G., Thatcher, R., Sato, S., Kufta, C., and Hallett, M. (1994). Event-related desynchronization and movement-related cortical potentials on the ECoG and EEG. Electroencephalogr. Clin. Neurophysiol. 93, 380–389.
Viceconti, M., Clapworthy, G., and Jan, S. (2008). The Virtual Physiological Human - a European initiative for in silico human modelling. J. Physiol. Sci. 58, 441–446. doi: 10.2170/physiolsci.RP009908
Viceconti, M, Edwin Morley-Fletcher, M., and Henney, A. (2016). In silico Clinical Trials: How Computer Simulation will Transform the Biomedical Industry. Available online at: http://avicenna-isct.org/wp-content/uploads/2016/01/AvicennaRoadmapPDF-27-01-16.pdf (Accessed May 2, 2016).
Wang, J., Chen, S., Nolan, M. F., and Siegelbaum, S. A. (2002). Activity-dependent regulation of HCN pacemaker channels by cyclic AMP: signaling through dynamic allosteric coupling. Neuron 36, 451–461. doi: 10.1016/S0896-6273(02)00968-6
Winograd, M., Destexhe, A., and Sanchez-Vives, M. (2008). Hyperpolarization-activated graded persistent activity in the prefrontal cortex. Proc. Natl. Acad. Sci. U.S.A. 105, 7298–7303. doi: 10.1073/pnas.0800360105
Keywords: dystonia, multiscale modeling, computer simulation, motor cortex, beta oscillations, corticospinal neurons, multitarget pharmacology, support vector machines
Citation: Neymotin SA, Dura-Bernal S, Lakatos P, Sanger TD and Lytton WW (2016) Multitarget Multiscale Simulation for Pharmacological Treatment of Dystonia in Motor Cortex. Front. Pharmacol. 7:157. doi: 10.3389/fphar.2016.00157
Received: 17 February 2016; Accepted: 30 May 2016;
Published: 14 June 2016.
Edited by:Hugo Geerts, In Silico Biosciences, USA
Reviewed by:Gaute T. Einevoll, Norwegian University of Life Sciences, Norway
Christian Meisel, University Hospital Carl Gustav Carus Dresden, Germany
James B. Aimone, Sandia National Laboratories, USA
Copyright © 2016 Neymotin, Dura-Bernal, Lakatos, Sanger and Lytton. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Samuel A. Neymotin, email@example.com