We developed a coupled social–climate network model that links opinion dynamics and the climate system, where opinions directly translate into actions and actions reflect underlying opinions, to examine how this feedback shapes collective behavior and global temperature trajectories. In contrast to previous social–climate models that discretized opinions, we assumed opinions on climate change form a continuum, and were thereby able to capture more nuanced interactions. The model shows that resistance to behavior change, elevated mitigation costs, and slow response to climate events can result in a global temperature anomaly in excess of 2[Formula: see text]C. However, this outcome could be avoided by lowering mitigation costs and increasing the rate of interactions between individuals with differing opinions (social learning). Our model is one of the first to demonstrate the emergence of opinion polarization in a human-environment system. We predict that polarization of opinions in a population can be extinguished, and the population will adopt mitigation practices, when the response to temperature change is sensitive, even at higher mitigation costs. It also indicates that even with polarized opinion, an average promitigative opinion in the population can reduce emissions. Finally, our model underscores how frequent and unexpected social or environmental changes, such as policy changes or extreme weather events, can slow climate change mitigation. This analysis helps identify the factors that support achieving international climate goals, such as leveraging peer influence and decreasing stubbornness in individuals, reducing mitigation costs, and encouraging climate-friendly lifestyles. Our model offers a valuable new framework for exploring the integration of social and natural sciences, particularly in the domain of human behavioral change.
Population dynamics in fields such as molecular biology, epidemiology, and ecology exhibit highly stochastic and nonlinear behavior. In gene regulatory systems in particular, oscillations and multistability are especially common. Despite this, none of the currently available stochastic models for population dynamics are both accurate and computationally efficient for long-term predictions. A prominent model in this field, the linear noise approximation (LNA), is computationally efficient for tasks such as simulation, sensitivity analysis, and parameter estimation; however, it is only accurate for linear systems and short-time predictions. Other models may achieve greater accuracy across a broader range of systems, but they sacrifice computational efficiency and analytical tractability. This paper demonstrates that, with specific modifications, the LNA can accurately capture nonlinear dynamics in population processes. We introduce a new framework based on center manifold theory, a classical concept from nonlinear dynamical systems. This approach enables the identification of simple, system-specific modifications to the LNA, tailored to classes of qualitatively similar nonlinear dynamical systems. With these modifications, the LNA can achieve accurate long-term simulations without compromising computational efficiency. We apply our methodology to classes of oscillatory and bistable systems and present multiple examples from molecular population dynamics that demonstrate accurate long-term simulations alongside significant improvements in computational efficiency.
This paper proposes a sparse regression strategy for discovery of ordinary differential equations from incomplete and noisy data. Inference is performed over both equation parameters and state variables using a statistically motivated likelihood function. Sparsity is enforced by a selection algorithm which iteratively removes terms and compares models using statistical information criteria. Large scale optimization is performed using a second-order variant of the Levenb erg-Marquardt method, where the gradient and Hessian are computed via automatic differentiation. The proposed method is illustrated and tested on several systems with varying levels of noisy and incomplete data. Comparisons are made to a state-of-the-art algorithm for system identification, demonstrating competitiveness of the proposed approach.
We introduce a tool, DiscretePhasePortrait.jl, written in Julia, that generates an augmented phase portrait for the analysis of two-dimensional discrete-time mappings. The augmented phase portrait is a generalization of that presented by Streipert and Wolkowicz [Math. Biosci., 355 (2023), 108924]. The purpose is to augment the phase portrait with additional information to allow for the determination of the global dynamics of planar discrete systems with complex behaviors. The first generalization is to consider isoclines for any given set of directions rather than simply those associated with the coordinate directions. We show that choosing directions aligned with eigenvectors of the Jacobian at a fixed point aids in making conclusions about global dynamics. The second generalization is to add a curve indicating where the determinant of the Jacobian is zero and the image of that curve. This helps identify the range of the mapping, effectively reducing the space one needs to consider when determining global dynamics. This tool also allows the user to generate phase portraits for multiple iterations of the map; thus, a phase portrait for a twice-iterated map can be used to eliminate the complicating oscillatory behavior of orbits near fixed points that have negative eigenvalues for their Jacobian. We illustrate the use of the tool throughout the paper, applying it to several different systems, showing how it can be used to establish the global stability of fixed points of interest. We also use the tool and some additional analysis to provide alternative proofs to four previously open problems about the local stability of some bifurcating fixed points for certain maps.
Abstract. The resultant is a well-known tool in number theory and algebraic geometry used to find simultaneous zeros of two polynomials. Here we show how to use the resultant for the analysis of polynomial or rational systems of ordinary differentials to find saddle-node, Hopf, Takens–Bogdanov (codimension two), and codimension three bifurcations. In the process, we find new features of classical dynamical systems including the Goodwin feedback oscillator, the Oregonator, and the Brusselator, among others. In addition, we prove a new theorem (Theorem 4A) relating the resultant and the calculation of Hopf bifurcation curves.
We analyze models of minimal gene regulatory networks (GRNs) introduced in the context of X chromosome inactivation (XCI), an epigenetic process occurring in placental female mammals. A core module in these GRNs is a 2D toggle-switch model resulting from a mutual inhibition. We first perform a bifurcation analysis of the toggle-switch model, where the main difficulty arises from the unknown coordinates of the fixed points. In the symmetric case, we prove the occurrence of a pitchfork bifurcation and compute explicitly the bifurcation lines. In the asymmetric case, we design a constructive numerical approach, providing one simultaneously with the bifurcation lines and fixed-point coordinates, and illustrate the occurrence of a saddle-node bifurcation using numerical continuation tools. Bistability in the toggle-switch accounts for the possible choice between an activated (Xa) and inactivated (Xi) state in the dynamics of a single X chromosome. We then proceed to a thorough analysis of 4D and 6D GRN models representing the joint dynamics of the X chromosome pair to (i) ensure the stability of the mono-inactivated (XiXa) state and then (ii) exclude the stability of the bi-inactivated (XiXi) and bi-activated (XaXa) states. Studying the 6D GRN involves the analysis of toggle-switch models coupled through a state-dependent parameter. Combining proper changes of variables and reparameterization, we manage to study both the transient and asymptotic behavior from a 2D phase plane analysis. We further restrict the parameter space to meet quantitative specifications on the relative gene expression level between the Xiand Xa states.
The Kuramoto model is a classical nonlinear ODE system designed to study synchronization phenomena. Each equation represents the phase of an oscillator, and the coupling between them is determined by a graph. There is an increasing interest in understanding the relation between the graph topology and the spontaneous synchronization of the oscillators. Abdalla, Bandeira, and Invernizzi [SIAM J. Appl. Dyn. Syst., 23 (2024), pp. 779--790] considered random geometric graphs on the d-dimensional sphere and proved that the system synchronizes with high probability as long as the mean number of neighbors and the dimension d go to infinity. They posed the question about the behavior when d is small. In this paper, we prove that synchronization holds for random geometric graphs on the two-dimensional sphere, with high probability as the number of nodes goes to infinity, as long as the initial conditions converge to a smooth function. We conjecture a similar behavior for more general simply connected closed Riemannian manifolds, but we expect global synchronization to fail if the manifold is not simply connected, as was shown in De Vita, Bonder, and Groisman [SIAM J. Appl. Dyn. Syst., 24 (2025), pp. 1--15] and suggested in Cirelli et al. [SIAM J. Appl. Math., 85 (2025), pp. 1719--1748].
The usual method to analyze the stability of hetero clinic cycles and networks is with transition matrices that are derived from return maps. In this paper, we introduce an extension to this methodology, the projected map, which we define by identifying trajectories that have, in a certain sense, qualitatively the same dynamics. The projected map is a discrete, piecewise-smooth map of one dimension fewer than the order of the transition matrix. We use these maps to describe the dynamics of trajectories near three hetero clinic networks in R4 with four equilibria. We find in all three cases that the onset of trajectories that switch between cycles of the network is caused by either a fold bifurcation or a border-collision bifurcation, where fixed points of the map no longer exist in the corresponding function's domain of definition. We are able to show that a given initial condition near any of these three networks is asymptotic to only one sub cycle and cannot switch between sub cycles multiple times, resolving a 30-year-old claim by Brannath. We are also able to generalize certain results to all quasi-simple networks, proving that a border-collision bifurcation of the projected map, corresponding to a condition on the eigenvectors of certain transition matrices, causes a cycle to lose stability.
Abstract. The usual method to analyze the stability of heteroclinic cycles and networks is with transition matrices that are derived from return maps. In this paper, we introduce an extension to this methodology, the projected map, which we define by identifying trajectories that have, in a certain sense, qualitatively the same dynamics. The projected map is a discrete, piecewise-smooth map of one dimension fewer than the order of the transition matrix. We use these maps to describe the dynamics of trajectories near three heteroclinic networks in [Formula: see text] with four equilibria. We find in all three cases that the onset of trajectories that switch between cycles of the network is caused by either a fold bifurcation or a border-collision bifurcation, where fixed points of the map no longer exist in the corresponding function’s domain of definition. We are able to show that a given initial condition near any of these three networks is asymptotic to only one subcycle and cannot switch between subcycles multiple times, resolving a 30-year-old claim by Brannath. We are also able to generalize certain results to all quasi-simple networks, proving that a border-collision bifurcation of the projected map, corresponding to a condition on the eigenvectors of certain transition matrices, causes a cycle to lose stability.
We explore the potential of feedback control to stabilize and detect unstable steady propagation modes for confined, deformable bubbles moving through a fluid-filled Hele-Shaw channel. We use a depth-averaged model as a numerical experiment to develop and test the capabilities of control-based continuation (CBC) in this system, with our work being the first use of CBC in free-surface fluid dynamics. We construct a suitable stabilizing feedback gain, with control actuation delivered via fluid injection at the sides of the channel and determined from real-time observations of the bubble shape as viewed from above. We show how this feedback gain can be used to stabilize, detect, and observe unstable states, without prior knowledge of their behavior. As the steady states in fact involve steady bubble propagation along the channel length, we investigate both the idealized case of comoving actuators, as well as the more plausible scenario when actuation is delivered through an array of equally-spaced injection points.
In this paper, a new theoretical approach of gait generation is presented for the legged robots by using delay-coupling oscillators. First, a general theory of the central pattern generator (CPG) is constructed through the time delay coupled oscillators with a unidirectional ring network structure. Based on the symmetric Hopf bifurcation, the parameter condition for rhythm generation is obtained. The parameter group is divided into different regions by Hopf bifurcation curves, where the periodic solutions exhibit distinct spatiotemp oral patterns in each region. Furthermore, to investigate the impact of these spatiotemp oral patterns on the gait generation of multilegged robots, the proposed theoretical approach is applied to simulate various gaits for hexapod and biped robots using six FitzHugh-Nagumo (FHN) oscillators and two Stuart-Landau (SL) oscillators, respectively. The findings reveal that the rhythm signals generated by the CPG controller can prompt the robot to exhibit multiple stable gaits. The CPG controllers, composed of four and eight Van der Pol (VDP) oscillators, are designed to control the hip and knee joint movements of quadruped robots, which achieves the coordinated movement of a robot's leg joints. To this end, a large number of numerical simulations are provided to validate the correctness of the proposed theoretical approach.
We investigate the classical model of competition of two populations in the chemostat when a linear coupling between the populations is taken into account and the removal rates of the populations are distinct from the dilution rate and their yield coefficients are also distinct. This model extends a model of wall growth and a model of lateral gene transfer, previously studied in the literature. We show the existence and uniqueness of the coexistence equilibrium at which the populations coexist, provided that the input concentration of the chemostat exceeds a critical value, or, equivalently the dilution rate does not exceed a critical value that can be computed explicitly. In contrast with the particular cases of this model, previously studied in the literature, the positive equilibrium can be unstable with the appearance of Hopf bifurcations and sustainable oscillations. We construct the operating diagram of the system, which is the two-parameter bifurcation diagram with respect to the operating parameters, that are the dilution rate of the chemostat and its input nutrient concentration. This study reveals a rich variety of dynamical behaviors, including the emergence and disappearance of stable and unstable limit cycles through Hopf bifurcations and limit point of cycles bifurcations. Furthermore, codimension-two bifurcations such as cusp and generalized Hopf points are identified, highlighting complex transitions in the system's dynamics.
Phase resetting, involving shifting the timing of oscillations by a prespecified amount, is a fundamental control problem arising in the investigation of oscillatory dynamical systems. This work investigates the relationships between phase resetting of the collective rhythm and desynchronization in populations of coupled oscillators in the context of an optimal control problem. Conditions on the emergence of a supercritical Hopf bifurcation are considered for a general model of coupled oscillators, and the Hopf normal form is used to simplify the analysis. Two solution archetypes for phase resetting emerge: weak resetting, which does not influence interoscillator synchrony, and strong resetting, whereby the network is first desynchronized before resynchronizing with the correct phase. Analytical expressions for the energy expenditure associated with these solutions yield general conditions for which desynchronization is expected during the course of phase resetting. In particular, strong phase resetting is more efficient than weak resetting when the following conditions are met: 1) the relaxation rate of the collective rhythm to its limit cycle is sufficiently slow relative to its natural frequency, and 2) the required phase shifts are sufficiently large.
In the present work, we investigate the canard explosion occurring in a neuro dynamics model-the van der Pol system with a time-delayed feedback. We assume that the feedback gain is small (on the same order as the time scale difference) such that the delay differential equation (DDE) can be reduced to a planar ordinary differential equation (ODE) on a non-local center manifold. We then expand the ODE flow on this center manifold in the small parameter and apply a nonlinear time transformation to obtain high-order approximations of the critical value and the critical manifold simultaneously. A significant advantage of the present method is that the tedious computations of the center manifold reduction with normal form usually involved in solving delay differential equations are avoided. We prove that each perturbation order of the critical manifold can be expressed as a polynomial of a spatial variable. This greatly simplifies the computations and makes high-order computations possible. We also compare the proposed approach with the existing small-delay expansion of which DDEs are reduced to ODEs by assuming a small delay value. While the latter cannot predict the critical manifold due to the discontinuity, our approach provides accurate and continuous approximations. More significantly, the analytical results obtained by the nonlinear time transformation method go far beyond the existing first-order results, showing an excellent agreement with the numerical simulations even for large delay values. The procedure developed in this work is efficient but simple. Because our method predicts both the critical value and the critical manifold accurately, it has great potential for practical applications in dynamical systems with time-delayed coupling and small coupling coefficients, such as controlling the occurrence of the neuronal spiking and its amplitude, in the future.
We develop an early-warning signal for bifurcations of one-dimensional random difference equations with additive bounded noise, based on the asymptotic behaviour of the stationary density near a boundary of its support. We demonstrate the practical use in numerical examples.
The transition from rotational to discontinuous behavior of the return map of the perturbed oscillators-step system, a paradigm model for a perturbation of a pseudo-integrable Hamiltonian impact system, is studied. The form of the return map is derived, and a truncated form of this map is simulated and analyzed. For a set of parameters the existence of a hovering set, a set of non-resonant orbits that pass sometimes above the step and sometimes to its side, without ever impacting it, is established and quantified. Its destruction as the sign of the perturbation term is reversed is established.
Complex patterns emerge across a wide range of biological systems. While such patterns often exhibit remarkable robustness, variation and irregularity exist at multiple scales and can carry important information about the underlying agent interactions driving collective dynamics. Many methods for quantifying biological patterns focus on large-scale, characteristic features (such as stripe width or spot number), but questions remain on how to characterize messy patterns. In the case of cellular patterns that emerge during development or regeneration, understanding where patterns are most susceptible to variability may help shed light on cell behavior and the tissue environment. Motivated by these challenges, we introduce methods based on topological data analysis to classify and quantify messy patterns arising from agent-based interactions, by extracting meaningful biological interpretations from persistence barcode summaries. To compute persistent homology, our methods rely on a sweeping-plane filtration which, in comparison to the Vietoris–Rips filtration, is more rarely applied to biological systems. We demonstrate how results from the sweeping-plane filtration can be interpreted to quantify stripe patterns (with and without interruptions) by analyzing in silico zebrafish skin patterns, and we generate new quantitative predictions about which pattern features may be most robust or variable. Our work provides an automated framework for quantifying features and irregularities in spot and stripe patterns and highlights how different approaches to persistent homology can provide complementary insight into biological systems.
In this article, we extend the framework developed previously to allow for rigorous proofs of existence of smooth, localized solutions in semi-linear partial differential equations possessing both space and non-space group symmetries. We demonstrate our approach on the Swift-Hohenberg model. In particular, for a given symmetry group 𝒢, we construct a natural Hilbert space H^l_𝒢 containing only functions with 𝒢-symmetry. In this space, products and differential operators are well-defined allowing for the study of autonomous semi-linear PDEs. Depending on the properties of 𝒢, we derive a Newton-Kantorovich approach based on the construction of an approximate inverse around an approximate solution, u_0. More specifically, combining a meticulous analysis and computer-assisted techniques, the Newton-Kantorovich approach is validated thanks to the computation of some explicit bounds. The strategy for constructing u_0, the approximate inverse, and the computation of these bounds will depend on the properties of 𝒢. We demonstrate the methodology on the 2D Swift-Hohenberg PDE by proving the existence of various dihedral localized patterns. The algorithmic details to perform the computer-assisted proofs can be found on Github.
We develop a method for computing the stochastic wave speed of pulse solutions in kinematic equations subject to small stochastic forcing based on the isochronal phase reduction. These kinematic equations arise as the singular limit of sharp pulse solutions in the FitzHugh-Nagumo system, and our approach contributes a new perspective and method to the growing body of work on stochastic wave propagation in excitable media. The method yields an effective Itô process for the wave's position. The coefficients of the Itô process can be computed deterministically allowing for efficient computation. We demonstrate the efficiency and accuracy of our method through numerical demonstrations.
We present an efficient and validated method for approximating the stationary measures of random dynamical systems with smooth additive noise. The approach leverages the strong regularizing properties of the associated transfer operator through a finite-dimensional reduction based on Fourier approximation. Explicit error bounds make the method suitable for use in computer-assisted proofs and rigorous numerical investigations; in particular, its efficiency enables systematic exploration of parameter space. The method provides access to the stationary measure and supports the analysis of key statistical properties of the system. As an application, we study noise-induced phenomena, focusing on the transition from positive to negative Lyapunov exponent---one of the telltale signs of a physical phenomenon known as noise-induced order---in families of random unimo dal maps with Gaussian additive noise. By analyzing the Lyapunov exponent as a function of the system parameters, we identify transitions along a hypersurface in parameter space. The parameters we consider include the standard deviation (intensity) of the Gaussian noise and the shape of the unimo dal map.
