Deciphering in silico the Role of Mutated NaV1.1 Sodium Channels in Enhancing Trigeminal Nociception in Familial Hemiplegic Migraine Type 3

Familial hemiplegic migraine type 3 (FHM3) is caused by gain-of-function mutations in the SCN1A gene that encodes the α1 subunit of voltage-gated NaV1.1 sodium channels. The high level of expression of NaV1.1 channels in peripheral trigeminal neurons may lead to abnormal nociceptive signaling thus contributing to migraine pain. NaV1.1 dysfunction is relevant also for other neurological disorders, foremost epilepsy and stroke that are comorbid with migraine. Here we used computer modeling to test the functional role of FHM3-mutated NaV1.1 channels in mechanisms of trigeminal pain. The activation of Aδ-fibers was studied for two algogens, ATP and 5-HT, operating through P2X3 and 5-HT3 receptors, respectively, at trigeminal nerve terminals. In WT Aδ-fibers of meningeal afferents, NaV1.1 channels efficiently participate in spike generation induced by ATP and 5-HT supported by NaV1.6 channels. Of the various FHM3 mutations tested, the L263V missense mutation, with a longer activation state and lower activation voltage, resulted in the most pronounced spiking activity. In contrast, mutations that result in a loss of NaV1.1 function largely reduced firing of trigeminal nerve fibers. The combined activation of P2X3 and 5-HT3 receptors and branching of nerve fibers resulted in very prolonged and high-frequency spiking activity in the mutants compared to WT. We identified, in silico, key determinants of long-lasting nociceptive activity in FHM3-mutated Aδ-fibers that naturally express P2X3 and 5-HT3 receptors and suggest mutant-specific correction options. Modeled trigeminal nerve firing was significantly higher for FHM3 mutations, compared to WT, suggesting that pronounced nociceptive signaling may contribute to migraine pain.


INTRODUCTION
The generation of disabling migraine pain involves the activation of the meningeal trigeminovascular system (Moskowitz, 2008;Messlinger, 2009), but the underlying pro-nociceptive mechanisms remain largely unknown. Current data suggest participation of mainly Aδ-fibers of the trigeminal nerve densely innervating the meninges (Melo-Carrillo et al., 2017;Haanes and Edvinsson, 2019).
To unravel molecular mechanisms of migraine pathophysiology, the study of monogenic subtypes of migraine, foremost familial hemiplegic migraine, has been instrumental . Familiar hemiplegic migraine type 3 (FHM3), which is caused by specific missense mutations in the SCN1A gene encoding the α1 subunit of voltage-gated Na V 1.1 sodium channels (Dichgans et al., 2005;Tolner et al., 2015), allows the specific interrogation of the role mutant Na V 1.1 channels in migraine pathophysiology. The observation that incubation at lower temperature and expression in neurons rescued folding/trafficking issues now firmly established that FHM3 is caused by a gain of Na V 1.1 function (Dhifallah et al., 2018), whereas earlier studies suggested foremost loss-offunction effects of FHM3 mutations when overexpressed in heterologous expression systems (Dichgans et al., 2005;Kahlig et al., 2008). As Na V 1.1 channels are strongly expressed in peripheral Aδ-fibers of the trigeminal nerve (Ho and O'Leary, 2011;Osteen et al., 2016), their modified activity may underlie the activation of peripheral trigeminal neurons leading to migraine pain. It can be expected that a gain-of-function enhances excitability of peripheral nerve fibers expressing modified Na V 1.1 channels, providing a high pro-nociceptive activity delivered to second order neurons.
Investigating Na V 1.1 channel dysfunction, which are expressed also in central neurons (Ogiwara et al., 2007;Favero et al., 2018;Sakkaki et al., 2020), in relation to changes in neuronal excitability, is also of relevance to various neurological disorders other than migraine, foremost childhood epilepsy (Menezes et al., 2020), autism spectrum disorder (Scheffer and Nabbout, 2019), Alzheimer's disease (Sakkaki et al., 2020), and perhaps less known, transient cerebral ischemia (Zhan et al., 2007). In fact, changes in neuronal hyperexcitability may therefore, at least to certain extent, underlie the comorbidity of several of the disorders with migraine, for instance when they lead to spreading depolarizations, as observed for migraine and stroke (Dreier and Reiffurth, 2015). Unlike for the other disorders, the evidence for a specific role of Na V 1.1 channels in cerebral ischemia is limited although voltage-gated cation channels, including sodium channels have been targets for the treatment of stroke (Gribkoff and Winquist, 2005).
Recently, we presented a mathematical model of the nociceptive neuro-immune synapse in meninges that, by activation with ATP and 5-HT, generates neuronal firing (Suleimanova et al., 2020). Meninges, which are densely innervated by trigeminal nerve fibers, are currently considered a main source of migraine headache (Moskowitz, 2008;Messlinger, 2009;Olesen et al., 2009). The two algesic substances, ATP and 5-HT selected to be modeled in this study produce a powerful and long-lasting activation of meningeal afferents (Yegutkin et al., 2016;Kilinc et al., 2017;Koroleva et al., 2019).
The data are consistent with the purinergic hypothesis of migraine suggested earlier by Burnstock (1981). There is strong evidence that serotonin (5-HT) is involved in migraine already because "triptans, " serotonin (5-HT1B/1D) agonists, are effective in treating migraine patients (Goadsby, 2007). However, the sources of the two neurotransmitters, the time spent in the extracellular space, the receptor kinetics, and most notably, the speed of desensitization of the transmitters are different, as we presented in our previous model (Suleimanova et al., 2020). Notably, the respective P2X3 and 5-HT3 receptors activated by ATP and 5-HT, respectively, are expressed in Aδ-fibers (Ford, 2012;Kilinc et al., 2017;Sato et al., 2018), thus in the same neurons that express Na V 1.1 sodium channels. Combining, in the mathematical model, effects of these triggers with dysfunction of Na V 1.1 channels due to mutations identified in patients may serve as a useful platform to explore mechanisms of activation of peripheral nociception relevant to migraine headache. Therefore, we here assessed, by using in silico modeling, how gain-and loss-of-function mutations in the α1 subunit of Na V 1.1 channels, as they occur in patients with FHM3 and childhood epilepsy, respectively, might affect peripheral trigeminal nociception in meningeal afferents. To this end, the WT vs. various mutants Na V 1.1 channels were in silico "coexpressed" along with other types of sodium (Na V 1.6, Na V 1.7, Na V 1.8) and several potassium channel types (K V 1, K V 3, K V 4, and calcium-activated potassium channel K Ca ) according to published profile of these channels in Aδ-fibers (Tigerholm et al., 2014;Mandge and Manchanda, 2018;Zemel et al., 2018;Zheng et al., 2019). Our data show a large amplification of nociception in meningeal afferents exclusively with gain-of-function mutations providing a scientific framework that a peripheral mechanism of generating migraine pain in patients with FHM3.

The Mathematical Model for Testing of FHM3 Mutations
To model the function of the trigeminal nerve in meninges we used the NEURON environment version 7.8 (Hines and Carnevale, 2003). Our model describes the activity of Aδ-fibers induced by single, or repetitive ATP and 5-HT release events from abundant meningeal mast cells (Theoharides et al., 1995;Levy, 2009;Kilinc et al., 2017). Thus, the Aδ-fiber coupled to a mast cell represent a model of the nociceptive "neuroimmune synapse" (Theoharides et al., 1995;Koroleva et al., 2019;Suleimanova et al., 2020; Figure 1A). Within such synapse, locally released ATP or 5-HT activates, at nerve terminals, P2X3 or 5-HT3 receptors, respectively. Both ATP and 5-HT are potent triggers of nociceptive firing in meningeal afferents (Yegutkin et al., 2016;Kilinc et al., 2017;Koroleva et al., 2019). Although P2X2 receptors are co-expressed with P2X3 subunits in rodents (Simonetti et al., 2006), the P2X3 subtype is the predominant ATP receptor subtype in human sensory neurons (Serrano et al., 2012). Therefore, in our model, we used P2X3 receptors as the main target for fast action of ATP on trigeminal meningeal afferents. We explored the action of different concentrations of two algogens but for most model trials we used 1 µM ATP and 2 µM 5-HT concentrations to activate receptors, since these values are close to EC 50 of the respective receptors (Sokolova et al., 2006;Corradi et al., 2009). The kinetics of P2X3 receptors was based on the previously published model of this receptor (Sokolova et al., 2006), whereas the kinetic of 5-HT3 receptors was taken from the study of Corradi et al. (2009). The lifetime of extracellular ATP is determined by fast hydrolysis via multiple ecto-enzymes (Yegutkin, 2008), whereas the profile of 5-HT is controlled by a relatively slow uptake via SERT transporters (Wood et al., 2014). Therefore, we included in our model the partial hydrolysis for ATP and uptake for 5-HT making the model consistent with experimental data on the action of ATP and 5-HT (Suleimanova et al., 2020). The 3D diffusion model suggested by Saftenku (2005) was used to reproduce the time course of transmitters in the nociceptive synapse. We also conducted in the current model experiments in which we varied ATP and 5-HT concentrations to explore the dependence of spike firing on the profile of the neurotransmitters in the meningeal neuroimmune synapse.
The diameter of the Aδ-fiber is 5 µm (West et al., 2015), each segment of the fiber consists of a paranode and a node of Ranvier (McIntyre et al., 2002) that is 2 × 2 µm (width and length; Figure 1A). Aδ-fibers highly expressed tetrodotoxin-sensitive (TTX-sensitive) Na V 1 channel type (Na V 1.1, Na V 1.6, Na V 1.7), while Na V 1.8 was at low level compared to C-fibers (Pinto et al., 2008;Zhang et al., 2013). Most of the input parameters on channel kinetics and voltage dependence of Na V 1.6, K V 1, K V 3, and K V 4 ion channels were obtained from somatic recordings (Zheng et al., 2019), whereas the data for the density of ion channels were obtained from axonal measurements (Waxman and Ritchie, 1993;McIntyre et al., 2002). Likewise, we used data from somatic recordings of the contribution of Na V 1 channels to I Na current (Zhang et al., 2013) and Na V channels mRNA levels (Ho and O'Leary, 2011). The WT Na V 1.1 channel and FHM3-associated mutations in this channel were modeled based on functional properties obtained in neurons or tsA201 cells (for details see Table 1).
To match the natural profile of sodium channels in the Aδ-fiber, we revised our previously used mathematical model (Suleimanova et al., 2020) by adding Na V 1.1 and Na V 1.6, since V 1/2 is the voltage of half-maximal activation or inactivation; k is a slope factor; tau slow is kinetics of the development of slow inactivation. Voltage indicated in brackets shows values at which the recordings were performed in the quoted papers and in model studies.
these fibers express both type of these channels in nodes of Ranvier (Duflocq et al., 2008;Leterrier et al., 2011). We also added Na V 1.7 and Na V 1.8 channels to our model, as they also support the generation and propagation of nociceptive signals from the periphery (Duflocq et al., 2008;Black et al., 2012; Figure 1B and Supplementary Table 1). The maximum conductance densities of the TTX-sensitive channels are based on single-channel conductance and channel density values in the range of 1,000 to 2,000 channels/µm 2 (McIntyre et al., 2002). The density of TTX-resistant Nav1.8 channel is significantly lower since TTX completely blocked spiking activity in Aδ-fiber (Pinto et al., 2008), whereas specific blockers of TTX-resistant channels did not significantly reduce spiking (Tsuchimochi et al., 2011). Of note, we tuned the values to provide a closer similarity to action potential and firing patterns observed in the experimental studies (Cestele et al., 2008(Cestele et al., , 2013 taking into account the contributions of Na V 1.1, Na V 1.6, and Na V 1.7 to TTX-sensitive current (Zhang et al., 2013) and therespective Na V channels mRNA level (Ho and O'Leary, 2011). When modeling, we also took into consideration a fiber-specific difference in expression of potassium channels (Zemel et al., 2018;Chien et al., 2007). Thus, our Aδ-fiber model includes voltage-gated potassium channels K V 1, Kv3, and Kv4 mediating A-type currents, which are found in axons and nerve endings (Zemel et al., 2018). We recently showed that 4-Aminopyridine, a blocker of A-type current (determined by voltage-gated potassium Kv1 and Kv3 channels) dramatically enhanced the firing of trigeminal afferents in rat meninges promoting appearance of fast large spikes (Andreou et al., 2020), typical for the phenotype of trigeminal Aδ-fibers (MacIver and Tanelian, 1993). We also added to the model calciumactivated potassium large-conductance K Ca channels (BK Ca ) (Mandge and Manchanda, 2018). The full set of maximum conductance densities of the channels and the kinetics underlying the single action potential is shown in Supplementary Table 1 and Figure 1B, where the current time course for each mutant (L263V, L1670W, and M145T) was compared with the respective WT to fit with experimental results from Mantegazza et al. (2005); Kahlig et al. (2006), and Dhifallah et al. (2018). Next, we modeled various SCN1A mutations by using steadystate voltage dependence and kinetics of channel activation and inactivation as presented in Spampanato et al. (2004). First, we fitted the kinetics of the channels using the threeparameter Gaussian exponential function: τ = a × exp (((V-V 0 )/b) 2 ), where V is the membrane voltage, V 0 is the mean voltage around which the Gaussian is positioned, a is the height, and b is the standard deviation. Second, the steady-state voltage dependence of activation was calculated by a Boltzmann function in the form: m ∞ = 1/(1 + exp(-(V-V 1/2 )/k) and h ∞ /s ∞ /r ∞ = 1/(1 + exp((V-V 1/2 )/k) for fast (h ∞ ), slow (s ∞ ) inactivation and recovery (r ∞ ), where V 1/2 is the voltage of halfmaximal activation or inactivation, and k is a slope factor. The constants for these calculations are presented in Table 1. Third, we used the Hodgkin-Huxley formalism where activation (m) and inactivation (h) states were in the range between 0 and 1 dependent on voltage and time. We used the following equation to define the channel activation state: dm/dt = (m ∞ -m)/τ m , where m ∞ is referred to as steady-state voltage dependence and τ m is the kinetics. In a similar way, we determined the fast and slow inactivation states and recovery. We assume that the slow inactivation and recovery states are independent from other channel gating states since development of slow inactivation and recovery was recorded separately in experimental papers, we added s-and r-variable to the model. Thus, the sodium current is described by following equation: where g Na is the maximum conductance; E Na is the equilibrium potential for sodium; m and h are the voltage-dependent activation and inactivation gating states; s and r are the additional slow inactivation development and recovery voltage-dependent variables.
We modified the activation voltage and the activation speed of WT Na V 1.1 channels (Zheng et al., 2019) to recapitulate better the characteristics of a FHM3 mutation. To this end, we decreased the activation voltage of FHM3associated Na V 1.1 channels (V 1/2 of WT = −21 mV, V 1/2 of FHM3 = −25 mV; Figures 2Aa,Ab and Table 1) since L263V-mutated Na V 1.1 (Na V 1.1-L263V) channels are activated at a lower voltage than WT Na V 1.1 channels, although the difference was not significant (Kahlig et al., 2008). The FHM3 L263V mutation shows slower kinetics than WT for both the activation and the inactivation state and an increased "window current" (Figure 2Aa and Supplementary Figure 2A) due to the shift of the voltage dependence of activation curve to lower voltages and the inactivation curve to higher voltages compared to WT (Kahlig et al., 2008), thus the current peak of Na V 1.1-L263V is wide (Figure 1Ca and Supplementary Figure 3). Second, we modeled the FHM3 L1670W mutation (Figure 1Cb), which also exhibits a gainof-function effect (Dhifallah et al., 2018). This mutation increases the persistent sodium current and shows a positive shift of the inactivation curve and faster recovery from fast inactivation in comparison to WT (Figures 2Ba,Bb and Supplementary Figures 1B, 2B). Third, the FHM3 L1649Q mutation also enhances persistent sodium current, recovery from slow inactivation, and slightly increased the  Kahlig et al. (2008) for the L263V mutation, Dhifallah et al. (2018) for the L1670W mutation, Cestele et al. (2008) for the L1649Q mutation, Cestele et al. (2013) for the Q1478K mutation, and Mantegazza et al. (2005) for the M145T mutation. Mutations were compared with the respective WT data from the publication that presented experimental data for the mutations.  Figures 1C, 2C), which can be the reason for the prolonged spiking activity (Cestele et al., 2013). Fourth, the FHM3 Q1478K mutation shows faster recovery, a positive shift of the voltage dependence of inactivation (Figures 2Da,Db and Supplementary Figures 1D, 2D), and a higher amplitude persistent current (Cestele et al., 2008). Next, we modeled to the model the R1648H mutation (Figures 2Ea,Eb and Supplementary Figures 1E, 2E) that is associated with childhood epilepsy but exhibits a persistent sodium current as well (Kahlig et al., 2006). Finally, for comparison, we modeled the loss-of-function mutation M145T (Figures 1Cc, 2Fa). The comparison of activation and inactivation properties of the WT and mutant Na V 1.1 channels are shown in Figure 2 and Supplementary Figures 1, 2. Using experimental data from Mantegazza et al. (2005); Kahlig et al. (2006Kahlig et al. ( , 2008; Cestele et al. (2008Cestele et al. ( , 2013, and Dhifallah et al. (2018) for the mutant and respective WT channel characteristics, we validated our model data with experimental results describing the voltage-dependence of channel activation and inactivation by voltage steps (with increment 10 mV) from −100 to 20 mV. The time constant of development of fast inactivation are shown in Supplementary Figure 2. The distribution of tau in our model was compared with time constant from experimental results at the different potentials in the range (−65 to 30 mV). We used median values of normalized conductance at voltage from −100 to 20 mV with a 10-mV step increment. Based on these values, we plotted steady-state curves, using interpolation, and compared results with simulated voltagedependent steady-states. This comparison indicating a high similarity (p-values with Kolmogorov-Smirnov test (K-S test) higher than 0.1 (Boyerinas, 2016) between experimental results and simulation data (Figure 2, notice dotted lines for experimental results).
Finally, we took into consideration that the meningeal nerve has a branched structure (Schueler et al., 2014;Barkai et al., 2020;Suleimanova et al., 2020). Therefore, we modified our model to a tree structure with two 1.1 cm long axon branches. In this version of the model, we added two varying concentrations of ATP or 5-HT to activate each of the axon branches. As the junction of branches can block spikes from a branch in the refractory state, we calculated that an interval of 15 ms between applications was sufficient for the propagation of spiking activity.

Statistical Analysis
We used Kolmogorov-Smirnov test (K-S test) (Boyerinas, 2016) to compare the distribution of inter-spike intervals in experiments with simulation. The ks_2samp function from SciPy library that is implemented in the K-S test statistics was used to compare two samples. With this approach, the p-value higher than 0.05, allowed to accept the null hypothesis indicating that the distribution of data in two samples are similar. Experimental data for validation of the model were taken from our publication describing the action of ATP and 5-HT on mouse trigeminal nerves (Koroleva et al., 2019).

RESULTS
Role of Na V 1.1 and Na V 1.6 Channels in ATP-Induced Activation of WT and FHM3 Aδ-Fibers Aδ-fibers are thought to play an important role in the generation of migraine pain (Melo-Carrillo et al., 2017;Haanes and Edvinsson, 2019). As Na V 1.1 and Na V 1.6 channels are the major Na V channel types in Aδ-fibers, we first explored their role in ATP-induced activation of WT and FHM3 nerve fibers. We first used a simplified model of the single trigeminal nerve fiber with one release site for ATP. To address the pure Na V 1 phenotype we started with in silico "knocking down" the Na V 1.6 channel subtype in WT, whereas the conductance of Na V 1.1 was set to 0.35 S/cm 2 . In this case, ATP only induced local receptor potential (Figures 3Aa,Ba black lines). Adding Na V 1.6, with a conductance of 0.35 S/cm 2 was sufficient to generate propagating spikes (Figures 3Ab,Bb black lines). Raising the activity of Na V 1.1 (conductance 0.5 S/cm 2 ) enhanced the firing with repetitive spikes (Figures 3Ac,Bc black lines). In contrast to WT, in the FHM3 model with gainof-function mutations L263V and L1670W, even Na V 1.1 alone was sufficiently effective to generate propagating spikes (Figures 3Aa,Ba red lines). Even higher activity was obtained when Na V 1.6 (Figures 3Ab,Bb red lines) was added or when the conductance of the Na V 1.1 channel was increased to 0.5 S/cm 2 (Figures 3Ac,Bc red lines).
Thus, even in this simple linear model, the activity of FHM3 gain-of-function mutations produced an increased number of nociceptive spikes, suggesting that excitability of the terminals of meningeal afferents was increased. Notably, the L263V mutation caused more spikes than the L1670W mutation when compared with the respective WTs.
Specific Role of Na V 1.1 Channels in the Spiking Activity of WT and FHM3

Aδ-Fibers
As it is known that meningeal afferents are branched (Schueler et al., 2014;Suleimanova et al., 2020), we next modeled the trigeminal nerves with two branches and two ATP release events activating these separate nerve branches (Figure 4A). Such model allowed us to explore the role of simultaneous and shifted in time activation of distinct branches and whether this activity may propagate to higher pain centers. We found that in case of the L263V mutation, two simultaneous ATP release events did not change the outcome number of spikes coming to TG (Figure 4B red line). However, the final number of spikes was increased by adding a second ATP release with the time interval being 5 ms ( Figure 4C). With an interval of 10 ms, the number of repetitive spikes further increased from four in WT to eight spikes in L263V (Figure 4D). Similar FIGURE 3 | Testing the role of Na V 1.1 and Na V 1.6 in firing of Aδ-fiber with one branch and a single ATP release event. (Aa) The receptor potential in a single WT (Kahlig et al., 2008) fiber expressing Na V 1.1 (black line) and firing in the fiber without Na V 1.6. with L263V mutation (red line) of Na V 1.1 with conductance 0.35 S/cm 2 . Of note, there is a membrane potential in the nerve compartment, where ATP is acting, whereas in other panels we present the appearance of the signal to the trigeminal ganglion. (Ab) WT (Kahlig et al., 2008) fiber (black line) and the L263V mutation (red line) with Na V 1.1 and Na V 1.6 (both with conductance 0.35 S/cm 2 ) and; (Ac) WT (Kahlig et al., 2008) fiber (black line) and the L263V mutation (red line) with Na V 1.1 conductance 0.5 S/cm 2 and Na V 1.6 conductance 0.35 S/cm 2 . (Ba) The receptor potential in a single WT (Dhifallah et al., 2018) fiber (black line) and L1670W mutation (red line) of Na V 1.1 conductance 0.35 S/cm 2 and without Na V 1.6. (Bb) WT (Dhifallah et al., 2018) fiber (black line) and the L1670W mutation (red line) with Na V 1.1 conductance 0.35 S/cm 2 and Na V 1.6 conductance 0.35 S/cm 2 . (Bc) WT (Dhifallah et al., 2018) fiber (black line) and the L1670W mutation (red line) with Na V 1.1 conductance 0.5 S/cm 2 and Na V 1.6 conductance 0.35 S/cm 2 . ATP signaling is limited by partial hydrolysis.
increase was observed with interval 15 ms ( Figure 4E). Then, for modeling of other mutant channels, we used a 15-ms interval between ATP release events since a further increase of the interval in the range of 15-100 ms did not change the number of produced spikes. We compared the firing for FHM3 gain-of-function mutations L263V and L1670W, and loss-of-function mutation M145T with the respective WTs. The latter, even modeled from distinct experimental studies (Mantegazza et al., 2005;Kahlig et al., 2008;Dhifallah et al., 2018) generated similar (five for L263V and L1670W and four for M145T, see Figures 4E-G) repetitive spikes. We found that the mutation L263V was much more effective in promote firing (eight spikes, Figure 4E red line) than WT (five spikes, Figure 4E black line). Firing in case of the L1670W mutation was only slightly higher than for WT ( Figure 4F red line, WT-black line). Instead, the M145T mutation led to a dramatically reduced firing ( Figure 4G). Thus, FHM3 mutations further increased spike frequency activated by multiple ATP release sites on branched Aδ-fibers. The branching of the single axon enhanced the spiking activity, only when it was combined with asynchronous ATP release (Figure 4) because spikes from different branches were efficiently summarized in the primary afferents.
As the concentration of extracellular ATP can vary in the neuro-immune synapse due to expression profile and location of the powerful ATP-degrading ectoenzymes (Yegutkin, 2008), we next explored, in a branched model, the role of different concentrations of ATP in generation of a single spike and on repetitive firing in the most prominent mutant L263V. The dependence of a single-spike probability from concentration of ATP and the role of ATP hydrolysis is shown in Supplementary Figure 4. In the range of the tested concentrations, the mutant L263V had a higher probability of the spike generation, with a most visible difference at lower ATP concentrations (Supplementary Figure 4A, exemplified in Supplementary Figure 4C). In case of a lack of ATP hydrolysis, the difference at low ATP concentrations was more noticeable (Supplementary Figure 4B). The repetitive firing  (Mantegazza et al., 2005) fiber (black line) with an ATP release interval of 15 ms. Note that branching and multiple ATP release increased the firing in the case of a gain-of function but reduced in the case of a loss-of-function. ATP signaling is limited by partial hydrolysis.
was less sensitive to the absence or presence of ATP hydrolysis (Supplementary Figures 4G,H).

5-HT Induced Activation of WT and FHM3 Aδ-Fibers
5-HT is the major neurotransmitter released from meningeal mast cells that can powerfully activate peripheral trigeminal nerve terminals through ligand-gated 5-HT3 receptors (Kilinc et al., 2017;Koroleva et al., 2019). Therefore, we conducted modeling experiments (two branches, two release sites) simulating 5-HT release from mast cells and its action on 5-HT3 receptor in Aδ-fibers meningeal afferents.
Like with ATP, we found that the excitatory action of 5-HT was twice more pronounced with the L263V mutation than with the respective WT ( Figure 5A L263V-red line, WT-black line). The dependence of a single-spike probability from a concentration of 5-HT in basal conditions or after removal of 5-HT upate for the mutant L263V resembled, in some tests, that of ATP (Supplementary Figures 4D,E). In the example shown in Supplementary Figure 4F, a concentration of 5-HT as low as 0.6 µM was already able to generate a single spike in the mutant  (Kahlig et al., 2006) fiber (black line) associated with childhood epilepsy gain-of-function of Na V 1.1. (G) The M145T Na V 1.1 mutant (red line) and WT (Mantegazza et al., 2005) fiber (black line) associated with a loss-of-function mutation of Na V 1.1. ATP signaling is limited by partial hydrolysis and 5-HT signals are reduced by specific uptake.
L263V, whereas the same concentration of neurotransmitter in the WT induced only a small receptor potential. We found also an increased repetitive firing with higher 5-HT concentrations for mutant L263V; in particular, when we switched from 2 to 10 µM 5-HT, not only the number of spikes but also the firing frequency rose (Supplementary Figures 5E,F). Interestingly, for mutant L263V, increasing the concentration of 5-HT in the micromolar range had a stronger effect on repetitive firing than similar changes in the concentration of ATP (Supplementary  Figures 5B,F vs. Supplementary Figures 5A,D).
As ATP and 5-HT can be co-released during migraine events in the meninges (Suleimanova et al., 2020), we next modeled the simultaneous action of ATP and 5-HT. The combined activation of P2X3 and 5-HT3 receptors in WT produced more spiking activity than ATP or 5-HT alone (Figures 5B-G black lines). Then, we extended our approach to compare the firing of SCN1A gain-and loss-of-function mutations with their respective control WTs (Figures 5B-G red lines). With simultaneous release of ATP and 5-HT, the neuronal firing of WTs was similar (six spikes for L263V, R1648H, M145T, and seven spikes for Q1478K, and L1649Q, Figures 5B-G). Again, like in previous modeling conditions, the most pronounced spiking activity was observed with the FHM3 mutations L263V (Figure 5B), Q1478K (Figure 5E), and L1649Q ( Figure 5D). The R1648H mutation ( Figure 5F) that is associated with childhood epilepsy produced less spikes than the L263V mutation but more than WT. Thus, the updated Aδ-fiber model with simultaneous action of 5-HT and ATP on the branched trigeminal nerve dramatically increased the spiking activity for the FHM3 mutations.

Modeling the Whole Nerve Activity
Finally, to explore the role of Na V 1.1 channels in a more physiological environment, we modeled a WT whole nerve comprising five Aδ-fibers, five C-fibers, and ten inactive (ATP-and 5-HT-insensitive) fibers (schematically presented in Supplementary Figure 6, left). In this case, to validate the modeling results with experimental results from rodent meningeal nerves comprised of Aδ-and C-fibers in trigeminal nerves (Levy et al., 2007;Haanes and Edvinsson, 2019), we added simulated C-fibers to the model of the whole nerve, although they do not express Na V 1.1 channels, while they express P2X3 and 5-HT receptors (Zeitz et al., 2002;Yegutkin et al., 2016;Kilinc et al., 2017). Details of the C-fiber model are presented in Supplementary Figure 6. We added co-expression of P2X2 and P2X3 receptors to the whole nerve model as they are naturally expressed in trigeminal neurons of rodents (Simonetti et al., 2006). A ratio of P2X3/P2X2 was set to 75/25%, since we found in our previous study (Suleimanova et al., 2020) that this ratio closely reproduced the experimental data. The current simulated data for the whole nerve were validated with our previously published experimental results on the activity of ATP and 5-HT in mouse meningeal afferents (Koroleva et al., 2019).
For validation, we compared the inter-spike intervals (ISI) (Figure 6B) obtained from experimental data when the spiking activity in meningeal afferents was induced by 100 µM ATP (Koroleva et al., 2019) with the simulated WT model neuronal activity in afferents, also induced by 100 µM ATP (Supplementary Figure 6Aa). Likewise, we also modeled the action of 2 µM 5-HT that was also compared with experimental data with 5-HT (Koroleva et al., 2019; Supplementary Figure 6Ba). During the simulation experiments, ATP application triggered spiking which spectral profile was almost the same as in the experimental approach (Figure 6Ba). The p-value of 0.104 in the K-S test indicates the high similarity of the model and experimental results with ATP application. Likewise, with 5-HT, the model approach reproduced the data from the experiment (Figure 6Bb). The p-value of 0.128 in the K-S test indicates the high similarity of the model and experimental spikes distribution.
To reproduce the role of the FHM3 mutation in the activity of the whole nerve, we selected the L263V mutation as it produced more spikes than any other tested FHM3 mutations. In contrast to WT (Figure 6Aa), the model of the nerve with the Aδ-fiber with the L263V mutation showed a much higher neuronal activity with ATP (Figure 6Ba and Supplementary  Figure 6Ab), with 5-HT (Figure 6Bb and Supplementary  Figure 6Bb), and with co-application of ATP together with 5-HT (Figures 6Ab,Bc). Notably, with the L263V mutation most spikes were distributed with frequency higher than 10 Hz (ISI 0.1s). Thus, with a FHM3 mutation an intensive spiking activity was produced in a whole nerve model with shorter intervals which suggested a more powerful nociceptive signaling.

DISCUSSION
The main result of this study is to provide a mechanistical explanation of peripheral mechanisms of enhanced nociceptive firing activity in trigeminal neurons and the whole nerve when in silico modeling the effect of FHM3 mutations. Our computational approach provides a scientific framework for the underlying molecular mechanisms of trigeminal nociceptive firing implicated in this migraine subtype. Our in silico approach allowed us to correct ("virtually treat") the abnormal voltage characteristics of mutated channels that resulted in a significant reduction of nociceptive firing. Thus, our data suggest that compounds affecting Na V 1.1 channels in Aδ-fibers of peripheral nerves in a mutation-specific manner, may be a promising avenue for novel type analgesic anti-migraine therapy.

Role of Na V 1.1 Channels in the Activation of Aδ-Fibers
Peripheral sensory nerves express a wide range of sodium channels needed for the generation and propagation of nociceptive spikes evoked by mechanical or chemical stimulation of nerve terminals (Basbaum and Woolf, 1999;Julius and Basbaum, 2001;Giniatullin, 2020). Among other Na V subtypes, the Na V 1.1 sodium channel subtype is highly expressed in Aδ-fibers of peripheral nerves (Ho and O'Leary, 2011;Osteen et al., 2016) in addition to in central cortical interneurons (Ogiwara et al., 2007). Their abundance predicts that dysfunction of Na V 1.1 channels in Aδ-fibers should affect the transmission of nociceptive signals in a pronounced way. As Aδ-fibers are implicated in the generation of migraine pain (Melo-Carrillo et al., 2017;Haanes and Edvinsson, 2019), this provides a rationale for how pain signals are generated in trigeminal nerves in meninges where migraine pain is generated (Moskowitz, 2008;Messlinger, 2009;Zakharov et al., 2015).

Trigeminal Neuron Firing in Na V 1.1 Channels With FHM3 Mutations
Familial hemiplegic migraine type 3 (FHM3) is caused by gain-of-function mutations in the SCN1A gene that encodes the α1 subunit of voltage-gated Na V 1.1 sodium channels (Castro et al., 2009). Therefore, in our model, we implemented and tested the functional role of different gain-of-function mutations (L1670W, L263V, L1649Q, and Q1478K), which were previously shown as a genetic cause of disease in patients with FHM3 (Kahlig et al., 2008;Dhifallah et al., 2018). For comparison (as a "negative control"), we also modeled nociceptive firing in trigeminal neurons with the loss-of-function mutation (M145T) in the same α1 subunit of Na V 1.1 channels (Mantegazza et al., 2005). The strongest increase of nociceptive firing, a predictive of severe migraine pain, was observed for the L263V mutation, likely due to the increased "window current" (Figure 2Aa) because of the shift of the activation in the voltage-dependence curve to lower voltages and the inactivation curve to higher voltages comparing to WT Na V 1.1 channels (Kahlig et al., 2008). Although, experimentally, a clear trend, albeit not significant, for increase of voltage dependence of activation (V 1/2 = −21 mV for WT vs. V 1/2 = −25 mV for the L263V mutation) was reported by Kahlig et al. (2008), the mutation has a clear effect in our model (Figure 2A). It is worth noting that this property was detectable only for the L263V mutant so not with the other tested mutants (Figure 2). Our L263V finding maybe not that surprising as even a small shift in activation parameters may provide a sufficiently strong functional effect on excitability, especially near threshold voltage levels when nerve terminals are activated by physiologically relevant low levels of endogenous agonists. Consistent with this, we showed that mutant L263V required less ATP or 5-HT for induction of spiking activity (Supplementary Figure 4C), which may translate to a lower pain threshold. Of note, recently, it has been shown in a knock-in mouse model that this mutation facilitated spontaneous cortical spreading depolarization (CSD) events, the likely underlying cause of the migraine aura (Jansen et al., 2020a). However, the effect of the mutation in peripheral trigeminal neurons has not been investigated. Moreover, the observation that mutated Na V 1.1 channels can facilitate spontaneous CSD events may also have relevance to, for instance, stroke as relative peri-infarct depolarizations (PID) known to circle around the infarct core (Sukhotinsky et al., 2010), when occurring more easily, may lead to an increase in the infarct size and worse stroke outcome. Although, this has not been demonstrated in the case of FHM3, it was shown that in knock-in mice expressing Ca V 2.1 calcium channels with a FHM1 mutation in the α1 subunit (van den Maagdenberg et al., 2004), which results in a gain-of-function and hyperexcitability phenotype , the number of PIDs was greater and the infarct size was larger when infarcts were introduced experimentally (Eikermann-Haerter et al., 2012). A hyperexcitability phenotype in FHM3 is also suggested from overexpression studies in neurons (Cestele et al., 2008;Dhifallah et al., 2018). In the case of the L1670W mutant, when overexpressed in mouse neocortical neurons, a faster recovery of the mutated Na V 1.1 channel activity was observed (Dhifallah et al., 2018). Also, for the Q1478K mutant, when expressed in rat neocortical neurons, the gain-of-function leads to hyperexcitability, which, remarkably, was self-limiting meaning that for only a short time high-frequency neuronal discharges were maintained before a depolarizing block occurred (Cestele et al., 2008).
As trigger for receptor potentials, we used ATP and 5-HT, which are among the most powerful agonists of peripheral nociception in meningeal afferents (Yegutkin et al., 2016;Kilinc et al., 2017;Koroleva et al., 2019). The respective P2X3 and 5-HT3 receptors are expressed in Aδ-fibers (Ford, 2012;Kilinc et al., 2017;Sato et al., 2018). For all gain-of-function mutants, which show an increased persistent sodium current (L263V, L1670W, L1649Q, and Q1478K), we found, in our modeling, the activation of Aδ-fibers by ATP or 5-HT was enhanced. This may explain why migraine pain develops in carriers of SCN1A gain-of-function mutations. In contrast, the lossof-function M145T mutation (Mantegazza et al., 2005), which is associated with febrile seizures dramatically reduce spiking activity of peripheral Aδ-fibers of the trigeminal nerve. The finding of opposite functional outcomes on neuronal activity of the Na V 1.1 mutations depending on whether they cause FHM3 or childhood epilepsy indicates that the modeling is a useful discriminating tool to predict disease outcome.

Role of Factors Amplifying the Effect of SCN1A Mutations
Our mathematical model allows a "knock down" or artificial induction of a Na V 1 channel subtype that can be used to explore molecular mechanisms that are challenging to achieve in an experimental setting. Moreover, our model allowed dissection in a stepwise approach the contribution of several of these factors, which is highly relevant as they can potentially modify the nociceptive effect of Na V 1.1 mutations.
As a first step, we "knocked down" the contribution of the Na V 1.6 subtype to leave only the Na V 1.1 subtype to study its putative role in nociceptive signaling. Subsequently, we added the Na V 1.6 to find out that Na V 1.6, which is naturally accompanying Na V 1.1, essentially supports the pro-nociceptive role of Na V 1.1 in nerve terminals with a FHM3 SCN1A mutation.
For the second step, we added a tree structure of nerve branches to reproduce the composition of meningeal afferents more realistically, as shown experimentally (Schueler et al., 2014;Suleimanova et al., 2020) or in model (Barkai et al., 2020). The branching of a single axon combined with multiple ATP and 5-HT release events from mast cells contacting different branches can largely increase the probably of repetitive firing (Suleimanova et al., 2020). A similar amplifying role of branching was shown in the current study (Figures 4, 5). Furthermore, the slow, due to high impedance (Goldstein et al., 2019;Barkai et al., 2020) ATP/5-HT-induced receptor potential in the fine nerve terminal can initiate the persistent current through slowly inactivating Na V 1.1. and Na V 1.8 channels. A contributing role of Na V 1.8 channels to the multiple firing of primary afferents was shown previously in our model of the neuro-immune synapse in meninges (Suleimanova et al., 2020). Moreover, recently we found ATP-gated (but not 5-HT-induced) nociceptive firing not only at peripheral nerve terminals but also in more central parts of trigeminal nerve fibers in rat meninges (Gafurov et al., 2021). The latter suggests that ATP can activate nerve fibers not only at the end points of an axon but also along the fiber and that this probability is higher for ATP than for 5-HT. These factors could be essential for the increased agonist-induced multiple firing of meningeal afferents. Together, our study demonstrated that branching, multiple asynchronous ATP/5-HT release events together with persistent sodium Na V 1.1 and the presence of Na V 1.8 currents collectively enhance repetitive nociceptive firing in Aδ-fibers, most notable in FHM3-associated mutants.
For the third step, we conducted modeling experiments with 5-HT, which is a major nociceptive mediator in meninges, likely released from mast cells during migraine attacks, and that can strongly activate meningeal fibers through 5-HT3 receptors (Kilinc et al., 2017;Koroleva et al., 2019). However, unlike seen for the combined action of ATP and 5-HT, which induces very pronounced and long-term spiking activity, 5-HT alone induced firing that was not as strong as seen with ATP, likely due to a lower amplitude of 5-HT3 receptor-mediated generator potential at nerve terminals.
Notably, the concentration of ATP and 5-HT released within the meningeal neuro-immune synapse has a different time profile (Suleimanova et al., 2020). In the case of ATP, the time course is primarily determined by enzymatic hydrolysis of this endogenous purinergic transmitter by numerous ectoATPases (Yegutkin, 2008). Therefore, ATP-induced excitation depends on the spatial distribution of these enzymes in the neuro-immune synapse. Any mismatch between the ATP release site and the presence of ectoATPases in the synapse would increase the excitatory nociceptive action of this purinergic transmitter due to reduced transmitter hydrolysis as shown in the current study for the case of absence of ectoATPases activity in the meningeal synapse. Interestingly, in experimental migraine-like conditions induced by the migraine mediator CGRP, the level of endogenous ATP in the rat meninges is significantly enhanced (Yegutkin et al., 2016) by increasing purinergic drive for nociceptive excitation.
In contrast to ATP, released 5-HT is removed from the neuroimmune synapse by a relatively slow uptake (Suleimanova et al., 2020). This is specific for the 5-HT transport system and can be blocked by selective serotonin reuptake inhibitors (SSRIs), commonly used in patients with depression. The inhibition of specific transporters can increase the level of serotonin in the extracellular space and potentially amplify the pronociceptive action of 5-HT in meninges. In opposite, an increased expression of serotonin transporters is expected to provide an anti-nociceptive effect, specifically when switching off the serotonergic drive for excitation of nerve terminals via 5-HT3 receptors (Kilinc et al., 2017).
We also explored the effect of ATP and 5-HT concentrations as they could vary in migraine-relevant conditions. In the meningeal neuro-immune synapse, the concentration and the time profile for these two potential triggers of nociceptive firing is very different as determined by distinct elimination mechanisms such as enzymatic degradation and uptake, respectively (Yegutkin, 2008;Wood et al., 2014;Suleimanova et al., 2020). Indeed, we found that the shift from basal conditions with partial ATP hydrolysis and functional 5-HT uptake (Supplementary  Figures 4G,I) to complete prevention of ATP hydrolysis and removal of 5-HT uptake (Supplementary Figures 4H,J) significantly increased repetitive firing of meningeal afferents. Consistent with this, the number of spikes was also increased when enhancing the concentration of both these algogens from 2 to 10 µM, but it was more noticeable with 5-HT ( Supplementary  Figures 4C-F).
In our next step, we simulated the more realistic multifiber complex nerve model to compare the computational outcome with results from experiments performed previously (Koroleva et al., 2019). Notably, in the model of the whole nerve, ATP and 5-HT produced a powerful long-lasting firing in in analogy to results observed with firing that was induced by ATP and 5-HT in the experimental condition in mouse meninges. Most relevant to this study, nociceptive spiking of the whole nerve activity largely increased when a FHM3 mutation was introduced in the model of the Aδ-fibers, indicating that our model properly reproduces an enhanced nociceptive firing in FHM3.

Virtual Treatment in silico of Carriers of SCN1A Mutations
Since different mutations exhibit distinct changes in the voltage dependence such as activation, fast or slow inactivation (Figure 2), we next set out to correct in silico the pathological phenotype. To this end, we used as an indicator of "treatment efficiency" the reduction of spiking in meningeal afferents and compared the result with the firing in the WT model ( Figure 7G).
We observed for mutant L263V that its nociceptive phenotype was associated with abnormal changes of activation, fast and slow inactivation (Figures 7Aa,Ab,Da). The best improvement was obtained with corrected activation (Figure 7Db), since mutant L263V is activated at a more negative voltage than WT that showed in Figures 2Aa, 7Aa as a shift of the voltage-dependence activation curve to the left and shows a delayed entry into fast inactivation state. Besides, this mutation accelerated recovery, and reduced depolarizing shift in the voltage dependence of both fast and slow inactivation. However, the correction of fast and slow inactivation states in L263V less efficiently reduced spiking activity (Figures 7Dc,Dd,G). In the case of mutants L1649Q and Q1478K (Figures 7Ea,Fa), which profoundly exhibit a positive shift in fast inactivation (Figures 7Ba,Ca), the number of spikes, as expected, was most strongly decreased when we corrected the fast inactivation (Figures 7Ec,Fc,G). Since there were no significant alterations in the slow inactivation state for mutants L1649Q and Q1678K (Figures 7Bb,Cb) correction of this state did not significantly change the neuronal spiking activity (Figures 7Ed,Fd).
Our data suggest that treatment of FHM3 patients that exhibit headache as a main complaint, but due to a different Na V 1.1 mutation, may be achieved through two distinct effects of possibly different drugs, as in some cases, the most efficient correction was based on improvement of activation (L263V), whereas in other cases, the best results was obtained with correction of the inactivation state (L1649Q and Q1478K). In addition to direct correction of the genetic defect, there might be alternative pharmacological mechanisms to limit the FHM3associated enhanced repetitive firing in trigeminal neurons. One possibility may be the activation of potassium-mediated K V 7/Mcurrents, which contribute to the stabilization of the membrane potential and can be activated by the analgesic drugs, such as paracetamol (Ray et al., 2019).
In our modeling approach, we did not simulate the outcome of changes in Na V 1.1 in central neurons where these channels are expressed in interneurons (Favero et al., 2018;Sakkaki et al., 2020). Recently, computer modeling of various neuron types, similar to our whole nerve paradigm, was used to predict functional outcome of an FHM3 mutation on the central neuronal network of transgenic Dravet mice, in which Na V 1.1 channels were ablated in hippocampus and cortex (Jansen et al., 2020b). One could model the role of Na V 1.1 channels in CSD events associated with FHM3 or analogous PID in stroke to assess brain recovery after ischemia (Sukhotinsky et al., 2010). The proposed modeling approach can also actively guide the development of selective activators or inhibitors of Na V 1.1 channels to treat diseases such as hemiplegic migraine, seizures and stroke.
In summary, our data demonstrated a long-term intensive spiking activity in meningeal afferents, with various gain-offunction mutations of sodium channels associated with FHM3 in patients. Improvement of high nociceptive activity was obtained in a mutation-specific manner, being in some cases based on correction of the activation process (Figures 7Db,Eb,Fb), whereas, in others, a better result was obtained after adjustment of both the fast and slow inactivation of mutated Na V 1.1 channels (Figures 7Dc,Dd,Ec,Ed,Fc,Fd). Such modeling provides a new tool for the exploration of peripheral mechanisms in trigeminal pain and suggests that molecules reducing, probably in a mutation-specific manner, the excessive activity of Na V 1.1 channels could present a novel type of analgesic therapy for patients with migraine.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
AS and MT contributed to modeling, analysis and writing the manuscript. AMJMvdM contributed to design of the study, interpretation, writing the manuscript, and the final editing. RG contributed to the study design and supervision, writing the manuscript, and the final editing. All authors approved the final version of the manuscript.

FUNDING
The study was supported by RFBR KOMFI (Grant No. 17-00-00053) and Kazan Federal University (Grant No. 0671-2020-0059). Kazan Federal University was supported by the Russian Government Program of Competitive Growth.