Coupling of the LC-NE and cardio-respiratory systems
Figure 8 illustrates the model used in this study. We modelled the amygdala, RTN, RVLM, caudal ventrolateral medulla (CVLM) and AMB using 20, 40, 30, 30 and 30 HH-type spiking neurons, respectively. Baroreceptor-related neurons in the NTS (NTS Baro) and IE include 20 and 30 sets of excitatory and inhibitory HH-type neurons, respectively. Chemoreceptors in the NTS (NTS Chemo) include two groups of 30 HH-type neurons: one receives CO\(_2\) signals from the lungs and excites or inhibits the amygdala55,56, and the other comprises second-order peripheral chemoreceptors, which have a direct projection to the RTN, preI/I (pre-BötC) and RVLM57. Since the number of HH-type spiking neurons may influence the simulation results58, we performed parametric simulations using different number of HH-type spiking neurons. The default numbers of HH-type spiking neurons for the amygdala, RTN, RVLM, CVLM, AMB, NTS Baro, NTS chemo and IE were set to 20, 40, 30, 30, 30, 20, 30 and 30, respectively. We simulated 0.5, 1.5 and 2 times of the default numbers of HH-type spiking neurons, for instance, in case of 2 times, they were set to 40, 80, 60, 60, 60, 40, 60 and 60, respectively. The parametric simulation results demonstrate that smaller numbers of HH-neurons had different values in BR and HR in comparison with the default numbers (Supplementary Fig. S4).
Each neuron can be fundamentally represented using the following equations.
$$\begin{aligned} C \frac{dV}{dt }= & {} -\overline{g}_K m_K^4 (V-E_K) - \overline{g}_{Na} m^3_{Na} h_{Na} (V-E_{Na}) \nonumber \\{} & {} {} -g_L (V - E_L) - I - \sqrt{2 D} \xi (t), \end{aligned}$$
(1)
where C is the membrane capacitance as \(C = 36\) [pF], while \(\overline{g}_K, \overline{g}_{Na}\) and \(g_L\) are the peak conductances of potassium and sodium as well as the leak conductance (Leak), respectively. \(\overline{g}_K = 250\) [nS], \(\overline{g}_{Na} = 400\) [nS] and \(g_L = 6\) [nS], while \(E_K, E_{Na}\) and \(E_L\) represent the reversal potentials of potassium, sodium and Leak, respectively. \(E_K = -94\) [mV], \(E_{Na} = 55\) [mV] and \(E_L = -60 \pm 1.2\) [mV], provided a Gaussian distribution with a mean of − 60 and variance of 1.2. \(m_K\), \(m_{Na}\) and \(h_{Na}\) are gating parameters that depend on the membrane variable V.
For \(z = {K, Na}\),
$$\begin{aligned}{} & {} \tau _{m z} d m_z/dt = m_{\infty z}(V) - m_z, \end{aligned}$$
(2)
$$\begin{aligned}{} & {} \tau _{h z} d h_z/dt = h_{\infty z}(V) - h_z \end{aligned}$$
(3)
The constants were determined based on59, as follows:
$$\begin{aligned}{} & {} m_{\infty K} = \alpha _{\infty K} / (\alpha _{\infty K} + \beta _{\infty K}), \end{aligned}$$
(4)
$$\begin{aligned}{} & {} \alpha _{\infty K} = 0.01 (V + 44.0) / (1.0 - \exp (-(V + 44.0)/5.0)), \end{aligned}$$
(5)
$$\begin{aligned}{} & {} \beta _{\infty K} = 0.17 \exp (-(V + 49.0)/40.0), \end{aligned}$$
(6)
$$\begin{aligned}{} & {} m_{\infty N_a} = 1.0/(1.0 + \exp (-(V + 34.0)/7.8)), \end{aligned}$$
(7)
$$\begin{aligned}{} & {} h_{\infty N_a} = 1.0/(1.0 + \exp ((V + 55.0)/7.0)). \end{aligned}$$
(8)
The time constants were set as follows:
$$\begin{aligned}{} & {} \tau _{mK} = 3.5 / \cosh ((V + 40.0)/40.0), \end{aligned}$$
(9)
$$\begin{aligned}{} & {} \tau _{h Na} = 8.456 / \cosh ((V + 67.5)/12.8), \end{aligned}$$
(10)
where \(\tau _{mK}\) and \(\tau _{h Na}\) are in [ms]. \(m_{Na}\) was assumed to be \(m_{\infty N_a}\), and the time constant was assumed to be \(\tau _{m Na}=0\) based on a previous result60. \(\xi (t)\) is Gaussian noise, and D is noise strength, with \(D=20\) [pA]. Equation (1) was updated using the Euler method. In this study, the time step was set as 1.0e−4 [s].
The LC was modelled using a \(15 \times 15\) excitatory neural network with an HH-type equation, as follows:
$$\begin{aligned} C \frac{dV}{dt}= & {} -\varphi g_L(V-E_L) - \overline{g}_K m_k^4 (V-E_K) \nonumber \\{} & {} {} - \overline{g}_{Na} m_{Na}^3 h_{Na} (V-E_{Na}) - I - \sqrt{2 D} \xi (t) \nonumber \\{} & {} {} -\varphi I_{syn}^E -I_{NE}^I - I_{gap}, \end{aligned}$$
(11)
where \(I_{syn}^E\) is the excitatory synaptic input, and \(I_{NE}^I\) represents the NE interaction, as described later. \(I_{gap}\) is the electrical current based on the gap-junction from the eight neighbouring neurons and is computed as \(I_{gap}=\sum _j g_{gap} (V_i - V_j)\). \(\varphi\) is the modulating term reflecting CO\(_2\) concentration (pCO\(_2\)) in the blood, and was described according to a previous study61, as follows:
$$\begin{aligned} \varphi = 1- \frac{1}{1+(s_{1/2}/s)^{h_s}}, \end{aligned}$$
(12)
where we assumed \(s_{1/2} = 10, h_s = 5\) and \(s=pCO_2/760.0\times 10^{-2}\). We set \(C = 25\)[pF], \(\tau _{mK} = 1.75/cosh((V+40))/40)\) [ms] and \(g_{gap} = 0.005\) [nS] to reduce the firing rate of spontaneous ignition to around 1 [Hz]. Although Quintero et al.61 used an HH-type model that included \(Ca^{2+}\) channels, we used one without them for simplicity. NE emitted by the excitatory firing of LC neurons is diffused to cognitive neurons and modulates the activation gain \(g_{NE}\). The modulatory dynamics for \(g_{NE}\) are described using the second-order delay system, as follows:
$$\begin{aligned}&\tau _0 \dot{y} = -y + \sum _j \delta (t-t_j), \end{aligned}$$
(13)
$$\begin{aligned}&\tau _1 {\dot{g}}_{NE} = - g_{NE} + y, \end{aligned}$$
(14)
where \(t_{j}\) is the jth firing spike, and \(g_{NE}\) corresponds to the NE-modulated conductance released into the surroundings. \(I_{NE}^I = g_{NE}^0 g_{NE} (V - V_{syn}^I)\). \(g_{NE}^0 = 0.03\times 10^{-3}\)[nS], \(V_{syn}^I = -75\)[mV]. The rise and decay time constants of \(\tau _0\) and \(\tau _1\) are 100 [ms] and 300 [ms], respectively, as previously reported62.
The LC in the pons supplies NE input to many regions of the cortex and the respiratory centre. The diffusion of NE in the cortex is associated with arousal; an increase in NE is accompanied by an increase of arousal level. Chemosensitivity in LC neurons is activated by an increase in CO\(_2\) concentration in the blood, as the LC receives excitatory information from the respiratory CPG (central pattern generator). We found that the controller gain of the chemical feedback from the lung to the LC can influence the time-to-respond for an emotional fear- or surprise-related stimulus63. LC neurons also receive excitatory input from preI/I (pre-BötC) neurons, as previously shown in adult mice11. Since LC neuron activity is synchronized with inspiratory bursts, we modelled the connectivity from preI/I to LC with all-to-all connections. However, there is no evidence of their synaptic strengths. We set their synaptic strengths to 0.1 based on synaptic strengths for the other connectivity in respiratory CPG. In addition, the direct pathway from the LC to the RTN/pFRG (parafacial respiratory group) induces active inspiration and expiration64. Therefore, NE diffusion, that is, the activity of LC, indirectly influences the control of breathing. We included the amygdala-LC interaction in the model because fear- and stress-induced activity in the amygdala influences sensory brain regions through connections with the LC65. Furthermore, increased tonic activity in LC neurons induces anxiety-like and aversive behaviour66. However, we did not include the cortical and subcortical controls of breathing associated with emotion and cognition. Although the amygdala-LC interaction includes GABAergic inhibitory connections67 and corticotropin-releasing factor, which has an excitatory effect on LC neurons68 and adjusts the activity as well as reactivity of the LC-NE system69, we only included the excitatory interactions between the amygdala and LC in our model.
The respiratory centre and respiratory system were developed based on a mathematical model of the closed-loop system for the control of breathing proposed by41. Breathing or lung ventilation represents the exchange of air between the lungs and the surrounding environment, which takes place according to the rhythmic contraction of the inspiration and expiration neurons. The firing activities of the phrenic and abdominal motor neurons control the inspiration and expiration neurons, respectively, representing major outputs of the brainstem respiratory CPG, which generates respiratory oscillations. The closed-loop respiratory system includes two major feedback pathways from the lungs to the respiratory CPG, namely mechanical feedback and chemical feedback. Mechanical feedback is provided by pulmonary stretch receptors (PSR) located in the lungs, which transmit information regarding lung volume to the brainstem via the vagus nerve. Chemical feedback is assumed to be provided by changes in the peripheral CO2 concentration as a signal from the peripheral chemoreceptors to NTS Chemo57, which is sensitive to the levels of oxygen (O\(_2\)) and CO\(_2\) in the blood and promotes RTN activity, thereby causing an increase in HR and BR (see Fig. 8). We adjusted the NTS Chemo parameters so that active expiration was not induced by the activity of abdominal muscles during normal breathing; however, active expiration was induced during hypercapnia according to previous results57. The arterial CO2 partial pressure of around 40 mmHg was provided to central chemoreceptors of the RTN during normal breathing. The central chemoreceptors of the RTN, which are highly sensitive to the CO\(_2\) concentration in the brainstem, and thus modulate respiratory CPG in a CO\(_2\)-dependent fashion. The parameters of the respiratory centre and system were originally adjusted to the adult rat respiratory system41. We replaced mechanical system of the lung with a mathematical model of respiratory muscles, chest wall and lungs based on adult human subject data previously proposed by70. Then, we modified the parameters to adjust the adult human respiratory system as listed in Table 1.
The cardiac centre was developed based on the pathways via which the central nucleus of the amygdala (CeA) influences blood pressure during mental stress or anxiety, as proposed by72, and autonomic chronotropic control of the heart, as previously depicted73. During stressful conditions, the CeA may inhibit baroreceptive neurons in the NTS. This inhibition might switch off inhibitory inputs from the CVLM to the RVLM, which could activate neurons in the RVLM, leading to an increase in sympathetic outflow with NE modulation and an increase in blood pressure and HR. The AMB also receives excitatory information from the NTS, which increases parasympathetic outflow via acetylcholine (Ach) modulation of heart activity, thus decreasing HR. The cardiac system was represented using a cardiac muscle activation model with electrical activity based on the FitzHugh–Nagumo (FHN) neuron model74, as well as the neuromodulation of NE and Ach, which are neurotransmitters in sympathetic and parasympathetic nerves, respectively. The electrophysiological model of cardiac muscle activity is described by the following diffusion equation:
$$\begin{aligned}{} & {} \dot{v} = D \nabla ^2v + c(v(v-\alpha )(1-v)-w), \end{aligned}$$
(15)
$$\begin{aligned}{} & {} \dot{w} = b(v-dw), \end{aligned}$$
(16)
$$\begin{aligned}{} & {} c = \left\{ \begin{array}{ll} c_{1S} &{} (\dot{v} \ge 0), \\ c_{2S} &{} (\dot{v} < 0), \\ \end{array} \right. \end{aligned}$$
(17)
where v is the membrane potential of cardiac muscle cells, and w is the recovery variable. \(c_{1S} = G_c \overline{c}_{1S}, c_{2S} = G_c \overline{c}_{2S}, \overline{c}_{1S} = 1\), and \(\overline{c}_{2S} = 0.22\). \(\alpha\), b and d are parameters that define the shape of v. We hypothesised that \(G_c\) and b would be determined by the neuromodulation of NE and Ach, as shown in the following Eqs. (20, 21), based on Hill’s Eq.75.
$$\begin{aligned}{} & {} G_c = (c_{Max} - c_{Min})*H(p_{NE}) + c_{Min}, \end{aligned}$$
(18)
$$\begin{aligned}{} & {} b = (b_{Max} -b_{Min})*H(q_{Ach}) + b_{Min}. \end{aligned}$$
(19)
$$\begin{aligned}{} & {} H(p_{NE}) = 1/(1+(p_{NE_{mid}}/p_{NE})^{\alpha _{NE}}), \end{aligned}$$
(20)
$$\begin{aligned}{} & {} H(q_{Ach}) = 1-1/\left( 1+\left( q_{Ach_{mid}}/q_{Ach}\right) ^{\alpha _{Ach}}\right) . \end{aligned}$$
(21)
The HR increases with the value of parameter b and decreases with the value of b. \(H(q_{Ach})\) significantly contributes to the adjustment of HR. The average HR ranges from 64 to 69 [bpm] in adult men aged 20–60 years26, while the average breathing rate (BR) in adults at rest is 12–20 [bpm]27. We set \(p_{NE_{mid}}=350.0\), \(q_{Ach_{mid}}=20.0\), \(\alpha _{NE}=10.0\), and \(\alpha _{Ach}=5.0\) to reproduce the average values of BR and HR.
In this study, we introduced inspiratory-to-expiratory phase-spanning neurons (IE), including excitatory and inhibitory neurons, into the pons to combine the activity of the respiratory centre with that of the cardiac centre. The excitatory IE was proposed by57,76 to receive an excitatory signal from ramp-I of the respiratory centre and project to the RVLM, providing a pontine-dependent inspiratory modulation of thoracic sympathetic nerves according to a suggestion of respiratory-sympathetic coupling based on previous experimental observation77. Further, the role of the excitatory IE was confirmed with additional computational and experimental studies78. The inhibitory IE was described to receive excitatory signals from the respiratory centre, projecting into the AMB and providing respiratory-related modulation of parasympathetic nerves based on recent neurophysiological studies79,80. Here, we added an excitatory projection from ramp-I neurons in the respiratory centre to the inhibitory IE in order to control how the flow of information from the NTS to the AMB is cut off during inspiration, as per73, as well as an excitatory projection from late-E neurons to inhibitory IE to represent hypoxia-induced transient tachycardia, as reported by81.
This model has several parameters; we mostly determined the parameters based on previous studies41,60,61,62,70,71, However, we used four parameters \(p_{NE_{mid}}\), \(q_{Ach_{mid}}\), \(\alpha _{NE}\) and \(\alpha _{Ach}\) to adjust the BR and HR to be 16 and 68 [bpm], respectively, which are almost the same as those in healthy adult males in the resting state, by approximately five times of trial and error.
Experimental design
We used the developed computational model for autonomic breathing simulation under the following two types of external body vibration conditions (Fig. 1). Condition 1 included a total of 155 inputs to simulate human body vibrations via sinusoidal waves with five amplitudes (− 3, − 6, − 9, − 12 and − 15 [pA]) and 31 frequencies (0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 and 20.0 [Hz]) in a time period of 80 [s]. Condition 2 included a total of 24 inputs to simulate those with four amplitudes (− 3, − 6, − 9 and − 12 [pA]) and six frequencies (0.15, 0.25, 0.3, 1, 10 and 20 [Hz]) in a time period of 60 [min], during which simulation results were obtained at 1 [min] and every 15 [min], hereafter referred to as extraction points. As mentioned previously, vibrational inputs to mammalian bodies are known to stimulate vestibular sensory organs and then induce LC activity that is physically synchronised with the original vibrational inputs15. However, how the amplitudes of vibrational inputs are transmitted to the LC is unknown. Therefore, we used current values inputted to HH-type neurons to represent the amplitudes of the stimulations inputted to the LC. All amplitudes represent excitatory stimulation, which increased with the absolute value.
We calculated the BR, HR and NE-modulated conductance from the LC to the cortex, which represent the arousal level, activity of the inspiration and expiration neurons, as well as neuronal activity of amygdala, NTS (Baro), RTN, IE, RVLM and AMB as model outputs. Activity of the inspiration and expiration neurons ranged from 0 to 1, while neuronal activity of amygdala, NTS (Baro), RTN, IE, RVLM and AMB was calculated as neuronal spike counts. The BR, HR and NE-modulated conductance were calculated as averaged values during the period ranging from 60 to 80 [s] in the autonomic breathing simulations for condition 1 and during the period of 20 [s] before each extraction point for condition 2. We investigated the relationship between external body vibration input and model outputs of arousal level, BR, and HR using the obtained simulation results. The model was implemented using C++ language, and all simulations were performed using a computer with an Intel Xeon Gold 6242 (16C/2.8Ghz) and 384 GB Memory.
