Data-driven discovery of governing equations for coarse-grained heterogeneous network dynamics

Data-driven modeling paradigms are an increasingly important tool for the characterization of behavior in complex, high-dimensional systems. For instance, Reduced order models (ROMs) [1] aim to exploit the ubiquitous observation from experiment and/or computation that meaningful input/output signals in high-dimensional systems are encoded in low-dimensional patterns of dynamic activity. However, many complex systems of interest exhibit behavior across multiple temporal or spatial scales, which poses unique challenges for modeling and predicting their behavior. It is often the case that while we are primarily interested in macroscale phenomena, the microscale dynamics must also be modeled and understood, as they play a central role in driving the macroscale behavior. This can make dealing with multiscale systems particularly difficult unless the time scales are disambiguated in a principled way. There is a significant body of research focused on modeling multiscale systems to produce coarse-grained descriptions: notably the equation-free (EF) method[2] for linking scales and the heterogeneous multiscale modeling (HMM) framework [3]. Additional work has focused on testing for the presence of multiscale dynamics so that analyzing and simulating multiscale systems is more computationally efficient [4, 5]. Many of the same issues that make modeling multiscale systems difficult can also present challenges for model discovery and system identification. These challenges motivate the development of specialized methods for performing model discovery on problems with multiple time scales, taking into account the unique properties of multiscale systems. Specifically, we integrate dimensionality-reduction, sparse regression, and robust statistics to discover governing evolution equations for coarse-grained heterogeneous network dynamics as well as to explicitly disambiguate and model multiscale temporal phenomena (See Fig. 3).

The discovery and pairing of coordinates and dynamics is the hallmark feature of scientific modeling. For instance, the second century Ptolemaic description of celestial mechanics featured an earth-centric coordinate systems with the retrograde planetary dynamics expressed as linear superposition of circular motion of different frequencies, i.e. the doctrine of the perfect circle. By Newton’s time in the seventeenth century, the heliocentric coordinate system was established allowing for an ${\bf F}=m{\bf a}$ description of planetary motion with the inverse square law of gravitation. Thus it is not surprising that methods for jointly finding coordinates and dynamics have emerged as a central theme in modern data-driven science. Indeed, there is significant diversity in mathematical approaches for pairing a coordinate representation with an accompanying dynamical model that characterizes the evolution. Reduced order models (ROMs) accomplish this pairing by constructing coordinates from the most dominant correlated activity, or proper orthogonal decomposition (POD) [6, 1], and then Galerkin projecting on the governing equations [1]. More broadly, modal decomposition techniques, such as POD and dynamic mode decomposition (DMD) [7], approximate linear subspaces using dominant correlations in spatio-temporal data [8]. Linear subspaces, however, are highly restrictive and ill-suited to handle parametric dependencies. Attempts to circumvent these shortcomings include using multiple linear subspaces covering different temporal or spatial domains [9, 10, 11, 12], multi-resolution DMD [13], diffusion map embeddings  [14, 15, 16, 17], or more recently, using deep learning to compute underlying nonlinear subspaces which are advantageous for dynamics, both linear and nonlinear [18, 19, 20, 21, 22]. These techniques represent data-driven architectures for extracting order-parameter descriptions of the underlying spatio-temporal dynamics observed [23].

Such flexibility and diversity in algorithms for the joint discovery of coordinates and dynamics has allowed for many new representations of physical systems, including in the networked dynamical systems of interest here. Specifically, we are interested in modeling the emergent coarse-grained heterogenous behavior that arises in high-dimensional, networked dynamical systems whose underlying dynamics are strongly nonlinear. In particular, we consider networks of coupled nonlinear oscillators whose collective behavior approaches relaxation-oscillation limit cycle dynamics. Recently, oscillatory phenomena have been coarse-grained into learned and/or collective coordinates (intrinsic coordinates) for modeling their evolution [24, 25, 26, 27, 28, 29]. Gottwald originally introduced collective coordinates to describe coupled phase oscillators of the Kuramoto model [25], and then extended the mathematical framework to reduce infinite-dimensional stochastic partial differential equations (SPDEs) with symmetry to a set of finite-dimensional stochastic differential equations which describe the shape of the solution and the dynamics along the symmetry group [26]. Shlizerman et al [27] coarse-grained using coordinates associated with the dominant correlated activity, or singular value decomposition, and projected into a framework of neural activity measures. Kemeth et al [24] instead learn collective dynamics on a slow manifold, after initial transients have died out, which can be approximated through a learned model based on local spatial partial derivatives in the emergent coordinates. These works are closely related to the present work, as their aims are well-aligned with the goals of our method. However, there are key differences. Specifically, we consider coarse-graining on a system for which regions of fast activity are critical for characterizing the dynamics, i.e. the dynamics produce temporal boundary layers. In particular, we wish to learn the dominant balance dynamics that occur in both fast and slow temporal regimes on network level dynamics.

In this work, we leverage dimensionality-reduction, sparse regression, and robust statistics to discover coarse-grained models of heterogeneous networked dynamical systems. We consider networks consisting of several types of oscillators and explore how changing the network structure influences the collective dynamics. Nishikawa et al [29] and Motter et al [30] consider synchronization of networks, including scale free and small world, as a function of the homogeneous or heterogeneous distribution of connectivity, showing that homogeneous networks are more easily synchronized. Here, we focus on aspects of the emergent dynamics that arise from the heterogeneity of the oscillators themselves. Specifically, in many regimes, the emergent dynamics exhibit two key features we aim to model: (i) nontrivial coarse-grained dynamics (relaxation-oscillations), and (ii) temporal boundary layers (rapid transient dynamics). Specifically, the proposed mathematical architecture discovers parsimonious governing evolution equations for coarse-grained network dynamics and, if desired, explicitly disambiguates emergent temporal boundary layer phenomena. The developed algorithm starts simply by projecting onto coordinates associated with the dominant correlated activity, i.e. by taking the singular value decomposition. We then identify rapid transition regions associated with transients in the time dynamics of the dominant modes by applying the sparse identification of nonlinear dynamics (SINDy)) [18] algorithm with data-trimming [31], which is a standard analysis tool from robust statistics. The governing equations produced by this stage of the method can be used to describe the full limit cycle. However, if a more detailed model is desired, the results of trimming are used to segment the coarse-grained trajectory, separating regions of rapid transition from the slow-temporal field, which is equivalent to the slow center-manifold considered previously [24, 25]. We then utilize the SINDy framework again to discover parsimonious governing equations for each region of the coarse-grained trajectory, constructing a hybrid model for the overall dynamics. Here we apply the full method to a network of coupled Fitzhugh-Nagumo oscillators, which produce rapid (spiking) behavior. The strongly nonlinear nature of the underlying agents leads to the existence of temporal boundary layers in the collective dynamics. The overall architecture allows for a data-driven approach to characterizing dominant-balance physics [23, 32] and singularly perturbed problems [33, 34, 35], which is a classic technique for understanding, for instance, boundary layer formation in fluid dynamics.

Coarse-graining is of broad scientific interest across a diversity of disciplines. Coarse-graining seeks to determine a set of variables that successfully summarize the detailed state of a network of interacting agents. Analytical and computational approaches have long been sought for achieving this end, with the former leveraging the dynamical systems theory of center manifold reductions, normal forms, bifurcation theory and perturbation analysis. Indeed, the seminal review of Cross and Hohenberg [23] detail the many mathematical architectures available for deriving order-parameter descriptions, or coarse-grained models, that characterize the dynamics. More recent computational approaches include statistical and machine learning algorithms [36, 37] for network clustering, principal component analysis, and, for instance, algorithms for detecting density (k-core decomposition). The diversity of methods developed allow for a new coordinate representation for the collective dynamics. This representation can then be leveraged by model discovery methods such as SINDy to determine the dynamics on the coarse-grained coordinates.

SINDy recovers parsimonious representations of the dynamics from measurement data by sparse regression to a library of candidate models [18]. Thus the goal of SINDy is to promote a parsimonious representation of the complex system by traditional and interpretable governing equations. Consider the nonlinear dynamical system, which will ultimately be our coarse-grained representation,

$$\frac{d{\bf x}}{dt}=f({\bf x},t)$$

(1)

which is measured at time points ${\bf t}=[t_{1},t_{2},\cdots,t_{m}]$. From measurements, we construct the matrix ${\bf X}({\bf t})$ = $[{{\bf x}_{1}({\bf t})\,\,{\bf x}_{2}({\bf t})\,\,\dots\,\,{\bf x}_{n}({\bf t})}]$ $\in\mathbb{R}^{m\times n}$. The method introduced in  [18] seeks to identify ${\bf f}$ via sequential threshold least-squares, which is a proxy for the sparsifying zero-norm. The set of $n$ state measurements are used to populate a library of candidate nonlinear terms ${\bf\Theta(X)=[1^{\top}\;X^{\top}\;(X\otimes X)^{\top}\;\cdots\;\text{sin}(X)^{\top}]}$, where ${\bf x\otimes y}$ defines the vector of all product combinations of the state components. Each candidate term should be unique, as a suitable library is crucial in the SINDy algorithm. A common strategy is start with polynomials and increase the complexity of the library with other terms, such as trigonometric functions. Thus, the system in (1) is approximated by:

$${\bf\dot{X}=\Theta(X)\Xi}.$$

(2)

The time derivatives ${\bf\dot{X}}({\bf t})=[{\bf\dot{x}}_{1}({\bf t})\,\,{\bf\dot{x}}_{2}({\bf t})\,\,\dots\,\,{\bf\dot{x}}_{n}({\bf t})]$, if not measured directly, can be found via numerical differentiation and should be appropriately de-noised, if necessary (low pass Butterworth, total variation regularization, etc.) [38]. The coefficients ${\bf\Xi}$ are the sparse weightings of the corresponding candidate library terms. Therefore, our regression relies on sparse regularization to enforce a parsimonious $\bf\Xi$ corresponding to the fewest nonlinear terms in our library that describe our dynamics well:

$${\bf\Xi}=\text{arg}\,\min\limits_{\hat{\bf\Xi}\,}\|{\bf\Theta(X)\hat{\bf\Xi}-\dot{X}}\|_{\text{2}}+\lambda\|{\bf\hat{\bf\Xi}}\|_{\text{0}}$$

(3)

Regressing to the zero-norm is often achieved by relaxing the the one-norm. However, modern optimization frameworks are allowing for computationally tractable proxies for the zero-norm that are superior to the one-norm relaxation [31, 39]. SINDy is modular and adaptible, allowing for modifications which include the discovery of spatio-temporal systems [40, 41], multiscale dynamics [20], parametric dependencies [42], hybrid dynamical systems [43], systems subject to control [44, 45], and implicit dynamics [46, 47, 48]. Further, it can be used in the low-data limit [49, 50, 51] and with denoising architectures [52, 53]. The SINDy algorithm, which has an accompanying python package [54], is the workhorse algorithm for our coarse-graining method.

To disambiguate fast scale dynamics, our coarse-graining method exploits a key extension of SINDy: data trimming. Modifying Eq. (3) to include trimming yields:

$$\begin{split}{\bf\Xi},{\bf v}~{}=~{}&\text{arg}\,\min\limits_{\hat{\bf\Xi},{\bf v}}\sum_{i=1}^{m}\frac{1}{2}v_{i}\|{(\bf\Theta(X)\hat{\bf\Xi}-\dot{X})_{i}}\|_{\text{2}}+\lambda\|{\bf\hat{\bf\Xi}}\|_{\text{0}}\\ &\text{s.t.}~{}~{}~{}~{}0\leq v_{i}\leq 1,~{}~{}~{}~{}~{}{\bf 1}^{T}\bf{v}=h,\end{split}$$

(4)

where $h<m$ is an estimate of the number of the $m$ input data points that are considered “inliers.” See the work of Champion et al. for a thorough explanation [31]. After the SINDy with trimming algorithm is applied, the vector ${\bf v}$ has lower entries where data points are harder to fit to a parsimonious model of system dynamics. In systems that exhibit multi-scale temporal dynamics, these troublesome points often correspond to spikes in the derivative. Thus, trimming effectively identifies regions of rapid change.

As an example, in Fig. 1a we illustrate how data-trimming partitions the limit cycle of a Rayleigh oscillator: the points which deviate from the slow-scale outer solution are trimmed, coinciding with the two boundary layers in the limit cycle. Trimming applied to Rossler dynamics on a strange attractor identifies rapid excursions in the z-component of the system (Fig. 1b). When applied to a low-rank representation of a network of spiking Fitzhugh-Nagumo (FHN) oscillators, trimming identifies regions of rapid spiking and relaxation (Fig. 1c). In each of these examples, the trimmed fraction, or $1-h/m$, is a hyperparameter of the method that must be tuned according to the proportion of the input measurements that appear to belong to regions of fast-scale dynamics. If the fraction is set too high, then points adjacent to the fast scale dynamics will also be trimmed. Conversely, if the trimmed fraction is set too low then only a subset of the fast-scale dynamics will be identified. Another factor to consider in applying this method is that in order to effectively repurpose an outlier-detection technique from robust statistics to partition our data, the data must be clean. If the input trajectory is noisy, then data trimming will pull out the noisy data points and their neighbors, whose derivative estimates are corrupted, and, possibly, also the desired regions of rapid scale dynamics.

The steps of our method for data-driven identification of temporal boundary layer phenomena in relaxation oscillations are outlined on the canonical Rayleigh oscillator in Fig. 2. First, SINDy with trimming is applied to the input trajectory to partition the data. Any reasonable library can be used at this point, as we only need the output from trimming, i.e. the vector ${\bf v}$ from (4), not the model coefficients. The entries of ${\bf v}$ are compared to a threshold value, and points with entries below the threshold are assigned to the fast-scale group, $X_{fast}$, while points with an entry above the threshold are assigned to the slow-scale group, $X_{slow}$. Next, the locations of the points belonging to $X_{fast}$ within the phase space are used to bound the region in which the fast scale dynamics will apply. The bounding box is formed around each stretch of consecutive trimmed points, and then overlapping boxes are joined resulting in one region for each segment of the limit cycle displaying fast-scale dynamics. Points outside of these bounding boxes are assigned to $X_{slow}$ during simulation. Then, SINDy is applied separately to data in $X_{fast}$ and $X_{slow}$. For the slow-scale group, the candidate library of functions is built using $x\in X_{slow}$. For the fast scale group, the candidate library of functions is built using $x\in X_{fast}$. With these two dynamic models in hand, and a partition on the phase space determining when each model applies, we can simulate the system dynamics. At each time step, before calculating an update to the system position, we check whether the point belongs to $X_{fast}$ or $X_{slow}$. If $x\in X_{fast}$, then the fast scale dynamics are used to calculate an update. If $x\in X_{slow}$, then the slow scale dynamics are used to calculate an update. Note that the hybrid model is only valid for points on the limit cycle of the system. In this example, trajectories initiated elsewhere in the phase space either remain on the line $y=0$, approach the line $y=1$ in the upper half plane, or approach the line $y=-1$ in the negative half plane.

We also demonstrate this hybrid model approach on the limit cycle of a Van der Pol oscillator exhibiting relaxation-oscillations, as shown in Fig. 3. Like the Rayleigh Oscillator example, we first partition the input trajectory using trimming. Note that the trimmed points, shown in red, are those points that deviate from the outer solution predicted by asymptotic approximation theory, shown by a dashed line. We then use the set of trimmed points to bound the regions of phase space governed by fast scale dynamics, and the remaining phase space is considered to be governed by slow scale dynamics. Next we apply SINDy to data points from the slow scale scale region, fitting the dynamics to a library that depends on $x$ and includes rational functions. In this case, we fit the dynamics of the two regions of fast scale dynamics separately, calling those in the lower half plane $X_{fast^{-}}$ and those in the upper half plane $X_{fast^{+}}$. For both subsets of data, we fit the dynamics to a library of cubic polynomials in $x$ and $y$. We treat the two segments of fast dynamics separately because fitting both regions at once results in an average between the two desired models that does not approximate either branch fully. We note, however, that the resulting models for $X_{fast^{-}}$ and $X_{fast^{+}}$ reported in Fig. 3, can easily be transformed to match each other by applying an absolute value sign to $x$. If such symmetries are clear from viewing the limit cycle of a system, then the library of candidate functions used in SINDy could be built using polynomials of $|x|$, and both regions of fast scale dynamics could be fit together. Again simulation of the system is accomplished by checking whether the current state belongs to $X_{slow},X_{fast^{+}},$ or $X_{fast^{-}}$, making an update with the appropriate dynamic rules, and then iterating these two steps until the desired time span is covered.

For our coarse-graining method, a reduced-order basis is first derived through the singular value decomposition [55, 56] (Fig. 4) and second, SINDy is deployed to discover the dynamics in this reduced order basis (Fig. 5). If the low-rank dynamics exhibit multiscale temporal behavior, then SINDy with trimming can be used instead to identify regions of fast-scale dynamics and fit a hybrid model to the data, facilitating analysis of the dominant balance physics in distinct dynamic intervals of the collective limit cycle. This method allows us to tackle coarse-graining for networks displaying heterogenous spatio-temporal dynamics that are not tractable with analytical approaches. For very high-dimensional problems, randomized algorithms [57] and compressed decompositions [58] can be exploited in order to compute the low-rank subspaces more efficiently.

The systems that we consider here are networks of coupled oscillators of the form

$$\frac{d{\bf x}}{dt}=f({\bf x},t),$$

(5)

where ${{\bf x}(t)}=(x_{1},x_{2},..,x_{n})$ is a vector of state variables, or nodes, and $f({\bf x})=(f_{1}({\bf x}),f_{2}({\bf x}),...,f_{n}({\bf x}))$ is a vector of functions encoding the dynamics on each node. The dynamics depend on intrinsic properties of each node as well as coupling to neighboring nodes. Neighbors are defined by a symmetric, unweighted adjacency matrix, $A\in\mathbb{R}^{n\times n}$, where the entries $A_{ij}=A_{ji}=1$ if node $i$ and node $j$ are connected and $A_{ij}=A_{ji}=0$ otherwise. Here we consider Erdos-Renyi adjacency matrices, so $A$ may be sparse or highly connected as we tune the edge probability, $p$. The functions $f_{i}({\bf x})$ could take on any arbitrary form, but we limit our attention to a selection of well-known oscillators.

Before applying our method to heterogeneous systems, we establish that it is producing reasonable results by considering a network consisting of Kuramoto oscillators with sinusoidal coupling. For the ubiquitous Kuramoto oscillator, the phase of node $j$, denoted $\theta_{j}$, is governed by

$$\dot{\theta}_{j}=\omega_{j}+\frac{K}{n}\sum_{i=1}^{n}A_{ij}\sin(\theta_{i}-\theta_{j}).$$

(6)

The key parameters besides the aforementioned adjacency matrix, are the connectivity strength $K>0$, and the intrinsic frequency $\omega_{j}$. For each node in the system, $\omega_{j}$ is drawn from a uniform distribution. Kuramoto oscillators have been widely studied and exhibit a well-known phase transition from asynchrony below a critical coupling threshold to partial synchrony above this threshold. The system continues to synchronize more fully as the coupling parameter is further increased. Panel (a) of Figure 4 illustrates the behavior of a single Kuramoto oscillator as well as the full system dynamics for a network of Kuramoto oscillators. Note that we are visualising the cosine of the phase. Taking the singular value decomposition and projecting onto the first few modes reduces the system dynamics to a stable elliptical limit cycle. Retaining two of these modes and applying the SINDy method to the 2D limit cycle produces a simple and accurate set of governing equations which can be used to simulate the system dynamics (Fig. 5a, final column). With this friendly, linear collective behavior no trimming is required to recover useful model equations.

Next, we consider a slightly more complicated system: Rayleigh oscillators with nonlinear coupling. In this case each node follows the second order governing equation

$$\epsilon_{j}\ddot{x}_{j}=\dot{x}_{j}-\frac{\dot{x}^{3}}{3}-x+\frac{K}{n}\sum_{i=1}^{n}A_{ij}(1+(x_{j}-x_{i})^{2})(\dot{x}_{i}-\dot{x}_{j})$$

(7)

where $\epsilon_{j}<<1$ is drawn from a uniform distribution and again $K>0$ is the connectivity parameter. The nonlinear coupling strategy, termed Haken-Kelso-Bunz (HKB) coupling, was introduced to study human bimanual experiments and applied to data capturing the synchronization of people rocking chairs by Alderisio et al. [59]. In numerical simulations, we observe partial synchrony in these systems when there is sufficient coupling (Fig. 4b). This collective behavior is captured by a nonlinear limit cycle in the low-dimensional representation. Applying SINDy to the 4D projection with trimming pulls out regions of rapid relaxation along the direction of $U_{4}$ and produces a model that is linear in $U_{1}$ and $U_{2}$ and cubic in $U_{3}$ and $U_{4}$. Here each oscillator has an intrinsic $\epsilon_{j}$ on the order of $10^{-3}$. The multi-scale nature of the collective system dynamics is evident in the fact that the coefficients learned for modes $U_{2}$ and $U_{4}$ are much higher than those learned for $U_{1}$ and $U_{3}$. With a distribution of $\epsilon$ values centered closer to $1$, we would see a smaller gap between scales as $1/\epsilon\rightarrow 1$. Comparing the trajectory predicted by the model against the true low-rank trajectory of the system (final column of Fig. 5), we see that after one lap of the limit cycle, the prediction gets out of phase with the true trajectory. The SINDy model is attracted to a smaller stable limit cycle parallel to the true dynamics. Applying this method to systems with larger values of $\epsilon$, i.e. systems for which the boundary layers are less sharp, results in model dynamics that stay on track for longer time.

We also demonstrate our method on a network of Rossler oscillators with sinusoidal coupling. In this case each node has 3 state variables governed by

$$\begin{split}\dot{x}_{j}&=-y_{j}-z_{j}+\frac{K}{n}\sum_{i=1}^{n}A_{ij}\sin(x_{j}-x_{i})\\ \dot{y}_{j}&=x_{j}+ay_{j}+\frac{K}{n}\sum_{i=1}^{n}A_{ij}\sin(y_{j}-y_{i})\\ \dot{z}_{j}&=b+z_{j}(x_{j}-c)+\frac{K}{n}\sum_{i=1}^{n}A_{ij}\sin(z_{j}-z_{i})\end{split}$$

(8)

The parameters $(a,b,c)$ are uniform across all nodes in the network for a given run. In the example shown in Fig. 4c, $(a=0.2,b=0.2,c=5.7)$. The trajectory of each node is drawn to a similarly shaped strange attractor, however the attractors may be translated in phase space depending on the initial condition of the node. The nodes of the network largely synchronize phases and frequencies in their $x$ component, synchronize frequencies of their $y$ component, and partially synchronize the timing of excursions in the $z$ component. Projecting onto the first three singular vectors, the system lands on an attractor that resembles a Rossler attractor with a few differences. Specifically, the attractor is tilted slightly with respect to the axes, so there is a background oscillation rate in $U_{3}$. Due to this tilt, the $U_{3}$ “excursions” also disrupt the background activity along $U_{1}$ and $U_{2}$ slightly. Applying SINDy with trimming to this low-dimensional trajectory results in a linear model for $U_{1}$ and $U_{2}$ and a dense cubic equation for the dynamics of $U_{3}$. Because the attractor for the collective dynamics is tilted in the SVD-space compared to an attractor for a single Rossler oscillator, these cross terms are necessary to capture the background frequency and the excursions. Observe that trimmed points are mostly spikes along the $U_{3}$ direction. In this case, the model predictions and the true trajectory match up very well for time steps that were used in the SVD projection initially. However, about time $80$ in Fig. 5c you can start to see that the model prediction is deviating from the true trajectory–capturing less of the $U_{3}$ excursions. Continuing to predict beyond the training data, the model still produces dynamics that resemble a collective Rossler oscillator but they lose fidelity to the true system dynamics.

The last type of oscillator that we considered were linearly coupled Fitzhugh-Nagumo oscillators. This simplified model of a neuron has two state variables: a voltage variable which displays rapid spiking and relaxation under an external stimulus and a slower recovery variable. Each FHN node is governed by the pair of equations

$$\begin{split}\dot{v}_{j}&=\alpha_{3}v_{j}^{3}+\alpha_{2}v_{j}^{2}+\alpha_{1}v_{j}-w_{j}+\frac{K}{n}\sum_{i=1}^{n}A_{ij}(v_{j}-v_{i})\\ \dot{w}_{j}&=cv_{j}-bw_{j}+\frac{K}{n}\sum_{i=1}^{n}A_{ij}(w_{j}-w_{i}).\end{split}$$

(9)

For these simulations we use parameters $\alpha=(-0.1,1.1,-1)$, $c=0.1$, and $b=0.1$ as studied in [27]. This combination yields semicoherent spiking behavior across the network (Fig. 4d). In the raster plot of the network, we see that subsets of the nodes are firing together. The limit-cycle that emerges from projecting this onto 3 SVD modes has a distinctive triangular shape. Each edge appears to result from the spiking of a different subset of nodes. Note that for different parameter values, the system may cease spiking and settle to a steady state or spiking may be more chaotic. For this system we projected to low dimensions by taking the SVD of the voltage and the recovery variables within the nodes separately. $U_{1},U_{2}$ and $U_{3}$ are the first three modes in the voltage variable. $U_{4}$, $U_{5}$, and $U_{6}$ are the first three modes in the recovery variable. Applying SINDy with trimming to this highly nonlinear, 6-dimensional limit cycle yields a model that is linear in the recovery variables and cubic in the voltage variables. The model predictions and the true dynamics of the system match up very well even over long time horizons (Fig. 5d).

The multiscale nature of the limit-cycle in the low-rank representation of our FHN network makes it a prime candidate for building a hybrid model. After the network dynamics are reduced to 6 dimensions, we can apply the same steps as outlined for the Rayleigh and Van der Pol oscillators. First, we apply SINDy with trimming. Trimming identifies regions of fast-scale dynamics along directions $U_{1}$ through $U_{4}$, which correspond to the three subsets of partially synchronized oscillators firing. In the low-rank dynamics these are the three corners of the triangular limit-cycle on the recovery modes and the three edges of the limit-cycle on the voltage modes (Fig. 6a). Along the limit cycle, alternating segments of fast and slow scale dynamics are used to partition phase-space into 6 regions, each of which has unique dominant-balance dynamics. So rather than just $X_{fast}$ and $X_{slow}$, we have $X^{1}_{fast}$, $X^{2}_{fast}$, $X^{3}_{fast}$, $X^{1}_{slow}$, $X^{2}_{slow}$, and $X^{3}_{slow}$. These six intervals are color coded in Fig. (6b). The derivatives of the low-rank state variables in each of the six regions are fit using a linear library. We can see that the slow dynamic regions (blue, orange and yellow) show great agreement in all variables. The fast dynamic regions are also well matched in the recovery modes, 4,5 and 6. The fit is not as clean at the ends of each fast region in the voltage modes, 1, 2 and 3. The coefficients for the 6 dynamic models can be used to assess the system behavior within each region (6c).

Modeling heterogeneous, multiscale physics remains a mathematically challenging proposition. Indeed, traditional computational techniques quickly become intractable when attempting to resolve such systems in both space and time. Emerging data-driven methods are enabling new mathematical architectures that can automate the discovery of coarse-grained coordinates and dynamics that are low-dimensional and parsimonious. This renders approximate models that allow for rapid evaluation of the heterogeneous, multiscale physics. More than that, these techniques can also be constructed to produce interpretable physics models.

We have proposed a set of algorithms for modeling heterogeneous, high-dimensional networked dynamical systems. Specifically, we have leveraged dimensionality-reduction, sparse regression, and robust statistics to discover interpretable, coarse-grained models of the underlying physics. The method further allows us to disambiguate between the resulting fast and slow scale dynamics of the system. Used in combination, these methods are able to: (i) identify low-dimensional embeddings (coordinates) on which to construct models, (ii) identify different time scales using the statistical robustification technique of data trimming, and (iii) identify interpretable parsimonious models of the coarse-grained dynamics at different timescales.

The mathematical architecture has been demonstrated with a series of numerical experiments on networked nonlinear oscillators. The oscillators, whether a single type of oscillator or a heterogeneous set of oscillators, evolve and interact to produce low-rank dynamics that often include relaxation oscillations. The literature is largely devoid of the consideration of such systems due to their complexity. However, the mathematical architecture developed here is able to extract meaningful models even with such high-dimensional heterogeneous interactions.

The authors acknowledge funding support from the Air Force Office of Scientific Research (AFOSR FA9550-19-1-0386) and from the National Science Foundation AI Institute in Dynamic Systems grant number 2112085.