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 “neuro-immune 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₅₀ 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 neuro-immune 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_V1 channel type (Na_V1.1, Na_V1.6, Na_V1.7), while Na_V1.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_V1.6, K_V1, K_V3, and K_V4 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_V1 channels to I_Na current (Zhang et al., 2013) and Na_V channels mRNA levels (Ho and O’Leary, 2011). The WT Na_V1.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_V1.1 and Na_V1.6, since these fibers express both type of these channels in nodes of Ranvier (Duflocq et al., 2008; Leterrier et al., 2011). We also added Na_V1.7 and Na_V1.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² (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, 2013) taking into account the contributions of Na_V1.1, Na_V1.6, and Na_V1.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_V1, 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 calcium-activated 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 steady-state 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 three-parameter Gaussian exponential function: τ = a × exp (((V–V₀)/b)²), where V is the membrane voltage, V₀ 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 half-maximal 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_Nais 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_V1.1 channels (Zheng et al., 2019) to recapitulate better the characteristics of a FHM3 mutation. To this end, we decreased the activation voltage of FHM3-associated Na_V1.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_V1.1 (Na_V1.1-L263V) channels are activated at a lower voltage than WT Na_V1.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_V1.1-L263V is wide (Figure 1Ca and Supplementary Figure 3). Second, we modeled the FHM3 L1670W mutation (Figure 1Cb), which also exhibits a gain-of-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 time constant for fast inactivation (Figures 2Ca,Cb and Supplementary 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_V1.1 channels are shown in Figure 2 and Supplementary Figures 1, 2. Using experimental data from Mantegazza et al. (2005); Kahlig et al. (2006, 2008); Cestele et al. (2008, 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 voltage-dependent 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.

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).

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_V1.1 and Na_V1.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_V1 phenotype we started with in silico “knocking down” the Na_V1.6 channel subtype in WT, whereas the conductance of Na_V1.1 was set to 0.35 S/cm². In this case, ATP only induced local receptor potential (Figures 3Aa,Ba black lines). Adding Na_V1.6, with a conductance of 0.35 S/cm² was sufficient to generate propagating spikes (Figures 3Ab,Bb black lines). Raising the activity of Na_V1.1 (conductance 0.5 S/cm²) enhanced the firing with repetitive spikes (Figures 3Ac,Bc black lines). In contrast to WT, in the FHM3 model with gain-of-function mutations L263V and L1670W, even Na_V1.1 alone was sufficiently effective to generate propagating spikes (Figures 3Aa,Ba red lines). Even higher activity was obtained when Na_V1.6 (Figures 3Ab,Bb red lines) was added or when the conductance of the Na_V1.1 channel was increased to 0.5 S/cm² (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.

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 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 was less sensitive to the absence or presence of ATP hydrolysis (Supplementary Figures 4G,H).

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 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.

Finally, to explore the role of Na_V1.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_V1.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.

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_V1.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_V1.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).

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_V1.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_V1.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_V1.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_V1.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_V2.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 (Ferrari et al., 2015), 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_V1.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 loss-of-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_V1.1 mutations depending on whether they cause FHM3 or childhood epilepsy indicates that the modeling is a useful discriminating tool to predict disease outcome.

Our mathematical model allows a “knock down” or artificial induction of a Na_V1 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_V1.1 mutations.

As a first step, we “knocked down” the contribution of the Na_V1.6 subtype to leave only the Na_V1.1 subtype to study its putative role in nociceptive signaling. Subsequently, we added the Na_V1.6 to find out that Na_V1.6, which is naturally accompanying Na_V1.1, essentially supports the pro-nociceptive role of Na_V1.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_V1.1. and Na_V1.8 channels. A contributing role of Na_V1.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_V1.1 and the presence of Na_V1.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 neuro-immune 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 pro-nociceptive 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 multi-fiber 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.

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_V1.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 FHM3-associated enhanced repetitive firing in trigeminal neurons. One possibility may be the activation of potassium-mediated K_V7/M-currents, 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_V1.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_V1.1 channels were ablated in hippocampus and cortex (Jansen et al., 2020b). One could model the role of Na_V1.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_V1.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-of-function 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_V1.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_V1.1 channels could present a novel type of analgesic therapy for patients with migraine.