Predictive modeling of hemodynamics during viscerosensory neurostimulation via neural computation mechanism in the brainstem

Lee, Jiho; Mun, Junseung; Choo, Minhye; Park, Sung-Min

doi:10.1038/s41746-025-01635-w

Download PDF

Article
Open access
Published: 23 April 2025

Predictive modeling of hemodynamics during viscerosensory neurostimulation via neural computation mechanism in the brainstem

Jiho Lee^1,2,
Junseung Mun^1,2,
Minhye Choo^1,2 &
…
Sung-Min Park ORCID: orcid.org/0000-0002-8359-8110^1,2,3,4,5

npj Digital Medicine volume 8, Article number: 220 (2025) Cite this article

616 Accesses
1 Altmetric
Metrics details

Subjects

Abstract

Neurostimulation for cardiovascular control faces challenges due to the lack of predictive modeling for stimulus-driven dynamic responses, which is crucial for precise neuromodulation via quality feedback. We address this by employing a digital twin approach that leverages computational mechanisms underlying neuro-hemodynamic responses during neurostimulation. Our results emphasize the computational role of the nucleus tractus solitarius (NTS) in the brainstem in determining these responses. The intrinsic neural circuit within the NTS harbors collective dynamics residing in a low-dimensional latent space, which effectively captures stimulus-driven hemodynamic perturbations. Building on this, we developed a digital twin framework for individually optimized predictive modeling of neuromodulatory outcomes. This framework potentially enables the design of closed-loop neurostimulation systems for precise hemodynamic control. Consequently, our digital twin based on neural computation mechanisms marks an advancement in the artificial regulation of internal organs, paving the way for precise translational medicine to treat chronic diseases.

The brain–heart axis: integrative cooperation of neural, mechanical and biochemical pathways

Article 03 March 2025

Vagal sensory neurons mediate the Bezold–Jarisch reflex and induce syncope

Article Open access 01 November 2023

Digital twin brain simulator for real-time consciousness monitoring and virtual intervention using primate electrocorticogram data

Article Open access 10 February 2025

Introduction

Recent progress in brain-computer interface (BCI) technologies has revolutionized neurostimulation practices, offering promising solutions for restoring neurological dysfunctions and facilitating brain functions. Indeed, advances made in the BCI for neurostimulation have achieved considerable success in different medical areas. In particular, applications for chronic disease treatments such as deep brain stimulation to mitigate brain disorders such as epilepsy^1,2 and spinal cord stimulation for rehabilitating motor inability following spinal cord injuries³ have been developed. Despite these achievements, the application of BCI technologies in viscerosensory neurostimulation therapies to artificially control functions of internal organs, such as hemodynamic functions of the cardiovascular system, remains unrealized to date. Although therapies have shown the potential to effectively modulate the hemodynamic functions⁴, their clinical translations have been hindered by high interpersonal variability in stimulus-driven outcomes. For example, baroreflex activation therapy designed to manage hypertension by stimulating barosensory afferents⁵ (see also Supplementary Note 1) has faced clinical limitations such as excessive blood pressure (BP) reduction and ineffectiveness, largely due to the use of open-loop stimulation modality and the heuristic selection of stimulus parameters⁶. This raised the need for a personalized computational model, a concept of digital twin (patient-specific virtual replica)⁷, as a facet of BCI technologies to individually predict stimulus-driven outcomes. From this viewpoint, the accurate predictive modeling provided as the digital twin can offer a robust framework for developing closed-loop feedback control systems⁸. However, the digital twin for the complex stimulus-driven responses, such as the nonlinear temporal evolution of blood pressure, observed during the baroreflex activation therapy^9,10, remains elusive.

A fundamental understanding of both the anatomical and computational mechanisms underlying how stimulus inputs modulate the autonomic nervous system provides an opportunity to design digital twins that can predict the stimulus-driven responses in internal organs, such as the cardiovascular system. Classical and recent studies have largely established the neuroanatomic atlas of the autonomic nervous system^11,12,13. The autonomic nervous system involves the hemodynamic regulation by relaying viscerosensory afferents from cardiovascular organs into the brainstem through vagal and sympathetic peripheral nerves¹⁴. The brainstem is an essential intermediary processor between afferent and efferent pathways. Especially, the nucleus tractus solitarius (NTS) in the brainstem plays an important role of integrating the viscerosensory information and attuning the pre-autonomic nodes such as rostral ventrolateral medulla (RVLM) and dorsal motor nucleus of the vagus (DMV)^11,13. The pre-autonomic nodes regulate the efferent pathway of the autonomic nervous system by innervating preganglionic neurons in spinal cord regions such as the intermediolateral nucleus (IML) to provide feedback on peripheral viscera¹⁵. Despite these neuroanatomic clarities, the computational mechanisms involving the stimulus-driven dynamics of neural activities and hemodynamic functions are of the uncharted unknowns. Therefore, we investigated the neural computation mechanisms embedded in the neuro-hemodynamic interactions to model stimulus-driven responses of hemodynamic functions during viscerosensory neurostimulation.

Among the aforementioned nodes in the autonomic nervous system, the NTS attracted our attention owing to its established correlation with changes in BP^9,16,17,18. This relationship between NTS neural activity and hemodynamic functions is attributed to the second-order neurons in the NTS. These neurons directly receive viscerosensory synaptic afferents from cardiovascular baroreceptors¹⁹, thereby underpinning the baroreflex and maintaining the hemodynamic homeostasis^20,21. Such predominant participation of NTS in baroreflex highlights the importance of the computational process within the NTS in transforming the barosensory inputs into the modulation of cardiovascular autonomic tones. However, the neural computation mechanisms involving how the NTS neurons orchestrate stimulus-driven hemodynamic outcomes remain unknown. While alterations in NTS neural activity induced by BP fluctuations are correlated with the BP changes owing to the role of NTS in integrating barosensory inputs^{19,22,23,24,25,26}, these insights for the afferent process do not adequately explain the efferent process from stimulus-driven changes in the NTS to the hemodynamic modulations during viscerosensory neurostimulation. This knowledge gap regarding the computational mechanism within the efferent pathway of the NTS influencing hemodynamic functions is largely due to the heterogeneous activation patterns observed among NTS neurons, even in response to identical stimuli^27,28,29. Notably, the NTS contains not only the second-order neurons directly connected to baroreceptors but also higher-order interneurons that form a recurrent network structure, capable of generating the diverse neuronal activities^19,30,31.

The network connectivity within the NTS raises a hypothesis regarding neural computation mechanism, proposing that the collective dynamics leveraging from the local circuit in the NTS neuronal population may significantly contribute to hemodynamic perturbations during the neurostimulation. This hypothesis is based on the established premise that the latent subspace of neuronal population activities, which unfolds the collective dynamics, is determined by the intrinsic network connectivity among neurons³². Both theoretical and empirical studies have demonstrated that the network connectivity can be captured by the latent space, given the close relationship between the eigenvectors of the connectivity matrix and the latent space^33,34. Furthermore, previous modeling studies have successfully replicated realistic neural circuits using this latent information^35,36,37. Consequently, this hypothesis introduces the potential for developing a digital twin to predict stimulus-driven hemodynamics by modeling realistic neural computational mechanisms within the NTS via collective dynamics embedded in latent space.

In this study, we dissect the neural computation mechanisms within the NTS by investigating the collective dynamics of baroreflex-relevant NTS neuronal populations involved in hemodynamic perturbations during viscerosensory stimulation. Specifically, we validate that the hypothesis holds under the conditions we tested, focusing on neuro-hemodynamic responses during the neurostimulation targeting barosensory afferents at the solitary tract with an electrical pulse train. We begin by demonstrating that the temporal evolution of the NTS collective dynamics exhibits a ring-shaped neural trajectory in a two dimensional (2D) latent space, which was derived through dimensionality reduction of single-neuronal activities from the NTS (Fig. 1a). These distinct trajectories for each rat could be normalized into a unit circle, and the normalized latent space was linearly coupled with stimulus-driven hemodynamics. Using these theoretical insights, we developed a digital twin framework for predictive modeling of viscerosensory stimulus-driven hemodynamic responses. Figure 1b illustrates the core concept of our modeling framework. It shows that, leveraging the neuro-hemodynamic coupling, hemodynamic perturbations are encoded into the neural trajectory within the normalized latent space. This neural trajectory can be predicted using a model of the neural circuit that reflects the NTS collective dynamics. Our approach enables us to predict changes in hemodynamics in response to stimulus inputs by decoding the predicted updates in the neural trajectory. Thus, our digital twin framework, incorporating these processes, represents the shape of the neural computation mechanism within the NTS by implementing both neuro-hemodynamic coupling and a biophysically plausible neural network model. Furthermore, the coupling between latent space and hemodynamic perturbation, which is bidirectionally switchable, enables us to individually optimize the model parameter solely using hemodynamic recordings. Therefore, this study can pave fundamental foundations for the BCI technologies toward the individually optimizable closed-loop system for viscerosensory neurostimulation to control functions of internal organs, as well as deepen our understanding of the neural bases of the autonomic nervous system to regulate the organs.

Results

Stimulus-driven collective dynamics of the NTS neuronal population in latent space

We recorded extracellular single-unit activities of NTS and simultaneously measured the femoral arterial BP of the rats (n = 10, Fig. 2a and Supplementary Fig. 1a). Electrical pulse train was delivered at solitary tracts in the brainstem projecting to the NTS as the viscerosensory stimulus inputs. The stimulation experiments consisted of 20 s of resting time (pre-stimulation), 60 s of stimulation (active stimulation), and 20 s of recovery time (post-stimulation) (Fig. 2a). During the stimulation, the temporal BP reduction was not sustained, highlighting a significant challenge in maintaining stimulus-driven BP reduction. This issue in hemodynamic control underscores the difficulties in achieving consistent and sustained BP modulation through neurostimulation, as described in “Introduction”. Only NTS neurons responsive to the neurostimulation were selected for recording (n = 192; 72% of 266 neurons from 10 rats; Fig. 1b and Supplementary Fig. 1e). The stimulus inputs resulted in temporal evolution of heterogeneous single-neuronal responses (Fig. 2b). The NTS neuronal population comprised 51% activated (n = 98), 29% deactivated (n = 55), and 20% complex-responding (n = 39) neurons of 192 total neurons (Fig. 1c and Supplementary Fig. 1f). The cross-correlations and shared variances among the observed neurons were analyzed to examine the effect of inter-neuronal dialogs on these heterogeneous stimulus-driven responses. The pairwise cross-correlations among the NTS neurons were higher than those of the dummy data. The dummy data used for control purposes consisted of randomly rearranged neuronal activities across the entire period of the stimulation experiment. The absolute values (log-scaled) of cross-correlations were as follows: 0.0533 (interquartile range [IQR]: 0.046–0.060) and 0.59 (IQR: 0.48–0.64) for the dummy and measured data, respectively, p < 1.83 × 10⁻⁴. Shared variances of each neuronal response were highly described with other neurons: 27% (IQR: 22–35%) and 83% (IQR: 76–87%) for the dummy and measured data, respectively, p < 1.83 × 10⁻⁴ (Fig. 2c). These couplings among the neurons suggest that the heterogeneous responses stem from interconnections formed within the NTS neuronal population, which supports that the NTS accompanies a self-modulating recurrent neural circuit through the interconnections. Moreover, our results are supported by previous studies that have established the existence of a local circuit within NTS shown^27,30.

**Fig. 2: Collective dynamics leveraging from the local neural circuit within the nucleus tractus solitarius (NTS).**

Our results validate that the local neural circuit within the NTS neuronal population derives the collective dynamics that simplify the complex neural activities of numerous single neurons into the dynamics of some effective variables, called latent state variables³³. This approach was based on a theoretical understanding of the dynamic properties of neural circuits^38,39. The neural circuit includes large-scale information stored in numerous single neurons⁴⁰ and thus unfolds the complex dynamics of neuronal population in a high-dimensional state space, where each neural activity of a single neuron is considered as one-dimensional (1D) dynamics. However, it has been established that the inter-neuronal couplings by the network structure in these circuits can constrain the dynamics of the population to a hidden dynamical state space with much lower dimension than that provided by the original number of neurons^{33,39,41,42,43}. These confined dynamics are referred to as the collective dynamics, and the hidden state space embedding the collective dynamics is referred to as the latent space. We confirmed that the stimulus-driven neural trajectory generated from the NTS neuronal population followed a ring topology (Fig. 2d), thus indicating that the effective latent dimension was bound to two dimensions (2D) (see detailed results in Supplementary Fig. 2, Supplementary Note 2, and Supplementary Note 4).

The collective dynamics in the latent space effectively summarize the original high-dimensional dynamics in the NTS neuronal population by effectively capturing both the information and variance inherent to the single-neuronal activities of NTS neurons. The NTS neural tuning to BP, as shown in Fig. 2e and Supplementary Fig. 3, provides insight into how activity patterns within the NTS neuronal population are influenced by changes in BP. This tuning shows that individual neurons respond to various levels of BP changes. We compared this BP tuning with the tuning to the neural state phase (θ), an internal variable representing phase shift along the ring-shaped neural trajectory. Our results show that the neural state phase offers a more nuanced discrimination of each tuning curve specific to NTS neurons, implying that collective dynamics effectively capture the characteristic actions of individual neurons. This comparison revealed that while each NTS neuron affects BP control at different levels, its activity is more closely linked to a specific collective dynamical state of the neuronal population. This finding suggests that the activity patterns of the NTS neuronal population are not directly tied to BP levels but are strongly associated with the collective dynamics. Notably, our findings indicate that the collective dynamics within the NTS neuronal population are crucial for linking the NTS neuronal activity to the resulting hemodynamic outcomes, with highlighting that the collective dynamics captures a high-resolution description of the diverse activity patterns within the NTS neuronal population. This observation was further supported by statistical analyses, which revealed that the neural state phase accounted for a significantly larger fraction of the variance and information in single-neuronal activities than did the $\Delta$BP (fraction of variance: 0.54 [IQR: 0.45–0.65] and 0.78 [IQR: 0.69–0.84] by the BP and neural state phase, respectively, p = 0.0113; mutual information: 0.82 [IQR: 0.70–0.87] and 1.0 [IQR: 0.91–1.1] by the BP and neural state phase, respectively, p = 5.83 × 10⁻⁴; Fig. 2f). Thus, our findings show that the dynamical and informational characteristics of NTS single neurons are more accurately represented through the collective dynamics rather than through external variable in physical domain such as the hemodynamic functions. This suggests that the collective dynamics may represent the properties of activity patterns of the NTS neuronal population while minimizing the loss of features of single-neuronal activities.

Neuro-hemodynamic coupling via collective dynamics in latent space

In our previous study⁹, we discovered a surprising and non-intuitive finding that the temporal dynamics of NTS neural activity have a similar tendency to the stimulus-driven $\Delta$BP during viscerosensory neurostimulation. This direct correlation between NTS activity and hemodynamic output suggested the possibility that the neural computation embedded within the NTS can importantly influence the temporal evolution of stimulus-driven hemodynamics during the neurostimulation. This hypothesis introduced a new viewpoint on the pivotal role of NTS quantitatively guiding baroreflex to modulate the circulatory system based on integrations of barosensory inputs, and thus, we elucidated the neural computational mechanism underlying the coupling between NTS neural activity and the viscerosensory stimulus-driven hemodynamic outcomes in the next series of experiments.

We first focused on the capability of single-neuronal activities to predict the $\Delta$BP during neurostimulation. The results indicated that most NTS neurons were unable to effectively predict the stimulus-driven $\Delta$BP (Fig. 3a, b). This observation is consistent with our prior results (Fig. 2) for the diverse activity patterns of NTS single neurons due to the local neural circuit, implying that only a small subset of these neurons can elicit a high correlation with the $\Delta$BP in trivial. Furthermore, such heterogeneity among NTS neurons reflects the nature of the NTS neuronal population in regulating baroreflex, where many neurons may be tuned to internal processes hidden in the NTS local circuit rather than hemodynamic functions in the external physical domain. These hidden processes suggest that while single neurons are not directly correlated to the external hemodynamics, their collective dynamics may be able to provide a more comprehensive understanding of the relationship between NTS and hemodynamics. Therefore, we analyzed the correlation between the latent space where the collective dynamics of the NTS neuronal population reside and the stimulus-driven hemodynamic perturbations (see Fig. 3, Supplementary Fig. 7, and Supplementary Note 5).

**Fig. 3: Relationship between the NTS and stimulus-driven hemodynamics.**

Our analyses were designed based on recent studies that suggested an emerging theory of “computation-through-dynamics”⁴³. This approach provided the usefulness of collective dynamics in elucidating how the brain encodes body functions by using the neural trajectory in latent space, including the rodents’ head directions encoded along the ring-shaped neural trajectory within the population of head-direction cells⁴⁴, spatial recognition encoded along the toroidal structure within the grid cells⁴⁵, and motor reactions cued by a three-dimensional neural trajectory within the cortex⁴⁶. Analogously, we assumed that the hemodynamic functions could be encoded by the latent space in the NTS neuronal population, as in cognitive and motor functions. We found that the 2D latent space was linearly coupled to the $\Delta$BP (Fig. 3c). Interestingly, each point along the ring-shaped neural trajectory (neural state phase) did not represent a unique value of $\Delta$BP, but the $\Delta$BP was represented in a linear subspace of the latent space that we refer to as “decoding space” (see “Disentanglement of neuro-hemodynamic axis” section). This linear decoding space suggests that the NTS latent space can be translated into the $\Delta$BP using a linear transformation protocol. The prediction of $\Delta$BP based on the linear protocol showed high prediction accuracy for the BP time course (89%; Fig. 3d; Supplementary Fig. 4 for the remaining animal subjects) and a lower mean squared error than did single neurons (log scale of mean squared error: −0.799 and −1.96 for single neurons and neural decoding, respectively; Fig. 3e). Therefore, our results indicate that the stimulus-driven $\Delta$BP is encoded within the collective dynamics in NTS latent space. Our statistical analyses provided an reliable evidence for this finding by presenting a much higher prediction accuracy with decoding latent space (termed the neural decoding) than the prediction accuracy with single neurons (prediction accuracy: 0.28 [IQR: 0.24–0.31] and 0.78 [IQR: 0.70–0.93] for single neurons and neural decoding, respectively, p < 1.83 × 10⁻⁴; Fig. 3f) and a smaller prediction error (log scale of mean squared error: −0.75 [IQR: −0.84 to −0.62] and −1.22 [IQR: −1.78 to −1.01] for single neurons and neural decoding, respectively, p = 2.46 × 10⁻⁴; Fig. 3g).

Inter-individual consistency of collective dynamics in NTS latent space

Although the ring-shaped topology and the linear protocol for neuro-hemodynamic coupling were common across rats as previously shown in Figs. 2 and 3, each rat subject displayed distinct neural trajectories, reflecting inter-individual variabilities both in the collective dynamics and the coupling protocol (Fig. 4a, b). Specifically, the geometrical heterogeneity of neural trajectories encompassed variations in the center point locations and the axis scales of the ring structure (Supplementary Fig. 8a). The heterogeneity in the dynamical properties encompassed differences in the positions of resting points (both resting state points before and after stimulation) (Fig. 4b) and the rotational directions across rats during neurostimulation (Supplementary Fig. 8b). Characteristic angles of the decoding vectors were also randomly distributed, implying the heterogeneity in the coupling protocol (mean ± STD of the angle of decoding vectors: 1.38${\rm{\pi }}$ ± 0.28${\rm{\pi }}$ rad; Supplementary Fig. 8c).

**Fig. 4: Normalization of different neural trajectories in NTS latent space across rats.**

Interestingly, despite the inter-individual variabilities, the rat-specific NTS latent spaces could be normalized across rats. Each neural trajectory in rat-specific latent space (${{\boldsymbol{z}}}_{{\boldsymbol{1}}}$, ${{\boldsymbol{z}}}_{{\boldsymbol{2}}}$) was normalized by aligning onto a unit circle centered at the origin (Supplementary Fig. 8d–f), and we termed the new space, including the unit circuit, as the normalized latent space (${{\boldsymbol{y}}}_{{\boldsymbol{1}}}$, ${{\boldsymbol{y}}}_{{\boldsymbol{2}}}$). Changes in the neural state phase ($\theta$) in the normalized latent space showed the increased commonality in the collective dynamics across rats. The normalization improved the commonality in the dynamical properties across rats by increasing each phase-locking value (PLV) of two distributions of stable points, one distribution for resting points before stimulation and the other one for recovering points after stimulation (resting before stimulation: from 0.26 to 0.83; recovering after stimulation: from 0.68 to 0.86), and decreasing their STD (resting: from 1.7 to 0.61; recovering: from 0.76 to 0.68; Fig. 4c). The normalized neural trajectories significantly enhanced the consistency across rats in the neuro-hemodynamic coupling by lowering the variability among rats in the latent space representation of $\varDelta$BP for each $\theta$ (n = 30 bins of $\theta$, STD of $\varDelta$BP representation for each $\theta$: 0.52 [IQR: 0.40–0.57] and 0.37 [IQR: 0.24–0.45] for non-normalized and normalized trajectories, respectively, p = 0.00258; Fig. 4d). Then, we evaluated the effect of normalization on the accuracy of $\varDelta$BP prediction by comparing the accuracy of decoding the normalized latent space to that of decoding rat-specific non-normalized latent space. We chose the across-rat common decoding space to decode the normalized latent space as the vertical axis for simple calculation (Fig. 4e and Supplementary Fig. 8e). There was no significant difference in the prediction accuracies between the non-normalized and the normalized latent spaces (p = 0.970; Fig. 4f), and this was statistically powered by an effect size of 0.0886 and a power of 0.97. These results indicate that our study had high power to detect even very small differences between normalized and non-normalized prediction accuracies. The effect size being close to zero enhanced the conclusion that there was no significant difference between the two conditions in prediction accuracy. This finding indicates that inter-individually variable hemodynamics during viscerosensory neurostimulation can be adequately predicted just with the common neuro-hemodynamic coupling protocol, which relies on the pre-determined normalized latent space.

The observed variability in stimulus-driven neural trajectories across individual rats, accompanied by the diverse prediction accuracies in Fig. 4f, raised an imperative ambiguity regarding the attributes that distinguish the individual (rat-specific) latent spaces from the normalized latent space. To address this ambiguity, we delved into the underlying factors of such variability in neural trajectory shapes and their associated prediction accuracies. One simple explanation could be that the variability might be caused by inter-individual differences in the network structure of the NTS neural circuit. However, the network structures, such as the rank of the connectivity matrix, are known to influence merely the dimensionality of the latent space³⁴. Our analyses of the NTS latent space demonstrated the consistent 2D latent space across rats, as shown in Supplementary Fig. 2, thereby rejecting the provided explanation.

Another plausible hypothesis is that while similar NTS neural circuits, corresponding to similar latent space, are shared across rats, the difference in recorded subpopulations for each rat, known as ‘recording instability’⁴⁷, might introduce the variability. To verify this hypothesis, we first showed that when subpopulations were extracted from a given neuronal population (n = 100 random subpopulations, generated by random sampling of 14 neurons from 28 neurons), these subpopulations exhibited different neural trajectories in latent space. The differences encompassed the shape distortion from the neural trajectory of original population (Supplementary Fig. 9a–c) and various prediction accuracies (Supplementary Fig. 9d). Then, we demonstrated that these heterogeneous neural trajectories could be normalized onto a unit circle while not changing their own prediction accuracy for each subpopulation (Supplementary Fig. 9e). Considering the similarity between the observed rat- and subpopulation-specific trajectories, we deemed it plausible that the heterogeneous neural trajectories did not originate from inter-individual distinct characteristics of NTS neural circuit but rather aroused from different subpopulation recordings for in each rat despite sharing similar NTS neural circuits among rats. In addition, we calculated a neural trajectory of the overall population, including all recorded neurons from different rats, and the overall trajectory still presented the 2D ring-shape (Supplementary Fig. 10a–d). This conservation of the ring topology in both the overall population and each rat-specific population is consistent with the above relationship between the original population and its subpopulations. Further, we validated that the neural trajectory in a given rat could be replicated by (properly) subsampling the overall population, not including the neurons from that rat (Supplementary Fig. 10e, f), indicating that similar single-neuronal responses were shared across rats. Therefore, these results reveal that the NTS neural circuit and latent space are consistent among individual rats, and the measured differences in neural trajectory of collective dynamics and their accuracy of neuro-hemodynamic coupling largely depend on the recording quality in the real world. In other words, we organized across-rat consistent properties of NTS collective dynamics during viscerosensory neurostimulation as follows: (1) the stimulus-driven collective dynamics follows ring-shaped neural trajectory in 2D latent space; (2) stimulus-driven hemodynamic outcomes can be represented on across-rat normalized latent space; (3) neural computation mechanisms embedded within the NTS neural circuit are similar across rats.

Predictive modeling of stimulus-driven hemodynamics based on neural computation mechanism

Our findings regarding the neural computation mechanism in the NTS enabled us to identify the dynamical system underlying the hemodynamic perturbation during viscerosensory neurostimulation and thus develop a digital twin to replicate the neuro-hemodynamic responses. Specifically, our digital twin was designed as a predictive modeling framework that predicts the stimulus-driven neuro-hemodynamic outcomes based on an NTS neural circuit model. The circuit model reflected collective dynamics in NTS latent space, and thus embedded the neuro-hemodynamic coupling protocol, enabling a linear transform from the collective dynamics into hemodynamics. We named this digital twin framework the hemodynamics model based on barosensory input-driven neural dynamics (H-BIND; Figs. 1b and 6a). The prediction of hemodynamic perturbation by the H-BIND starts with transforming a current hemodynamic value (BP_t) into a state in the modeled normalized NTS latent space (${\vec{{\boldsymbol{y}}{{^{\prime} }}}}_{t}$) (process of “Encoding”). Then, the H-BIND calculates the next latent state (${\vec{{\boldsymbol{y}}{{^{\prime} }}}}_{t+1}$) using the artificial neural network, reflecting the NTS neural circuit model embedding the collective dynamics (process of ‘Stimulus-driven dynamics calculation’). Finally, the H-BIND returns the updated hemodynamic value (BP_t+1) by inverse-transforming the updated latent state (process of ‘Decoding’). For this H-BIND design, the system identification process was subdivided to construct the NTS neural circuit model and the neuro-hemodynamic encoder/decoder.

The NTS neural circuit in H-BIND used the structure of a recurrent neural network (RNN)³⁵. This was based on previous modeling studies that have resembled the realistic network connectivity pattern by capturing the latent space with the RNN-based model, effectively mimicking the structure of actual neural circuit^35,36,37. This is with a rationale that the latent space results from the inherent network connectivity in neuronal population³², thus suggesting that characteristics of network connectivity are captured by the latent space. Specifically, the eigenvector space of the connectivity matrix of the neural population relates to the latent space^33,34. This aligns with our observation that an activity pattern simulated by a trained RNN in our neural circuit model could reproduce the activity pattern in NTS recordings (Fig. 5). We first posited that the collective dynamics leveraged by the RNN could be characterized as a linear combination of individual unit activities. This assumption of linearization was supported by our empirical finding, which demonstrated that linear regressions of NTS single-neuronal activities effectively captured stimulus-driven changes in each reduced dimension of NTS latent space (see Fig. 5a, b, and Supplementary Note 6). The RNN formulated with the linear relationship between single-neuronal activities and latent space successfully reproduced the stimulus-driven neural trajectory in NTS latent space, and the pattern of single-neuronal activities by RNN units could reflect the observed pattern of single-neuronal activities of NTS (Fig. 5c). Indeed, the trained RNN successfully modeled the collective dynamics of the real NTS neural circuit, replicating both the measured single-neuronal activities and neural trajectory (Fig. 5d–g). For the neuro-hemodynamic encoder/decoder, we used the bi-directional property of the linear neuro-hemodynamic coupling protocol, implying that the stimulus-driven hemodynamic perturbations can be transformed into the stimulus-driven neural trajectory in the normalized latent space and vice versa (see “H-BIND” section). Moreover, the bi-directionality of neuro-hemodynamic coupling enabled us to train the RNN without neural recordings by estimating the ground-truth neural trajectory from hemodynamics (Fig. 6c, e), while avoiding distortion of the neural trajectory in empirical neural recording (Supplementary Fig. 9b, c). This aligns with our observation that the neural trajectory reconstructed from hemodynamics could successfully capture the ground-truth neural trajectory (see Supplementary Note 7). Thus, the H-BIND can be individually optimized only with hemodynamic recordings, providing the advantages of not requiring severely invasive NTS recordings and minimizing the aforementioned recording quality issue that involves neural trajectory distortion and low prediction accuracy.

**Fig. 5: Modeling NTS neural circuit using neural trajectory measures.**

**Fig. 6: Predictive modeling of stimulus-driven hemodynamics using neural computation mechanisms within the NTS.**

We evaluated the performance of H-BIND to predict stimulus-driven hemodynamics by comparing its prediction accuracy to that of decoding the measured NTS activities (Fig. 6b and Supplementary Fig. 13). In three different rats each, the individually optimized H-BIND accurately predicted the hemodynamic perturbations, in contrast with the relatively lower and inter-subject variable accuracies with the neural decoding. The diminished performance in neural decoding is plausibly attributable to the poor recording quality (see “Inter-individual consistency of collective dynamics in NTS latent space” section). Further, the reconstructed neural trajectories from hemodynamic recordings, based on the bi-directional neuro-hemodynamic coupling, were substantially different with the trajectory from neural recordings but captured by the H-BIND (Fig. 6c). These highlight superiorities of the H-BIND over the direct neural decoding not only in the prediction performance but also in the reliability to reflect latent dynamics. Moreover, we performed statistical analyses after constructing a total of 10 H-BIND models individualized for each rat (n = 10). Our statistical analyses support the above representative results (prediction accuracy: 0.83 [IQR: 0.74–0.88] and 0.95 [IQR: 0.93–0.96] for empirical neural decoding and H-BIND, respectively, p = 0.00195; prediction error: 0.17 [IQR: 0.11–0.20] and 0.036 [IQR: 0.024–0.041] for neural decoding and H-BIND, respectively, p = 0.00195; Fig. 6d) and neural trajectory prediction (phase error: 0.13 rad [IQR: 0.075–0.17] and 0.054 rad [IQR: 0.040–0.087] for neural decoding and H-BIND, respectively, p = 0.0371; radial error: 0.017 [IQR: 0.013–0.024] and 0.00072 [IQR: 0.00041–0.0012] for neural decoding and H-BIND, respectively, p = 0.00195; Fig. 6e).

Next, utilizing a trained H-BIND, we propose a design of a closed-loop system for precise control of hemodynamic outcomes during viscerosensory neurostimulation, and provide in-silico simulations for its proof-of-concept. For the control strategy, we employed nonlinear model-predictive control (NMPC), a method that optimizes the stimulation by predicting the stimulus-driven effects in a few future steps based on precise nonlinear predictive models⁴⁸. Our controller design for NMPC was configured to receive the state in 2D latent space of RNN in H-BIND as input and provide feedback on stimulation as output (Fig. 7a). The NMPC requires the objective function (cost function) that quantifies the control performance to optimize the stimulation feedback. We chose the 2D Euclidean distance between the target and the current states in the latent space as our cost matrix. This approach of using latent space for monitoring the stimulus-driven responses significantly reduces computational burden, bypassing the need to consider the entire units in the RNN of H-BIND. Our closed-loop system design capitalizes on the ability of NTS neural computation mechanism to simplify the complex stimulus-driven system into manageable dynamics within 2D latent space, thus enabling to design of a computationally efficient control system. To pinpoint an appropriate target within the latent space for the controller, we analyzed the phase portrait of the latent space, selecting a stable node (not the resting point) among equilibria as the target (Fig. 7b). The developed control system successfully directed the state within the latent space towards the target, achieving precise control of stimulus-driven response of $\Delta$BP through the exploitation of bi-directionality in the neuro-hemodynamic coupling (Fig. 7c–e). In contrast, the uncontrolled open-loop stimulation was unable to guide neither the latent states nor the $\Delta$BP response towards the targets. Statistical analyses affirm the control performance of our developed system (n = 30 repeated simulations; mean ± STD of log scale cost in latent space: −1.61 ± 0.182; mean ± STD of cost in $\Delta$BP: 0.2596 ± 0.0884 mmHg; Fig. 7f). Further, as illustrated in Fig. 7d, our control simulation achieved approximate continuous decrease of −10 mmHg in the BP compared to the open-loop condition, suggesting its clinical potential for modulate the BP (see Supplementary Note 1).

**Fig. 7: Design of closed-loop system for precise control of hemodynamic response during viscerosensory stimulation and its proof-of-concept in in-silico simulation.**

We further compared the dynamical characteristics at both the stable and unstable targets (Supplementary Fig. 14a). At the unstable point, the dynamical velocity long the ${y{\prime} }_{1}$ axis remains consistently positive, regardless of the stimulus frequency, indicating that the dynamical system does not allow to maintain the target. In contrast, at the stable point, the velocity can be directed in any orientation depending on the stimulus frequency, allowing for sustained control to remain at the target. Moreover, when evaluating control performance in targeting the unstable versus the stable points (Supplementary Fig. 14b), our results show that targeting the unstable point not only failed to maintain the target but was also less effective in prolonging BP reduction compared to targeting the stable point. Thus, we anticipate that our approach prioritizing stable BP reduction aligns with crucial clinical needs in hemodynamic control via neurostimulation.

Discussion

This study first elucidated the neural computation mechanisms embedded within the NTS in the brainstem, pivotal for understanding hemodynamic perturbations during viscerosensory neurostimulation, then introduced a digital twin framework to predict the stimulus-driven hemodynamic outcomes by replicating the computational mechanisms. We discovered that the collective dynamics hidden in the NTS neuronal population produce a ring-shaped neural trajectory in 2D latent state space during the stimulation, with this latent space linearly representing the stimulus-driven hemodynamics. This neuro-hemodynamic coupling through the NTS latent space demonstrated that the computational mechanism, unfolded by the collective dynamics of the NTS neuronal population, played a dominant role in determining the stimulation outcomes. These findings provided motifs to develop the digital twin framework (H-BIND) to accurately predict and model the viscerosensory stimulus-driven hemodynamics by replicating the computational mechanisms within the NTS. The framework employed the RNN to realistically replicate the NTS neural circuit embedding the collective dynamics in latent space, and this implied that prediction for the stimulus-driven hemodynamic outcome was performed by translating changes in the latent space into hemodynamic perturbations, using the neuro-hemodynamic coupling protocol (also a linear transformation, see “Neuro-hemodynamic coupling via collective dynamics in latent space” section). Our digital twin approach based on neural computation mechanism in the autonomic nervous system, particularly in the brainstem, has the potential to facilitate the development of a closed-loop neurostimulation system for controlling hemodynamics, as shown in Fig. 7.

Closed-loop systems for precise neurostimulation have shown outperforming clinical efficacies over open-loop systems in treating various chronic disorders, such as drug-resistant epilepsy². However, we emphasize that their use of on-off controllers, while suitable for episodic symptoms such as seizure attacks, is inadequate for the continuous regulation of visceral functions, such as hemodynamics, which is essential for treating chronic disorders affecting internal organs. Modeling the dynamical system of ‘plant’ (equal term of dynamical system in control engineering) can provide a robust framework for designing precise and robust controllers⁸ for these visceral functions. The model of plant quantitatively predicts the complex responses of system variables within the plant to stimulus inputs⁴⁹, thereby enabling us to design a biomarker for feedback in a closed-loop system. These model-based control systems, such as NMPC, have been established to excel in enhancing or individually optimizing stimulus-driven responses^50,51,52. In this context, our digital twin approach could provide a successful solution to design a model-based control system for closed-loop neurostimulation. Our results show that each individualized H-BIND captured diverse stimulus-driven responses across rats (Fig. 6), and this demonstrated that our modeling framework satisfies the concept of the digital twin. Further, we show that the latent space embedded by the RNN in the digital twin framework could be utilized as a feedback biomarker for the NMPC system design. Our NMPC simulations based on the latent space of H-BIND effectively controlled the stimulus-driven BP responses (Fig. 7). These results provide evidence for the clinical promise of our digital twin framework in paving the way for developing closed-loop viscerosensory neurostimulation, suggesting a personalized solution for chronic cardiovascular disorders such as hypertension⁴. This digital twin framework may be extensively applicable to other various visceral functions, such as respiratory and digestive motility. This is due to neuroanatomical similarity across visceral functions in the aspect that the NTS receives its viscerosensory afferents^14,53,54. Additionally, neurostimulation modalities such as vagus nerve stimulation, which engages neuroanatomical pathways similar to those examined in this study, present a potential to seamlessly apply our approach for precise neurostimulation (for further details, see Supplementary Note 10). Therefore, we envision that this study, which provides a novel design of digital twin for precise neurostimulation based on neural computational mechanisms, suggests a new paradigm of the clinical BCI technologies for artificial regulation of internal organs, potentially advancing toward future personalized medicine treating chronic diseases.

Meanwhile, regarding the clinical relevance of our example closed-loop system design, we acknowledge that clinically optimal BP targets may be unstable, which is not supported by our design that targets only stable points. However, unstable BP targets fundamentally face a significant challenge in maintaining a reduced BP after reaching a target point. Achieving stability at such targets often requires more complex control strategies, and in some cases, may even be unachievable if the system’s dynamics do not support stability at the desired point. In contrast, our approach robustly elongates BP reduction by focusing on stable points. While this may represent a compromise in terms of not achieving the most ambitious clinical outcomes, our method prioritizes ensuring stable BP reduction, which is crucial for effective hemodynamic control via neurostimulation. Indeed, our simulation to compare stable and unstable targets show that at the unstable point, the dynamical system cannot maintain the target and the control performance was less effective in prolonging the BP reduction compared to targeting the stable point (Supplementary Fig. 14). Therefore, we expect that this focus on stability addresses a critical need and provides a practical solution for hemodynamics management.

The proposed H-BIND as a digital twin framework provides significant benefits for the real-world application of a closed-loop neurostimulation control system based on this model. First, the developed framework leverages the neural computation mechanisms in the NTS to reduce the high-dimensional neuro-hemodynamic response into 1D dynamics along the circumference within a 2D latent space. This reduced representation enables the 2D Euclidean distance as a biomarker for stimulation feedback, thus facilitating a computationally efficient cost function for the control system. Second, hemodynamic responses reliably capture the collective dynamics within NTS as shown in Supplementary Fig. 11, thus, the bi-directionality in neuro-hemodynamic coupling between the NTS latent space and hemodynamics enables us to train the H-BIND only with hemodynamics (especially, training with the neural trajectory reconstructed from the hemodynamic data). This training approach avoids the risk of neural trajectory distortions due to the recording instability (Supplementary Fig. 9), thereby effectively incorporating the neural computation mechanisms underlying stimulus-driven outcomes. These benefits of our digital twin approach present an advance over current methods for closed-loop neurostimulation. Traditional data-driven predictive modeling methods, such as linear state space model⁵⁰, have faced a challenge in elucidating the underlying computational mechanisms, complicating the understanding of actual neural dynamics. In the context of baroreceptor activation therapy, achieving precise control through predictive modeling still remains elusive, with existing controllers requiring intensive manual tuning for parameters in response to BP changes⁵⁵. In comparison, our modeling method not only provides insights into the inherent neural dynamics during neurostimulation but also supports efficient controllers tailored to the individually optimizable predictive model. Therefore, we anticipate that our digital twin approach can suggest a promising solution to the challenges of precise closed-loop neurostimulation.

However, the current H-BIND framework is a computationally heavy model owing to calculations of dynamics in each unit in the RNN to model complex dynamics in the NTS neural circuit. We expect that this limitation of computational burden may be improved by approximating the RNN structure into the combinations of multiple linear state space models⁵⁶, not only simplifying dynamics calculation but also enabling the design of computationally efficient controllers such as the linear quadratic regulator⁵⁰. In addition, modeling the 1D variable, such as BP, with the RNN that incorporates a 2D latent space can fundamentally provide more accurate predictions than directly modeling 2D variables, such as neural activities. This raises the need for further investigation into the capability of our framework to estimate naive neural activities effectively.

As for the clinical application of this study, in vivo validation for the closed-loop system using the developed predictive model is essential. This process will require overcoming significant technical challenges (see Supplementary Note 11). Further, the current modeling framework does not account for effects from complex interactions with other organ systems, such as the respiratory and endocrine systems, which may alter the dynamics depending on their state. For example, the cardio-respiratory interaction demonstrates that respiration frequency influences heart rate through the NTS, which integrates afferents from pulmonary stretch receptors, baroreceptors, and chemoreceptors to detect arterial CO₂ concentration^25,57. Renal fluid mechanisms offer slower but more substantial feedback on BP compared to the baroreflex¹⁴, implying the mechanical difference between short-term and long-term hemodynamic regulation. These physiological factors can significantly modulate the collective dynamics of the NTS and its relationship with hemodynamics. Another aspect is the physiological state of subjects under anesthesia, which could limit the interpretability of our results to the awake subjects. Furthermore, regarding the translation of our modeling framework from rodent models to humans, it needs to address the challenges associated with interspecies physiological differences. Potential strategies include employing large animal models such as macaque, which more closely resemble human physiology^37,50,58, and utilizing non-invasive recording techniques such as functional magnetic resonance imaging (fMRI) to capture NTS neural activities²⁴. These strategies are expected to enable calibration of model parameters based on data from large animal studies and human imaging studies to ensure accurate predictions in human subjects. Thus, further studies need to explore how to improve the developed model to cover these physiological factors (see also Supplementary Notes 8 and 9). This would broaden the applicability and efficacy of our modeling framework into diverse clinical scenarios.

Methods

Animals

All surgical and experimental protocols were approved by the Institutional Animal Care and Use Committee at Pohang University of Science and Technology (POSTECH-IACUC; approval numbers: POSTECH-2021-0017 and POSTECH-2022-0020) and followed the animal experiment guidelines from the NIH Guide for the Care and Use of Laboratory Animals. We used the adult male Sprague-Dawley rats (age: 8–10 weeks; weight: 270–320 g; provided by Orient Bio, Rep. of Korea, colony reference at Charles River Laboratories, and by Koatech, Rep. of Korea, colony reference at Envigo). We performed surgeries with the rat subjects for non-survival acute experiments, and thus, at the end of the experiments, the subjects were euthanized using CO₂ gas.

Surgery

Before the surgical procedures, the rat subjects were anesthetized with intraperitoneal injection of customized urethane solution (injection does: 1.2–1.4 g/kg; density: 0.2 g/mL diluted in saline solution; urethane: purchased from Merck (Sigma–Aldrich), product code U2500) after initial induction using isoflurane (5% in O₂ flow rat 0.4–0.6 L/min; 5–10 min). Supplemental doses of 0.3–0.6 g/kg were additionally injected if needed. After anesthesia, the left femoral artery was cannulated with a polyethylene tube (SP-10, Natsume Seisakusho Co.) for arterial BP recording. The subjects underwent the surgeries to expose the brainstem as described in previous studies^9,17,59. Tracheal intubation was administered only if the respiration of the subject was not stable. The subjects were placed on stereotaxic frame (68002, RWD Life Science), and the craniotomy was performed to plant the reference and ground anchor screws for electrophysiological recording. Thereafter, the placement of the subject was changed to permit sharp flexion of the neck for a surgery to approach the dorsal surface of the medulla in the brainstem. Dorsal neck muscles were incised and retracted to expose the atlanto-occipital membrane. The membrane was incised, and part of the occipital bone was removed to clearly expose the rhomboid fossa.

Solitary tract stimulation

Our stimulation protocol followed the method described in our previous study⁹. In brief, we targeted the solitary tracts in the NTS⁶⁰ to electrically simulate the barosensory afferents in the NTS. The electrical stimuli were delivered with an arbitrary waveform generator (IZ2-32, Tucker-Davis Technologies (TDT)), which is electrically isolated by using a battery (LM48M-250, TDT). We used the customized bipolar electrode fabricated with two enamel-encapsulated stainless steel wires (diameter: 150 µm, GoodFellow). The electrode was fixed at the stereotaxic arm with a supporter structure made in custom. Stimulation used the bipolar square waveform (pulse width: 100 µs; frequency: 20 Hz; intensity: 150–250 µA, differently chosen along individual subjects as inducing hemodynamic responses such as the decrease of blood pressure). Although the stimulus frequency is a significant factor in determining the properties in temporal dynamics of NTS neural activity^16,61 and even hemodynamics, especially the blood pressure¹⁰, we chose the 20 Hz based on a previous study showing that electrical stimulation varying around 20 Hz could effectively restore impaired baroreflex⁵⁵. The stimulation was applied for 60 s to evoke sufficient dynamic perturbations both in the NTS neural activities and hemodynamic functions. We ensured for 10 min before stimulation experiments to stabilize and thus prevent unknown fluctuations in neuro-hemodynamic signals. For the stimulation, we placed a bipolar electrode on the solitary tracts ipsilateral to the NTS to emulate viscerosensory afferent inputs. This approach was guided by a range of in vivo and in vitro studies^27,61,62,63. We opted to target solitary tract rather than peripheral baroreceptors due to the neuroanatomical evidence showing that the NTS is monosynaptically connected to peripheral baroreceptors through a single-neuronal fiber¹⁹. This choice allowed us to efficiently innervate the barosensory afferents near the NTS while avoiding the risk associated with additional surgery required to expose the peripheral baroreceptors. We specifically targeted the solitary tracts adjacent to the dorsomedial NTS, known to be primary sensory hubs for baroreflex afferents⁵⁹.

Meanwhile, natural stimuli to the NTS comprise the integration of diverse firing patterns from numerous baroreceptor afferent fibers, which also covary with BP change. This complexity results in precisely replicating the natural stimuli through artificial electrical stimulation being extremely challenging or nearly impossible. Nonetheless, our stimulation protocol was designed to induce a dynamic response in the NTS that approximates natural conditions, though in a more limited form. Our supplementary results support this approach, demonstrating that the stimulus-driven NTS dynamic responses were similar to those observed with peripheral baroreceptor activation during BP elevation, as shown in Supplementary Fig. 5 (see also Supplementary Note 3). This suggests that while our method may not perfectly reproduce natural stimuli, it effectively elicits a relevant and comparable NTS response similar to peripheral baroreceptor activation.

Electrophysiological neural recording

After exposing the surgery of the dorsal surface of the brainstem, we recorded the NTS neural activities following the methods shown in our previous study⁹. We visually chose the obex (the causal end point of quill-nib shaped calamus scriptorius, the inferior part of rhomboid fossa) as the anatomical reference for stereotaxic coordinates. The dorsomedial NTS, which is known as barosensory neurons gathers⁵⁹, was targeted with regard to the reference (rostral: 0.3–0.6 mm; lateral: 0.4–0.8 mm; ventral: 0.3–0.8 mm) based on the widely known anatomy of the rat brain⁶⁴. Extracellular single-neuronal spikes were observed to monitor neural activity of the NTS with a 16-channel silicone probe microelectrode (A1x16-Poly2s-5mm-50s-177-A16, NeuroNexus Technologies). Raw analog signals were converted into digital signals with a sampling rate of about 24 kHz (exactly 24414.0625 Hz) by a commercialized DAQ platform, which consists of a microprobe-headstage adapter (ZCA-DIP16, TDT), a headstage (ZC16, TDT), and a neuro-digitizer (PZ5-32, TDT). The DAQ platform was protected from electromagnetic noise, especially power-source 60 Hz noise and ground noise, by shielding our recording environment with a customized Faraday’s cage (copper-mesh and aluminum profiles). The electrical ground of the DAQ system was connected to the cage, and the DAQ output was transferred to a PC-controlled bioprocessor (RZ5, TDT) outside the cage through a shielded optical cable. Real-time digital filters were applied on the neural data using a software suite (OpenEx, TDT): High-pass filter (>0.5 Hz) was applied to reject any long-term baseline fluctuations; a Notch filter at 60 Hz was additionally applied to compensate for the analog canceling of electrical noise. The recorded neural activities were first saved as TDT tank format and then transformed into MATLAB (Mathworks) data format “.mat” by using TDT’s offline MatlabSDK.

Hemodynamic recordings

The real-time beat-to-beat arterial blood pressure was monitored through the intravenous catheter in the femoral artery using a commercialized recording system (IX-RA-834, iWorx), which was administered by its own software (LabScribe, iWorx), and a pressure transducer (BP-102, iWorx) placed at the level of the subject's heart. The catheter tubing was threaded out of the Faraday cage and connected to the transducer, and this was to avoid electrical noise from hemodynamic recording devices from interfering with the neural recording. The blood pressure was recorded in 100 samples/s, and the event points, such as stimulation on/offset, were manually marked on the time course. The transducer and catheter were filled with heparinized saline (heparin sodium: purchased from Merck(Sigma–Aldrich), product number: H5515-25KU; solution density: 100 IU/mL) to prevent blood coagulation that severely attenuates the pressure. Additionally, the blood diffused out into the catheter was perfused back into the artery using the solution about once an hour as well. The recorded raw pressure data was exported into the MATLAB data format “.mat”.

We calculated hemodynamic functions of the mean blood pressure, the heart rate, and the heart rate variability by post-processing the pressure recordings. The mean blood pressure (normally referred to as the arterial BP in previous sections) was calculated by bandpass filtering the time course of raw pressure (high: >0.001 Hz and low: <0.5 Hz). To synchronize with neural data, we discretized the mean blood pressure in time steps of 0.5–0.75 s, empirically chosen for animals (rat1:0.6 s; rat2: 0.5 s; rat3: 0.75 s; rat4: 0.6 s; rat5: 0.6 s; rat6: 0.6 s; rat7: 0.6 s; rat8: 0.6 s; rat9: 0.6 s; rat10: 0.6 s). The heart rate was calculated from the detected systolic peaks in the high-pass-filtered raw pressure signal (>1Hz). We counted the number of beats in each time sampling step (0.01 s) as the instantaneous heart rate and elicited the heart rate by low-pass filtering the instantaneous heart rate with a moving average filter (window: 3 s). As for the heart rate variability, we used the spectral power of the heart rate. We calculated the spectrogram of heart rate for each 15 s time step with the Hanning window and 80% overlap among adjacent time steps. The low-frequency power of heart rate variability was the cumulative spectral power in the frequency range <0.15 Hz and the high-frequency power of heart rate variability in the frequency range 0.15–3 Hz (we empirically chose the ranges to clearly distinguish the difference responses between low and high-frequency powers with considering various ranges reported in previous studies^65,66,67).

Stimulus artifact rejection

The stimulus artifacts were removed following the methodology outlined in our previous study⁹, which was guided by earlier research⁵⁰. First, template subtraction was used to address the long tail of stimulus artifact. The template was created by averaging the tails of hundreds of artifacts from 0.8 to 12.3 ms after the stimulus spike, which was then subtracted from each raw data. Subsequently, the ‘discard-and-treat method’ was employed to reject the short-term and high-intensity stimulus spike. This method involves discarding the instantaneous spike, recovering the missing data by pasting the segments immediately before and after the artifact. Specifically, we discarded ±1.2 ms of stimulus spike, and then filled in the missing with adjacent data. After removing artifacts, the data were bandpass filtered (300–6000 Hz).

Identification of single neurons

We post-processed the recorded NTS neural activities offline to detect the single-neuronal spikes and appropriately cluster them to identify each single neuron (the process called ‘sorting’^68,69). Before the sorting, we first removed stimulus artifacts from the neural data (see “Stimulus artifact rejection” section). We performed the sorting using an open-source Klusta software suite⁷⁰. The spikes were automatically detected from the filtered neural data (300–6000 Hz) and were clustered into putative single neurons. The clusters were manually corrected if needed, using an open-source Phy Python package⁷¹. The toolkit provided the autocorrelogram of clusters that were calculated as the histogram of temporal lags around each spike (1 ms bin-step within 20 ms window range) and thus represented the refractory period of a single neuron. We considered the cluster having peak values within the range of 2 ms as the capture of multi-unit activities rather than single-neuronal activity, and excluded the clusters. Additionally, we excluded the clusters that did not have the V-shaped spike waveform or spread out the spikes at more than 10 channels, thereby considering as a mis-detection of noise.

Classification of NTS neurons

We divided the recorded and sorted NTS neurons into the responsive and non-responsive neurons to the solitary tract stimulation. The responsive neurons were chosen as showing stably identical responses to repeated stimulation trials. We evaluated the repeatability by using the phase-locking value calculated with a firing rate of a neuron as $\left|\sum _{k}\frac{1}{N}{e}^{j({x}_{k}-{x}_{{ref}})}\right|$ where N is the time length, x_k is the firing rate at time k, and x_ref is the reference signal at time k. For simplicity of calculation, we used a cosine wave as the reference signal that has the same period as the time course of repeated stimulation. A threshold of PLV > 0.1 was empirically chosen to determine whether the neuron shows a repeatable response. Interspike-interval (ISI) validated that the stimuli evoke the spikes. We presented the probability distribution of ISI and found peaks synchronized to the stimulus interval of 50 ms. We excluded non-responsive neurons and utilized only responsive neurons in further analyses. Moreover, we classified the responsive neurons into three subtypes with consideration of their heterogeneous responses to stimulation input: activated, deactivated, and complex. The activated neurons were ones showing the increased firing rate during stimulation, the deactivated neurons were ones showing the decreased firing rate during stimulation, and the complex neurons were ones showing both the increase and the decrease of firing rate during the stimulation.

Analysis of single-neuronal interactions

To evaluate the joint activity induced by the interconnection among the NTS neurons, we calculated the pairwise cross-correlation and the shared variance⁷². The correlations were reported with the Pearson’s correlation coefficient by using MATLAB built-in function ‘corrcoef’ with inputs of firing rates of each neuron. Although the correlation coefficient elicited both positive and negative value due to heterogeneous responses of NTS neurons (see “Classification of NTS neurons” section), we took the absolute value of correlation coefficients to only assess the strength of inter-neuronal coupling and prevent the mean across neurons being 0 that cannot be distinguished to the case of no coupling.

The shared variances were reported from the Factor analysis by using MATLAB built-in function ‘factoran’:

$${\bf{X}}={{{\mathbf {M}}}}+{\boldsymbol{\Lambda }}{\bf{F}}+{\bf{E}}$$

(1)

where ${\bf{X}}$ is the k–by–n observed firing rate matrix (k: number of observations; n: number of neurons), ${\bf{M}}$ is the k–by–n observed mean matrix whereby column vectors are identical, ${\boldsymbol{\Lambda }}$ is the k–by–d factor loading matrix (d: number of common factors), ${\bf{F}}$ is the d–by–n common factor matrix, and ${\bf{E}}$ is the k–by–n independent specific factor matrix. The covariance matrix ${cov}({\bf{X}}-{\bf{M}})$ can be expressed as: ${cov}\left({\bf{X}}-{\bf{M}}\right)={\bf{L}}{{\bf{L}}}^{{\bf{T}}}+{cov}\left({\bf{E}}\right).$ This elicits that the total variance is equal to the summation of the specific variance and the factor loading, and thus the shared variances of neurons ${\boldsymbol{i}}$ equal to the fraction of the factor loading relative to the total variance was calculated as:

$${\rm{diag}}({\bf{L}}{{\bf{L}}}^{{\bf{T}}},{\boldsymbol{i}})/{\rm{diag}}({\bf{L}}{{\bf{L}}}^{{\bf{T}}}+{cov}\left({\bf{E}}\right),{\boldsymbol{i}})$$

(2)

where ${\rm{diag}}({{\cdot }},{\boldsymbol{i}})$ is ${{\boldsymbol{i}}}^{th}$ diagonal component of matrix. The shared variance stands for the fraction of variance of a neuron explained by the other neurons. Both the (absolute) correlation coefficient and the shared variance were compared to those of shuffled dummy data, which were expected to have no entanglement. The shuffled dummies for neurons were generated by a random permutation of the time courses of firing rate.

Dimensionality reduction for collective dynamics

The joint activities across single NTS neurons confined the degree of freedom in neuronal population dynamics, and this characterized that the n-dimensional state space of neuronal population (n is the number of neurons) could be compressed into a low-dimensional subspace that effectively embeds the neural states for collective dynamics. The inherent noise in neural state representations in high dimensions makes it hard to understand the properties of collective dynamics^45,73, and thus, the redundant subspaces were removed by reducing the dimensions. We used a nonlinear dimensionality reduction technique, called ‘Isomap’⁷⁴, as previously described⁴⁴ because it has been known to stably and effectively estimate the ‘neural trajectory’ in latent space compared to the linear approach, such as principal component analysis. Before the dimensionality reduction, the firing rates of neurons were estimated by counting the number of spikes in the time step (0.5–0.75 s, empirically chosen for animals, see also “Hemodynamic recordings” section) and then filtering with the Gaussian kernel (window: 3 s; standard deviation: equal to time step). The firing rate patterns were normalized across stimulation trials by dividing the original firing rates by the average value of maximums across stimulation trials (${{nFR}}_{i,t}=F{R}_{i,t}/\mathop{{\rm{avg}}}\limits_{i}[\mathop{\max }\limits_{t}F{R}_{i,t}]$ where $i$ is for the ${i}^{{th}}$ stimulation trial, $t$ is the time-point, ${FR}$ is the original firing rate, and ${nFR}$ is the normalized firing rate). We performed the dimensionality reduction by using the ‘isomap’ class in the Scikit-learn Python package (neighbor number: 5). The Isomap used the geodesic distance matrix of the NTS neuronal neural state representation ns in high-dimensional coordinates, and then achieved the reduced dimensions of the latent coordinate by calculating eigenvectors. In addition, the amplitude in reduced dimension basically reflects the firing rate, as the original high-dimensional space, prior to dimensionality reduction, represents the firing of individual neurons.

Persistent homology

We used the persistent homology^44,45 in topological data analysis to identify the geometric characteristics of neural data cloud in the latent space. This topological data analysis approach has been successfully used to reveal neural geometries such as the sphere in primary visual cortex⁷⁵, the ring of head-direction cells⁴⁴, and the toroid of grid cells⁴⁵. Thus, we applied it to discovering the topology embedded within the NTS. The Ripser Python package⁷⁶ was implemented to calculate the persistent homology. The calculation of persistent homology starts with defining a geometrical space, where an axis of the space represents changes in firing rate of an NTS neuron, meaning that the space is an N-dimensional space according to N NTS neurons. This space is generally called the neural (state) space. Then, the neural space enables us to present multiple activities of NTS neurons at a time-point as a spatial point in the space. Thus, the time courses of activities of NTS neurons are represented as dynamical trajectories in space, called the neural trajectory. This suggests that the dimensionality reduction is the process to extract a low-dimensional subspace, efficiently capturing the neural trajectory, from the original high-dimensional neural state space, and the extracted space is the latent space. Meanwhile, for real data, the neural recordings are digital due to discrete time steps, and thus the trajectory should also be a point cloud. In this context, the persistent homology is the process to quantitatively describe the topology of a data point cloud in the latent space. We considered each point of NTS neural activity in the latent space as the 0-diameter balls. If we gradually increase the radius of each ball, the overlaps among the balls occur. Drawing connection lines between each pair of points which their balls are overlapped, the connections could make closed-structures such as ring and hollow sphere at a range of radius scale (Supplementary Fig. 2a). These structures topologically designate ‘hole’ and the essence of calculation was to analyze the geometrical ‘holes’ (for example, if a structure shapes ring, this implies there exists 1-dimensional hole). In formal notations, the persistent homology resulted in the properties of topological holes in Vietoris-Rips complexes⁷⁷ (formal term for the closed structure we use here). The properties of holes included the dimension and the scale range where holes persist. The dimension of holes was sub-classified into three different homology groups^78,79: 0-dimensional group H₀ (no holes; ‘point’); 1-dimensional group H₁ (a hole made by single line; ‘ring’), and 2-dimensional group H₂ (a hallow space enclosed by single surface; including examples of ‘hallow sphere’ or ‘toroid’). The scale range is referred to as ‘lifetime’. The longer lifetime implies the stable existence of a hole, and thus confirms the topological characteristics. The Betti barcode graph visualized the lifetimes of holes and thus enabled us to determine the geometrical shape of the NTS neural trajectory.

Effective dimension of latent space

Despite the latent space embedding by the dimensionality reduction (see “Dimensionality reduction for collective dynamics” section), this providing of nonlinear eigenvectors (reduced dimensions) of neural data did not specify the optimal dimension that captures the neural trajectory in latent space; this is because the maximum number of eigenvectors is theoretically equal to the high dimension of original coordinate. We achieved the effective dimension of NTS neural trajectory in the latent space by analyzing the variance explained by each reduced dimension and the correlation dimension of the latent space. The explained variance along the order of reduced dimensions was computed as the percentage of variance of the neural trajectory being projected on a reduced dimension relative to the total variance of the neural trajectory. The correlation dimension was calculated, as previously described^80,81, by analyzing the statistics of nearest neighbors in the neural data points. We acquired the geodesic distance matrix of the latent space. The cumulative number of neighbors ‘$m$’ within a given distance threshold ‘$d$’ was counted while gradually increasing the distance threshold. The correlation integral was achieved as $C\left(d\right)\approx m/{N}^{2}$ where $N$ is the total number of data points. The correlation dimension was estimated with the slope of linear regression in the log-log graph between the correlation integral and the $d$ such that $C\left(d\right)=k{d}^{r}$ where $r$ is the correlation dimension and $k$ is a constant.

Single-neuronal dependence on internal variables

We analyzed the capture of single-neuronal activities by the external variable in the physical domain (especially, the change in blood pressure) and the internal variable tracking along the neural trajectory. We used the internal variable defined as the ‘neural state phase’: the angle of a point along the neural trajectory spanning on latent space, adopting the polar coordinate system with the center point of the trajectory as the reference origin. The capture of single-neuronal activities by the variables was quantitatively assessed with the fraction of variance explained by the tuning curve of neurons along the external and the interval variables and the mutual information between the single-neuronal activities and the variables. The tuning curves of single neurons along the external and the internal variable were calculated as the mean firing rate of the neuron at a given value in 25 bins for each variable. The fraction of variance for each neuron was calculated as the ratio of variance for the changes in tuning curve relative to the total variance of single-neuronal activities. The mutual information for each neuron was calculated as previously described⁸². We analogously used this spatial information formula to compute the mutual information (MI) of a neuron for the external and the internal variables:

$${\boldsymbol{MI}}=\mathop{\sum }\limits_{i=1}^{M}{p}_{i}\left(\frac{{f}_{i}}{\bar{f}}\right){\log }_{2}\left(\frac{{f}_{i}}{\bar{f}}\right)$$

(3)

where ${p}_{i}$ is the occupation ratio to be in ith bin of the variable, M is the total number of bins, ${f}_{i}$ is the firing rate of the neuron at the i^th bin, and $\bar{f}$ is the mean firing rate of the neuron for overall bins.

Disentanglement of neuro-hemodynamic axis

We analyzed the entanglement of input-driven hemodynamics into the latent space representations and compared it with the entanglement into the single-neuronal activities. To assess the entanglement, we performed the prediction of the input-driven hemodynamics with the neural data and used the prediction accuracy as the strength of entanglement. We used linear regression for the prediction in both cases (including a constant term as the bias) since decoding neuronal activities into the physical domain generally follows linear relationships³⁵. The regression was performed with a MATLAB built-in function ‘regress’ with a 95% confidence interval. The measure of prediction accuracy was assessed using the goodness-of-fit of regression, specifically employing the coefficient of determination to quantitatively provide how successfully the neural data fit into the hemodynamics. The coefficient of determination of linear regression is equal to the square of Pearson’s correlation coefficient, which has been widely used for the measure of accuracy in neuroscience and neural engineering applications, such as modeling the stimulus-induced brainwave⁵⁰ and decoding the behavior from cortex recording⁵⁸. Even the coefficient of determination quantifies the fraction of variance in ground-truth inputs explained by the prediction outputs, thereby representing the measure of robustness. Additionally, we quantified the error of prediction by calculating the squared errors between the measured hemodynamic perturbations and the regressed predictions.

For single-neuronal entanglement, the prediction accuracy of single neurons was calculated from a single-variable linear regression between the time course of firing rate and the hemodynamic perturbations.

As for the latent space representations, we used two-variable linear regression of the 1st and 2nd reduced dimensions into the hemodynamic perturbations (since we confined the effective dimension of latent space as 2D, as shown in Fig. 2 and Supplementary Fig. 2). This allowed to define ‘decoding vector’ consisting of the regression coefficients, and correspondingly ‘decoding space’, a 1D linear subspace (line subspace), consisting of the projections of collective neuronal neural state along the neural trajectory onto the decoding vector. The projection of the states onto the decoding vector designated the hemodynamic values. This could be expressed as follows:

$$\begin{array}{c}\Delta H=\left[\begin{array}{cc}a & b\end{array}\right]\cdot \left[\begin{array}{c}{Z}_{1}\\ {Z}_{2}\end{array}\right]+{c}_{{bias}}\\ =a{Z}_{1}+b{Z}_{2}+{c}_{{bias}}\end{array}$$

(4)

where $\Delta H$ is the change in hemodynamic functions such as blood pressure and heart rate, ${\left[\begin{array}{cc}a & b\end{array}\right]}^{{\bf{T}}}$ is the decoding vector (and the decoding space should be the $y=\frac{b}{a}{x}$ line subspace in latent space), ${\left[\begin{array}{cc}{Z}_{1} & {Z}_{2}\end{array}\right]}^{{\bf{T}}}$ is the neural state along the neural trajectory in 2D latent space, and ${c}_{{bias}}$ is a constant bias (The projection of neural state ${\left[\begin{array}{cc}{Z}_{1}^{0} & {Z}_{2}^{0}\end{array}\right]}^{{\bf{T}}}$ corresponding to $\Delta H=0$ onto the decoding vector is not the origin of decoding space, and this generally means that $a{Z}_{1}^{0}+b{Z}_{2}^{0}$ is not 0. Thus, we added the bias to adjust the origin).

Normalization of neural trajectory

The inter-individually heterogeneous neural trajectories, as shown in Fig. 4a, b and Supplementary Fig. 6, were normalized by aligning the neural trajectories with combinations of linear transformation; ‘sliding’, ‘scaling’, ‘flipping’, and ‘rotating’ (Supplementary Fig. 6d) (thus, the combined transform was defined as a linear mapping function; ${f}_{{align}}{{:}}{\mathbb{Z}}{{\to }}{\mathbb{Y}}$ where ${\mathbb{Z}}$ is original latent space spanning heterogeneous neural trajectories and ${\mathbb{Y}}$ is aligned(normalized) latent space).

The sliding adjusted the center points (${\vec{Z}}_{{center}}$) of ring-shaped neural trajectories into the origin point of latent space by subtracting the vector of the center point from the neural trajectory; $\vec{Z}-{\vec{Z}}_{{center}}$. We estimated the center points as the mean of the maximum and the minimum values along each horizontal and vertical axis of the latent space. The scaling weakly transformed the ring trajectories into the unit circle by normalizing each scale of the axes; $S\cdot \vec{Z}$ where $S=\left[\begin{array}{cc}1/{s}_{1} & 0\\ 0 & 1/{s}_{2}\end{array}\right]$, ${s}_{1}$ is the horizontal scale, and ${s}_{2}$ is the vertical scale. The scales of axes were estimated by the half of the length between the maximum and the minimum values. These established a consistent geometrical reference frame of trajectories. The flipping and the rotating were to normalize the dynamical traces and the neuro-hemodynamic coupling protocols. The flipping unified the directions of stimulus-driven temporal rotation along the trajectory into counter-clockwise, the direction of increasing the angle of neural state phase; $F\cdot \vec{Z}$ where $F=\left[\begin{array}{cc}-1 & 0\\ 0 & 1\end{array}\right]$. The rotating transformation involved in adjusting the decoding vector to be equal to the vertical axis; $R\cdot \vec{Z}$ where $R=\left[\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right]$ where $\theta =\frac{{\rm{\pi }}}{2}-{\tan }^{-1}\left(\frac{b}{a}\right)$ and $\frac{b}{a}$ is the slope of the decoding vector (see “Disentanglement of neuro-hemodynamic axis” section). Following these, the normalized neural trajectory (${\left[\begin{array}{cc}{{\boldsymbol{y}}}_{{\boldsymbol{1}}} & {{\boldsymbol{y}}}_{{\boldsymbol{2}}}\end{array}\right]}^{{\bf{T}}}{\mathbb{\in }}{\mathbb{Y}}$) possessed the across-rat common protocol of neuro-hemodynamic coupling:

$$\Delta H=\left[\begin{array}{cc}0 & 1\end{array}\right]\cdot \left[\begin{array}{c}{{\boldsymbol{y}}}_{{\bf{1}}}\\ {{\boldsymbol{y}}}_{{\bf{2}}}\end{array}\right]={{\boldsymbol{y}}}_{{\bf{2}}}$$

(5)

where ${\left[\begin{array}{cc}0 & 1\end{array}\right]}^{{\bf{T}}}$ is the normalized decoding vector, and the values in the vertical axis should be equal to the hemodynamic perturbations that were elicited as $a{Z}_{1}+b{Z}_{2}+{c}_{{bias}}$ with unaligned and rat-specific neural trajectory. Additionally, the normalized coupling protocol exhibited that the normalized decoding vector is a unit vector of the vertical axis and that the removal of bias. The sliding allowed the zero-point of ring-shaped neural trajectory that corresponds to no change in hemodynamics to be the origin point of latent space and this rejected bias; The scaling normalized the scale of neural trajectory along vertical axis in range from −1 to +1 and this matched in one-to-one relationship between the vertical value and the hemodynamic perturbations (which were normalized in the range from −1 to +1).

To assess the normalization across animal subjects, we measured the PLV and STD of the polar distribution of neural state phase ($\theta$). We made a polar distribution histogram spanning 30 bins, and compared the histogram between the unaligned and the aligned (normalized). The PLV and STD were measured for resting and recovery status that the neural state in the latent space should be static and stable. The PLV was calculated by:

$$\left|\frac{1}{N}\mathop{\sum }\limits_{k=1}^{N}{e}^{j\left({\theta }_{k}-\bar{\theta }\right)}\right|$$

(6)

where $N$ is total number of $\theta$ data, ${\theta }_{k}$ is the ${k}^{{th}}\theta$. The $\bar{\theta }$ is polar mean of $\theta$ distribution calculated by:

$$\frac{1}{N}\mathop{\sum }\limits_{k=1}^{N}{\tan }^{-1}\frac{\sin \left({\theta }_{k}\right)}{\cos ({\theta }_{k})}$$

(7)

and the circular STD was calculated as

$$\sqrt{-\log \left(\frac{1}{N}\mathop{\sum }\limits_{k=1}^{N}{(\sin \left({\theta }_{k}\right)}^{2}+{\cos \left({\theta }_{k}\right)}^{2})\right)}$$

(8)

Modeling the local neural circuit in NTS

To model the local neural circuit within the NTS, we designed the model structure with an artificial neural network, especially the RNN, which has been widely used to model the biophysically plausible modeling of dynamics in neuronal population^{36,38,39,42,83}. The RNN used in our study is designed to handle continuous time and employs a sigmoidal activation function, distinguishing it from discrete applications in machine learning³⁵. This design is fundamentally based on the Wilson–Cowan model, which itself mathematically represents a synaptic connection between neurons^84,85. Our model design is as follows:

$$\begin{array}{c}\tau \frac{d\vec{{\bf{x}}}\left(t\right)}{{dt}}=-\vec{{\bf{x}}}\left(t\right){\boldsymbol{+}}{\bf{W}}\sigma \left(\vec{{\bf{x}}}\left(t\right)\right)+\overrightarrow{{{\bf{w}}}_{{\boldsymbol{u}}}}u\left(t\right)+\vec{{\boldsymbol{\varepsilon }}}\left(t\right),\\ \vec{{\bf{r}}}(t)=\sigma (\vec{{\bf{x}}}\left(t\right))\end{array}$$

(9)

where $\vec{{\bf{x}}}\left(t\right){\boldsymbol{\in }}{{\mathbb{R}}}^{N}$ is the vector of membrane currents of N neurons (also called ‘activation’), $\vec{{\bf{r}}}$ is the firing rates of neurons, $\sigma \left(\cdot \right)$ is the activation function of $\tan h\left(\cdot \right)$, $\tau$ is the time constant of neuron, ${\bf{W}}$ is the N-by-N connectivity matrix, $\vec{{{\bf{w}}}_{{\boldsymbol{u}}}}$ is the synaptic weight vector of stimuli inputs on neurons, $u\left(t\right)$ is the time course of visceral sensation-mimicking stimulation, and $\vec{{\boldsymbol{\varepsilon }}}\left(t\right)$ is the vector of independent Gaussian noise with variance of 1/N. We used the difference equation of the RNN model by discretizing with the time step $\Delta t$ equal to $\tau =100$ ms:

$${\vec{{\bf{x}}}}_{t+1}={\bf{W}}\sigma \left({\vec{{\bf{x}}}}_{t}\right)+\overrightarrow{{{\bf{w}}}_{{\boldsymbol{u}}}}{u}_{t}+{\vec{{\boldsymbol{\varepsilon }}}}_{t},$$

(10)

(and equivalently derived as ${\vec{{\bf{r}}}}_{t+1}=\sigma \left({{\bf{W}}\vec{{\bf{r}}}}_{t}+\overrightarrow{{{\bf{w}}}_{{\boldsymbol{u}}}}{u}_{t}+{\vec{{\boldsymbol{\varepsilon }}}}_{t}\right)$ in the firing rate form⁸⁵).

We embedded the latent representation into the RNN model as previously described⁸⁶: $\vec{{\bf{y}}}={{\bf{W}}}_{{\boldsymbol{P}}}^{{\bf{T}}}\vec{{\bf{r}}}$ and ${\bf{W}}\vec{{\bf{r}}}={{\bf{W}}}_{{\boldsymbol{FB}}}\vec{{\bf{y}}}$ where $\vec{{\bf{y}}}$ is the vector of neural state long the neural trajectory in the normalized latent space, ${{\bf{W}}}_{{\boldsymbol{P}}}$ is the N-by-2 projection matrix from firing rate to the latent space, and ${{\bf{W}}}_{{\boldsymbol{FB}}}$ is the N-by-2 recurrent feedback matrix that transforms the latent space into the synaptic recurrent inputs of neurons. This derived our final model equations:

$$\begin{array}{c}{\vec{{\bf{x}}}}_{t+1}={{\bf{W}}}_{{\boldsymbol{FB}}}{\vec{{\bf{y}}}}_{t}+\overrightarrow{{{\bf{w}}}_{{\boldsymbol{u}}}}{u}_{t}+{\vec{{\boldsymbol{\varepsilon }}}}_{t},\\ {\vec{{\bf{r}}}}_{t}=\sigma \left({\vec{{\bf{x}}}}_{t}\right),\\ {\vec{{\bf{y}}}}_{t}={{\bf{W}}}_{{\boldsymbol{P}}}^{{\bf{T}}}{\vec{{\bf{r}}}}_{t}.\end{array}$$

(11)

Before training the RNN structure of the NTS circuit model, the activation vector ${\vec{{\bf{x}}}}_{t=0}$ and the N-by-2 matrix ${{\bf{W}}}_{{\boldsymbol{P}}}$ were initialized as 0. The vector $\vec{{{\bf{w}}}_{{\boldsymbol{u}}}}$ and the N-by-2 matrix ${{\bf{W}}}_{{\boldsymbol{FB}}}$ were pre-determined by random generation from a Gaussian distribution. Specifically, each element of $\overrightarrow{{{\bf{w}}}_{{\boldsymbol{u}}}}$ was drawn from a Gaussian distribution with zero mean and unit standard deviation, and each column of ${{\bf{W}}}_{{\boldsymbol{FB}}}$ was drawn independently from the Gaussian distribution. We chose not to optimize but pre-determine the $\overrightarrow{{{\bf{w}}}_{{\boldsymbol{u}}}}$ and ${{\bf{W}}}_{{\boldsymbol{FB}}}$ based on our practical need and empirical evidence. Previous studies have generally pre-determined these parameters rather than optimization based on training. For instance, Darshan et al.⁸⁶ used identical and deterministic sinusoidal distributions for $\overrightarrow{{{\bf{w}}}_{u}}$ and ${{\bf{W}}}_{{FB}}$, while Dubreuil et al. ⁸³ optimized ${{\bf{W}}}_{{FB}}$ but fixed $\overrightarrow{{{\bf{w}}}_{u}}$ as a Gaussian distribution with zero mean and unit standard deviation. Russo et al.⁸⁷ used ${{\bf{W}}}_{{FB}}$ of a uniform distribution ranging from −0.5 to 0.5. Motivated on these, we empirically set the $\overrightarrow{{{\bf{w}}}_{u}}$ and ${{\bf{W}}}_{{FB}}$ to provide a practical and effective solution for our model. In addition, including them as part of the model training would significantly increase the complexity of the model. Given the limited data size from our in vivo experiments, optimizing them could lead to overfitting or training failure, thus, we aimed to balance the model complexity and data availability.

For the training, we optimized ${{\bf{W}}}_{{\boldsymbol{P}}}$ using the recursive least squares (RLS) approach⁸⁸ (RLS hyperparameters: constant parameter = 50, epoch iteration = 20, number of RNN units = 500, noise addition = Gaussian distribution with 0 mean and 0.04 standard deviation). The procedures of training started to calculate ${\vec{{\bf{r}}}}_{t}$ as $\sigma \left({{\bf{W}}}_{{\boldsymbol{FB}}}{\vec{{\bf{y}}}}_{t-1}+\overrightarrow{{{\bf{w}}}_{{\boldsymbol{u}}}}{u}_{t-1}+{\vec{{\boldsymbol{\varepsilon }}}}_{t-1}\right)$, and then compared the estimation of neural trajectory ${\widehat{{\vec{{\bf{y}}}}}_{t}}={{{\bf{W}}}_{{\boldsymbol{P}}}^{{\bf{T}}}}_{\left({\boldsymbol{t}}-{\boldsymbol{1}}\right)}{\vec{{\bf{r}}}}_{t}$ with the ground-truth neural trajectory ${\vec{{\bf{y}}}}_{t}$. Finally, the ${{{\bf{W}}}_{{\boldsymbol{P}}}^{{\bf{T}}}}_{\left({\boldsymbol{t}}\right)}$ was updated using the estimation error ${\vec{{\bf{y}}}}_{t}-{\widehat{{\vec{{\bf{y}}}}}_{t}}$⁸². The training data was the time course of ${\vec{{\bf{y}}}}_{t}$ representing the neural trajectory in the 2D normalized latent space. We reformed the time course of a single stimulation trial (basically, 20 s resting time, 60 s stimulation, and 20 s recovery time) for the training data. During the training, the parameter fitting was repeated for at least 20 epochs. This helped the precise fitting since the RLS algorithm progresses along the time axis and could result in more accurate fitting with longer time courses. After the training finished, we validated the fitted RNN by performing the feedforward simulation and verifying the accurate match between the simulated neural trajectory and the training data. The simulation was serially repeated 5 times, and we verified the successful prediction for all repeats.

H-BIND

To infer the viscerosensory input-driven hemodynamic perturbations, we used the NTS neural circuit model and the coupling protocol of the neuro-hemodynamic axis. The NTS model calculated the input-driven collective dynamics within the NTS, and the dynamics were transformed into the hemodynamic representations using the coupling protocol. This procedure mathematically quantified the hemodynamics after visceral sensation, and thus, the H-BIND identified the hemodynamic system using the neural computation mechanisms within the NTS.

The model architecture of H-BIND consisted of ‘encoder’, ‘NTS neural circuit model’, and ‘decoder’ (Supplementary Fig. 12a). First, the encoder transformed the changes in hemodynamic functions into the drift within the normalized latent space as known that the hemodynamic perturbations correspond to the neural trajectory along the unit circle in the normalized latent space. For the encoding, we first calculated the current phase ($\theta$) of neural state in the latent space as ${\sin }^{-1}\Delta H$ (where $\Delta H$ is the change in hemodynamic, normalized from −1 to 1), and then obtained the horizontal value of neural state as $\cos \theta$; its vertical value was equal to the $\Delta H$. Although ${\sin }^{-1}\Delta H$ could generate two possible angles, the $\theta$ was selected to satisfy monotonic increase since the neural trajectory only turns in a counter-clockwise direction in our normalized latent space (see “Normalization of neural trajectory” section). This was based on our observations that the neural trajectories exhibited monotonic rotations as illustrated in Supplementary Fig. 8b. The monotonic rotations were normalized into a counter-clockwise direction without loss of generality. This counter-clockwise rotation implied the monotonic increase of $\theta$ in normalized latent space. Moreover, we used the minimum point of $\Delta H$, which is −1, as the reference of our encoding. This minimum point was mapped to [0, −1] in the normalized latent space, which corresponds to $\theta =-\frac{\pi }{2}$. Following this approach, we assigned $\theta$-values in the left half-plane for $\Delta H$ occurring before the minimum point and $\theta$-values in the right half-plane for $\Delta H$ occurring after the minimum point. Next, after encoding the hemodynamic state into the latent state, the neural trajectory was obtained by calculating the temporal drift of neural state using the NTS neural circuit model. The model received both the current neural state and the stimulation as inputs, and then computed the new neural state in the next time step. The neural trajectory was simply transformed back into the change in hemodynamics by the decoder since the vertical values of the neural trajectory are equal to the hemodynamic perturbations. Thus, this allowed us to predict the temporal dynamics of hemodynamics induced by stimulation (Fig. 5a).

The training of H-BIND identically means the training of the NTS neural circuit model (since the encoding/decoding protocols are based on the determined neuro-hemodynamic coupling rules), and this required the neural trajectory as training data (see “Modeling local neural circuit in NTS” section). However, we did not use the neural trajectory extracted from neuronal recordings, and this was due to their fundamental challenge in stably capturing the collective dynamics of neuronal populations. The neural recordings in the real world capture only a part of the neuronal population rather than the whole neurons in population (the ‘recording instability’ described in the “Inter-individual consistency of collective dynamics in NTS latent space” section). This capture of subpopulation largely distorted the neural trajectory as shown in Supplementary Fig. 9. The latent dynamics distortion would weaken the neuro-hemodynamic coupling and thus result in inaccurate prediction apart from actual hemodynamics. Indeed, our results showed that the direct decoding of neural recordings failed to estimate accurate hemodynamics (Supplementary Fig. 9d, e). Consequently, we decided to train the H-BIND with the reconstructed (encoded) neural trajectory from hemodynamic recordings. This decision was based on the bidirectionally switchable coupling between dynamics in (normalized) latent space and hemodynamics (see “Neuro-hemodynamic coupling via collective dynamics in latent space” section). Further, encoding hemodynamics into neural trajectory via the coupling accurately captured the latent dynamics by avoiding the risk of capturing only subpopulation (see Supplementary Fig. 11 and Supplementary Note 7), thus enabling quality training of the H-BIND.

We assessed the performance of H-BIND to predict the stimulus-driven hemodynamics using the k-fold cross-validation method⁵⁰. This approach enables a thorough assessment of model performance by partitioning the dataset into a train and test subset. For example, in a 4-fold cross-validation, the dataset is divided into 4 subsets. 3 segments are used for training, and the remaining one serves as the test subset. This process was repeated 4 times, with each segment taking turns as the test part. This method provides a robust estimate of model performance, leveraging the available data effectively while reducing the risk of overfitting. Our choice of the k-fold cross-validation was based on the small size of our dataset, which consisted of repeated stimulation trials with a maximum of 4 trials per rat. After completing the cross-validation, the final model accuracy was obtained by averaging the results across all folds. Additionally, we evaluated the H-BIND’s performance in predicting neural trajectories by examining both radial and phasic accuracies.

Nonlinear model-predictive controller

We aimed to control the collective dynamics within the NTS and correspondingly the stimulus input-driven BP response during viscerosensory neurostimulation. The collective dynamics reside in the 2D latent space, and the 2D Euclidean distance between the current state and the target state in the latent space was used as the objective (cost) function for the control. We chose a stable node but not a resting point as the target point. We defined an objective function to articulate the cost (or loss) between current and target points in the latent space as follows (Fig. 7a):

$$L\left(\overrightarrow{{\boldsymbol{y}}{{^{\prime} }}}\right)={{||}\Delta \overrightarrow{{\boldsymbol{y}}{{^{\prime} }}}{||}}_{2}$$

(12)

where $\vec{{\boldsymbol{y}}{\boldsymbol{{\prime} }}}=\left[\begin{array}{c}{{\boldsymbol{y}}{\boldsymbol{{\prime} }}}_{{\boldsymbol{1}}}\\ {{\boldsymbol{y}}{\boldsymbol{{\prime} }}}_{{\boldsymbol{2}}}\end{array}\right]$. This geometrical approach simplified the complex effect of stimulus inputs on the NTS neuronal population and lightened the computational complexity to determine the feedback on stimulation. Thus, the control problem could be simplified into solving the optimal solution of the stimulation pattern in a two-variable dynamical system with this cost metric of 2D distance.

To design a controller based on the dynamical system and the cost, we employed a nonlinear model-predictive controller (NMPC)⁴⁸. It enables us to handle the state variables in a closed loop by predicting the expected cost for a given input. For the NMPC, we chose the prediction horizon of 5 steps, indicating the controller determines an optimal input based on the expected cost in a 5-step future. This determined the total cost as the sum of each cost in the prediction horizon. Using the total cost, the NMPC primarily calculated the series of optimal inputs for each step in the prediction horizon. However, it actually applied the solution of the first step and resolved the optimization again for the next step. This means that the NMPC determined only the current input value with future predictions. Moreover, the constraints of the control problem were bound to the range of input, the stimulus frequency, from 0 to 1. The final optimization problem was as follows:

$$\begin{array}{c}{u}_{{opt}}\left(t\right)=\mathop{{\rm{argmin}}}\limits_{{\boldsymbol{u}}}\left(\mathop{\sum} \limits_{i}L\left(\vec{{\boldsymbol{y}}{{^{\prime} }}}\left(t+i\right)\right)\right)\\ {\bf{s}}{\boldsymbol{.}}{\bf{t}}{\boldsymbol{.}}\,0\le u\le 1\end{array}$$

(13)

The controller solved the optimization problem that minimizes the cost in the prediction horizon by using the interior point method⁴⁸, realized by the built-in function ‘fmincon’ in MATLAB (Mathworks).

Data analysis and statistics

Resources included all measured neurons that satisfied the classification criteria (see “Identification of single neurons” section) from 10 different rat subjects. No statistical standard was used to determine the sample size for subjects and neurons. The experimental protocols of the study did not require any experimental or control groups of animals, thereby not including random allocation and blinding tests. The statistical tests were all two-sided. This required the normality test with the Kolmogorov–Smirnov method. With resources rejecting the null hypothesis of the normality test that the resources follow normal distribution, the Mann–Whitney U-test was performed; actually, all cases we analyzed rejected the null hypothesis. When failing to reject the null hypothesis of the Mann–Whitney U-test, we performed statistical power analysis using G*Power (version 3.1.9.7). Two-sided resources, including parametric statistics, usually used boxplots with MATLAB built-in function ‘boxplot’ with raw data scatters: cross-correlation and shared variance (measure vs. shuffle dummy as in Fig. 2b); fraction of variance and mutual information (external vs. internal variable as in Fig. 2f); prediction accuracy (single neuron vs. neural trajectory as in Fig. 3f, g; non-normalized vs. normalized as in Fig. 4d, f; empirical decoding vs. H-BIND as in Fig. 6d, e); if not the case, we specified so.

Data availability

Data used in this project are available at https://osf.io/z7pj8/.

Code availability

Custom scripts in MATLAB and Python are available at our GitHub repository (https://github.com/Ez0-9606/NTS_Neuro-haemodynamic). Analyses were performed with custom scripts in MATLAB (R2022b) and Python (3.9.7) with Anaconda3 virtual environment (conda 23.1.0). We used open-source Python packages: numpy (1.21.5), scipy (1.7.3), ripser (0.6.1), scikit-learn (1.0.2), matplotlib (3.5.1), klusta (3.0.16), and phy (2.0b1). Python was used to detect and sort neuronal spikes and compute the latent space and its effective dimension. MATLAB was used for other analyses, including the preprocessing of electrophysiological neural and hemodynamics recording, the disentanglement of neuro-hemodynamic axis, the normalization of neural trajectory, and the construction of H-BIND.

References

Takeuchi, Y. et al. Closed-loop stimulation of the medial septum terminates epileptic seizures. Brain 144, 885–908 (2021).
Article PubMed Google Scholar
Kang, W. et al. Closed-loop direct control of seizure focus in a rodent model of temporal lobe epilepsy via localized electric fields applied sequentially. Nat. Commun. 13, 7805 (2022).
Article CAS PubMed PubMed Central Google Scholar
Wagner, F. B. et al. Targeted neurotechnology restores walking in humans with spinal cord injury. Nature 563, 65–71 (2018).
Article CAS PubMed Google Scholar
Victor, R. G. Carotid baroreflex activation therapy for resistant hypertension. Nat. Rev. Cardiol. 12, 451–463 (2015).
Article PubMed Google Scholar
Heusser, K. et al. Carotid baroreceptor stimulation, sympathetic activity, baroreflex function, and blood pressure in hypertensive patients. Hypertension 55, 619–626 (2010).
Article CAS PubMed Google Scholar
Bisognano, J. D. et al. Baroreflex activation therapy lowers blood pressure in patients with resistant hypertension: results from the double-blind, randomized, placebo-controlled Rheos Pivotal Trial. J. Am. Coll. Cardiol. 58, 765–773 (2011).
Article PubMed Google Scholar
National Academies of Sciences, Engineering, and Medicine. Foundational Research Gaps and Future Directions for Digital Twins (National Academies Press, 2023).
Gevers, M. Identification for control: from the early achievements to the revival of experiment design. Eur. J. Control 11, 335–352 (2005).
Article Google Scholar
Mun, J., Lee, J., Park, E. & Park, S. M. Frequency-dependent depression of the NTS synapse affects the temporal response of the antihypertensive effect of auricular vagus nerve stimulation (aVNS). J. Neural Eng. 19, 046039 (2022).
Article Google Scholar
Fan, W., Reynolds, P. J. & Andresen, M. C. Baroreflex frequency-response characteristics to aortic depressor and carotid sinus nerve stimulation in rats. Am. J. Physiol. Heart Circ. Physiol. 271, H2218–H2227 (1996).
Article CAS Google Scholar
Wehrwein, E. A., Orer, H. S. & Barman, S. M. Overview of the anatomy, physiology, and pharmacology of the autonomic nervous system. Regulation 37, 125 (2016).
Google Scholar
Nilsson, S. Comparative anatomy of the autonomic nervous system. Auton. Neurosci. 165, 3–9 (2011).
Article PubMed Google Scholar
Saper, C. B. The central autonomic nervous system: conscious visceral perception and autonomic pattern generation. Annu. Rev. Neurosci. 25, 433–469 (2002).
Article CAS PubMed Google Scholar
Guyenet, P. G. The sympathetic control of blood pressure. Nat. Rev. Neurosci. 7, 335–346 (2006).
Article CAS PubMed Google Scholar
Amendt, K., Czachurski, J., Dembowsky, K. & Seller, H. Bulbospinal projections to the intermediolateral cell column; a neuroanatomical study. J. Auton. Nerv. Syst. 1, 103–117 (1979).
Article CAS PubMed Google Scholar
Liu, Z., Chen, C. Y. & Bonham, A. C. Frequency limits on aortic baroreceptor input to nucleus tractus solitarii. Am. J. Physiol. Heart Circ. Physiol. 278, H577–H585 (2000).
Article CAS PubMed Google Scholar
Kolpakova, J. et al. Responses of nucleus tractus solitarius (NTS) early and late neurons to blood pressure changes in anesthetized F344 rats. PLoS ONE 12, e0169529 (2017).
Article PubMed PubMed Central Google Scholar
Rogers, R. F., Paton, J. & Schwaber, J. S. NTS neuronal responses to arterial pressure and pressure changes in the rat. Am. J. Physiol. Regul. Integr. Comp. Physiol. 265, R1355–R1368 (1993).
Article CAS Google Scholar
Andresen, M. C., Doyle, M. W., Jin, Y. H. & Bailey, T. W. Cellular mechanisms of baroreceptor integration at the nucleus tractus solitarius. Ann. N. Y. Acad. Sci. 940, 132–141 (2001).
Article CAS PubMed Google Scholar
Cowley, A. W. Jr, Liard, J. F. & Guyton, A. C. Role of the baroreceptor reflex in daily control of arterial blood pressure and other variables in dogs. Circ. Res. 32, 564–576 (1973).
Article PubMed Google Scholar
Kirchheim, H. R. Systemic arterial baroreceptor reflexes. Physiol. Rev. 56, 100–177 (1976).
Article CAS PubMed Google Scholar
Andresen, M. C. & Kunze, D. L. Nucleus tractus solitarius—gateway to neural circulatory control. Annu. Rev. Physiol. 56, 93–116 (1994).
Article CAS PubMed Google Scholar
Ran, C., Boettcher, J. C., Kaye, J. A., Gallori, C. E. & Liberles, S. D. A brainstem map for visceral sensations. Nature 609, 320–326 (2022).
Article CAS PubMed PubMed Central Google Scholar
Sclocco, R. et al. Stimulus frequency modulates brainstem response to respiratory-gated transcutaneous auricular vagus nerve stimulation. Brain Stimul. 13, 970–978 (2020).
Article PubMed PubMed Central Google Scholar
Accorsi-Mendonça, D. & Machado, B. H. Synaptic transmission of baro- and chemoreceptors afferents in the NTS second order neurons. Auton. Neurosci. 175, 3–8 (2013).
Article PubMed Google Scholar
Spyer, K. Annual review prize lecture. Central nervous mechanisms contributing to cardiovascular control. J. Physiol. 474, 1 (1994).
Article CAS PubMed PubMed Central Google Scholar
Bailey, T. W., Appleyard, S. M., Jin, Y. H. & Andresen, M. C. Organization and properties of GABAergic neurons in solitary tract nucleus (NTS). J. Neurophysiol. 99, 1712–1722 (2008).
Article PubMed Google Scholar
Baptista, V., Zheng, Z., Coleman, F., Rogers, R. & Travagli, R. Characterization of neurons of the nucleus tractus solitarius pars centralis. Brain Res. 1052, 139–146 (2005).
Article CAS PubMed PubMed Central Google Scholar
Seagard, J., Dean, C. & Hopp, F. Discharge patterns of baroreceptor-modulated neurons in the nucleus tractus solitarius. Neurosci. Lett. 191, 13–18 (1995).
Article CAS PubMed Google Scholar
Kawai, Y. & Senba, E. Organization of excitatory and inhibitory local networks in the caudal nucleus of tractus solitarius of rats revealed in in vitro slice preparation. J. Comp. Neurol. 373, 309–321 (1996).
Article CAS PubMed Google Scholar
Zhang, J. & Mifflin, S. W. Differential roles for NMDA and non-NMDA receptor subtypes in baroreceptor afferent integration in the nucleus of the solitary tract of the rat. J. Physiol. 511, 733–745 (1998).
Article CAS PubMed PubMed Central Google Scholar
Pandarinath, C. et al. Latent factors and dynamics in motor cortex and their application to brain–machine interfaces. J. Neurosci. 38, 9390–9401 (2018).
Article CAS PubMed PubMed Central Google Scholar
Ju, H. & Bassett, D. S. Dynamic representations in networked neural systems. Nat. Neurosci. 23, 908–917 (2020).
Article CAS PubMed Google Scholar
Mastrogiuseppe, F. & Ostojic, S. Linking connectivity, dynamics, and computations in low-rank recurrent neural networks. Neuron 99, 609–623.e629 (2018).
Article CAS PubMed Google Scholar
Yang, G. R. & Wang, X.-J. Artificial neural networks for neuroscientists: a primer. Neuron 107, 1048–1070 (2020).
Article CAS PubMed PubMed Central Google Scholar
Sussillo, D., Churchland, M. M., Kaufman, M. T. & Shenoy, K. V. A neural network that finds a naturalistic solution for the production of muscle activity. Nat. Neurosci. 18, 1025–1033 (2015).
Article CAS PubMed PubMed Central Google Scholar
Pandarinath, C. et al. Inferring single-trial neural population dynamics using sequential auto-encoders. Nat. Methods 15, 805–815 (2018).
Article CAS PubMed PubMed Central Google Scholar
Churchland, M. M. et al. Neural population dynamics during reaching. Nature 487, 51–56 (2012).
Article CAS PubMed PubMed Central Google Scholar
Shenoy, K. V. & Kao, J. C. Measurement, manipulation and modeling of brain-wide neural population dynamics. Nat. Commun. 12, 633 (2021).
Article CAS PubMed PubMed Central Google Scholar
Kriegeskorte, N. & Wei, X.-X. Neural tuning and representational geometry. Nat. Rev. Neurosci. 22, 703–718 (2021).
Article CAS PubMed Google Scholar
Gallego, J. A., Perich, M. G., Miller, L. E. & Solla, S. A. Neural manifolds for the control of movement. Neuron 94, 978–984 (2017).
Article CAS PubMed PubMed Central Google Scholar
Mante, V., Sussillo, D., Shenoy, K. V. & Newsome, W. T. Context-dependent computation by recurrent dynamics in prefrontal cortex. Nature 503, 78–84 (2013).
Article CAS PubMed PubMed Central Google Scholar
Vyas, S., Golub, M. D., Sussillo, D. & Shenoy, K. V. Computation through neural population dynamics. Annu. Rev. Neurosci. 43, 249–275 (2020).
Article CAS PubMed PubMed Central Google Scholar
Chaudhuri, R., Gerçek, B., Pandey, B., Peyrache, A. & Fiete, I. The intrinsic attractor manifold and population dynamics of a canonical cognitive circuit across waking and sleep. Nat. Neurosci. 22, 1512–1520 (2019).
Article CAS PubMed Google Scholar
Gardner, R. J. et al. Toroidal topology of population activity in grid cells. Nature 602, 123–128 (2022).
Article CAS PubMed PubMed Central Google Scholar
Sauerbrei, B. A. et al. Cortical pattern generation during dexterous movement is input-driven. Nature 577, 386–391 (2020).
Article CAS PubMed Google Scholar
Dabagia, M., Kording, K. P. & Dyer, E. L. Aligning latent representations of neural activity. Nat. Biomed. Eng. 7, 337–343 (2023).
Article PubMed Google Scholar
Grüne, L. & Pannek, J. Nonlinear Model Predictive Control: Theory and Algorithms (Springer London, 2011).
Lewis, F. L., Vrabie, D. & Syrmos, V. L. Optimal Control (John Wiley & Sons, 2012).
Yang, Y. et al. Modelling and prediction of the dynamic responses of large-scale brain networks during direct electrical stimulation. Nat. Biomed. Eng. 5, 324–345 (2021).
Article PubMed Google Scholar
Khambhati, A. N. et al. Functional control of electrophysiological network architecture using direct neurostimulation in humans. Netw. Neurosci. 3, 848–877 (2019).
Article PubMed PubMed Central Google Scholar
Zhang, Q. et al. A prototype closed-loop brain–machine interface for the study and treatment of pain. Nat. Biomed. Eng. 7, 533–545 (2023).
Article PubMed Google Scholar
Houghton, L. A., Lee, A. S., Badri, H., DeVault, K. R. & Smith, J. A. Respiratory disease and the oesophagus: reflux, reflexes and microaspiration. Nat. Rev. Gastroenterol. Hepatol. 13, 445–460 (2016).
Article PubMed Google Scholar
Travagli, R. A., Hermann, G. E., Browning, K. N. & Rogers, R. C. Brainstem circuits regulating gastric function. Annu. Rev. Physiol. 68, 279–305 (2006).
Article CAS PubMed PubMed Central Google Scholar
Tohyama, T. et al. Smart baroreceptor activation therapy strikingly attenuates blood pressure variability in hypertensive rats with impaired baroreceptor. Hypertension 75, 885–892 (2020).
Article CAS PubMed Google Scholar
Smith, J., Linderman, S. & Sussillo, D. Reverse engineering recurrent neural networks with Jacobian switching linear dynamical systems. Adv. Neural Inf. Process. Syst. 34, 16700–16713 (2021).
Google Scholar
Spyer, K. M. & Gourine, A. V. Chemosensory pathways in the brainstem controlling cardiorespiratory activity. Philos. Trans. R. Soc. B Biol. Sci. 364, 2603–2610 (2009).
Article Google Scholar
Sani, O. G., Abbaspourazad, H., Wong, Y. T., Pesaran, B. & Shanechi, M. M. Modeling behaviorally relevant neural dynamics enabled by preferential subspace identification. Nat. Neurosci. 24, 140–149 (2021).
Article CAS PubMed Google Scholar
Wang, W.-Z., Gao, L., Pan, Y.-X., Zucker, I. H. & Wang, W. AT1 receptors in the nucleus tractus solitarii mediate the interaction between the baroreflex and the cardiac sympathetic afferent reflex in anesthetized rats. Am. J. Physiol. Regul. Integr. Comp. Physiol. 292, R1137–R1145 (2007).
Article CAS PubMed Google Scholar
Mifflin, S. W. & Felder, R. B. Synaptic mechanisms regulating cardiovascular afferent inputs to solitary tract nucleus. Am. J. Physiol. Heart Circ. Physiol. 259, H653–H661 (1990).
Article CAS Google Scholar
Chen, C.-Y., Horowitz, J. M. & Bonham, A. C. A presynaptic mechanism contributes to depression of autonomic signal transmission in NTS. Am. J. Physiol. Heart Circ. Physiol. 277, H1350–H1360 (1999).
Article CAS Google Scholar
McDougall, S. J. & Andresen, M. C. Low‐fidelity GABA transmission within a dense excitatory network of the solitary tract nucleus. J. Physiol. 590, 5677–5689 (2012).
Article CAS PubMed PubMed Central Google Scholar
Mun, J., Lee, J. & Park, S.-M. Real-time closed-loop brainstem stimulation modality for enhancing temporal blood pressure reduction. Brain Stimul. 17, 826–835 (2024).
Article PubMed Google Scholar
Paxinos, G. & Watson, C. The Rat Brain in Stereotaxic Coordinates: Hard Cover Edition (Elsevier, 2006).
Thireau, J., Zhang, B., Poisson, D. & Babuty, D. Heart rate variability in mice: a theoretical and practical guide. Exp. Physiol. 93, 83–94 (2008).
Article CAS PubMed Google Scholar
Aubert, A. E. et al. The analysis of heart rate variability in unrestrained rats. Validation of method and results. Comput. Methods Prog. Biomed. 60, 197–213 (1999).
Article CAS Google Scholar
Kuwahara, M. et al. Power spectral analysis of heart rate variability as a new method for assessing autonomic activity in the rat. J. Electrocardiol. 27, 333–337 (1994).
Article CAS PubMed Google Scholar
Buzsáki, G. Large-scale recording of neuronal ensembles. Nat. Neurosci. 7, 446–451 (2004).
Article PubMed Google Scholar
Harris, K. D., Henze, D. A., Csicsvari, J., Hirase, H. & Buzsaki, G. Accuracy of tetrode spike separation as determined by simultaneous intracellular and extracellular measurements. J. Neurophysiol. 84, 401–414 (2000).
Article CAS PubMed Google Scholar
Rossant, C. et al. Spike sorting for large, dense electrode arrays. Nat. Neurosci. 19, 634–641 (2016).
Article CAS PubMed PubMed Central Google Scholar
Jun, J. J. et al. Fully integrated silicon probes for high-density recording of neural activity. Nature 551, 232–236 (2017).
Article CAS PubMed PubMed Central Google Scholar
Williamson, R. C. et al. Scaling properties of dimensionality reduction for neural populations and network models. PLoS Comput. Biol. 12, e1005141 (2016).
Article PubMed PubMed Central Google Scholar
Bellman, R. Dynamic Programming (Princeton Univ. Press, Princeton, NJ, 1957).
Tenenbaum, J. B., Silva, V. D. & Langford, J. C. A global geometric framework for nonlinear dimensionality reduction. Science 290, 2319–2323 (2000).
Article CAS PubMed Google Scholar
Singh, G. et al. Topological analysis of population activity in visual cortex. J. Vis. 8, 11–11 (2008).
Article Google Scholar
Tralie, C., Saul, N. & Bar-On, R. Ripser.py: a lean persistent homology library for Python. J. Open Source Softw. 3, 925 (2018).
Article Google Scholar
Bauer, U. Ripser: efficient computation of Vietoris–Rips persistence barcodes. J. Appl. Comput. Topol. 5, 391–423 (2021).
Article Google Scholar
Carlsson, G. Topology and data. Bull. Am. Math. Soc. 46, 255–308 (2009).
Article Google Scholar
Ghrist, R. Barcodes: the persistent topology of data. Bull. Am. Math. Soc. 45, 61–75 (2007).
Article Google Scholar
Rubin, A. et al. Revealing neural correlates of behavior without behavioral measurements. Nat. Commun. 10, 4745 (2019).
Article PubMed PubMed Central Google Scholar
Nieh, E. H. et al. Geometry of abstract learned knowledge in the hippocampus. Nature 595, 80–84 (2021).
Article CAS PubMed PubMed Central Google Scholar
Markus, E. J., Barnes, C. A., McNaughton, B. L., Gladden, V. L. & Skaggs, W. E. Spatial information content and reliability of hippocampal CA1 neurons: effects of visual input. Hippocampus 4, 410–421 (1994).
Article CAS PubMed Google Scholar
Dubreuil, A., Valente, A., Beiran, M., Mastrogiuseppe, F. & Ostojic, S. The role of population structure in computations through neural dynamics. Nat. Neurosci. 25, 783–794 (2022).
Article CAS PubMed PubMed Central Google Scholar
Wilson, H. R. & Cowan, J. D. Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1–24 (1972).
Article CAS PubMed PubMed Central Google Scholar
Chow, C. C. & Karimipanah, Y. Before and beyond the Wilson–Cowan equations. J. Neurophysiol. 123, 1645–1656 (2020).
Article PubMed PubMed Central Google Scholar
Darshan, R. & Rivkind, A. Learning to represent continuous variables in heterogeneous neural networks. Cell Rep. 39, 110612 (2022).
Article CAS PubMed Google Scholar
Russo, A. A. et al. Neural trajectories in the supplementary motor area and motor cortex exhibit distinct geometries, compatible with different classes of computation. Neuron 107, 745–758.e746 (2020).
Article CAS PubMed PubMed Central Google Scholar
Sussillo, D. & Abbott, L. F. Generating coherent patterns of activity from chaotic neural networks. Neuron 63, 544–557 (2009).
Article CAS PubMed PubMed Central Google Scholar

Download references

Acknowledgements

This research was supported by the Pioneer Research Center Program through the National Research Foundation of Korea, funded by the Ministry of Science, ICT & Future Planning (2022M3C1A3081294); the National Research Foundation of Korea(NRF) grant funded by the Korea government(MSIT) (No. 2020R1A2C2005385).

Author information

Authors and Affiliations

Department of Convergence IT Engineering, Pohang University of Science and Technology (POSTECH), Pohang, Republic of Korea
Jiho Lee, Junseung Mun, Minhye Choo & Sung-Min Park
Medical Device Innovation Center, Pohang University of Science and Technology (POSTECH), Pohang, Republic of Korea
Jiho Lee, Junseung Mun, Minhye Choo & Sung-Min Park
Department of Electrical Engineering, Pohang University of Science and Technology (POSTECH), Pohang, Republic of Korea
Sung-Min Park
Department of Mechanical Engineering, Pohang University of Science and Technology (POSTECH), Pohang, Republic of Korea
Sung-Min Park
Institute of Convergence Science, Yonsei University, Seoul, Republic of Korea
Sung-Min Park

Authors

Jiho Lee
View author publications
You can also search for this author inPubMed Google Scholar
Junseung Mun
View author publications
You can also search for this author inPubMed Google Scholar
Minhye Choo
View author publications
You can also search for this author inPubMed Google Scholar
Sung-Min Park
View author publications
You can also search for this author inPubMed Google Scholar

Contributions

Conceptualization, J.L.; Methodology, J.L., J.M., and M.C.; Software, J.L.; Validation: J.L.; Formal Analysis: J.L. and M.C.; Investigation, J.L., J.M.; Resources: J.L.; Data Curation: J.L.; Writing-Original Draft, J.L.; Writing-Review & Editing, J.L. and S.M.P.; Funding Acquisition, J.L., J.M., and S.M.P.; Project Administration: J.L.; Supervision, S.M.P. All authors contributed to the final manuscript and approved its final version.

Corresponding author

Correspondence to Sung-Min Park.

Ethics declarations

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Cite this article

Lee, J., Mun, J., Choo, M. et al. Predictive modeling of hemodynamics during viscerosensory neurostimulation via neural computation mechanism in the brainstem. npj Digit. Med. 8, 220 (2025). https://doi.org/10.1038/s41746-025-01635-w

Download citation

Received: 24 March 2024
Accepted: 11 April 2025
Published: 23 April 2025
DOI: https://doi.org/10.1038/s41746-025-01635-w

Subjects

Abstract

Similar content being viewed by others

The brain–heart axis: integrative cooperation of neural, mechanical and biochemical pathways

Vagal sensory neurons mediate the Bezold–Jarisch reflex and induce syncope

Digital twin brain simulator for real-time consciousness monitoring and virtual intervention using primate electrocorticogram data

Introduction

Results

Stimulus-driven collective dynamics of the NTS neuronal population in latent space

Neuro-hemodynamic coupling via collective dynamics in latent space

Inter-individual consistency of collective dynamics in NTS latent space

Predictive modeling of stimulus-driven hemodynamics based on neural computation mechanism

Discussion

Methods

Animals

Surgery

Solitary tract stimulation

Electrophysiological neural recording

Hemodynamic recordings

Stimulus artifact rejection

Identification of single neurons

Classification of NTS neurons

Analysis of single-neuronal interactions

Dimensionality reduction for collective dynamics

Persistent homology

Effective dimension of latent space

Single-neuronal dependence on internal variables

Disentanglement of neuro-hemodynamic axis

Normalization of neural trajectory

Modeling the local neural circuit in NTS

H-BIND

Nonlinear model-predictive controller

Data analysis and statistics

Data availability

Code availability

References

Acknowledgements

Author information

Authors and Affiliations

Contributions

Corresponding author

Ethics declarations

Competing interests

Additional information

Supplementary information

Supplementary information

Rights and permissions

About this article

Cite this article

Share this article

Search

Quick links