Javier Aguilar

Species populations often modify their environment as they grow. When environmental feedback operates more slowly than population growth, the system can undergo boom–bust dynamics, where the population overshoots its carrying capacity and subsequently collapses. In extreme cases, this collapse leads to total extinction. While deterministic models typically fail to capture these finite-time extinction events, we propose a stochastic framework, derived from an individual-based model, to describe boom–bust–extirpation dynamics. We identify a noise-driven, threshold-like behavior where, depending on initial conditions, the population either undergoes a “boom” or is extirpated before the expansion occurs. Furthermore, we characterize a transition between an excitable regime, where most trajectories are captured by the absorbing state immediately after the first bust, and a persistent regime, where most populations reach a metastable state. We show that this transition is governed by the noise strength and the ratio of environmental-to-population timescales. This framework provides a theoretical basis for understanding irreversible transitions in invasive species, plant succession, microbial dynamics, and the elimination of cancerous tumors.

Identifying the network of species interactions is a fundamental step toward understanding ecosystem stability and biodiversity. However, the interpretability of empirical interaction measures remains a major challenge. Experimental estimates frequently exhibit puzzling temporal fluctuations, including sign shifts typically interpreted as transitions between competition and facilitation. Here, we analyze the temporal behavior of pairwise interaction measures to demonstrate that these fluctuations—and apparent shifts in ecological roles—can emerge intrinsically from standard population dynamics, without any underlying change in the actual ecological relationships. We show that inferred interactions are heavily distorted by experimental protocol choices, particularly the duration of observation and microbial growth constraints. By systematically evaluating interactions across timescales, we uncover a principled mechanism to mitigate these biases: short-term measurements reliably isolate direct, pairwise species couplings, whereas longer-term observations inevitably absorb indirect community feedbacks and systemic experimental constraints. By disentangling direct couplings from indirect network effects, our framework provides a robust, timescale-aware approach to interpreting empirical interaction matrices, offering critical quantitative guidance for experimental design and predictive ecosystem modeling.

Robust inference for stochastic dynamical systems is often hampered by sparse sampling and the absence of closed-form likelihoods. We introduce a Monte Carlo path-inference framework that leverages full-path statistics and bridge processes to deliver reliable parameter estimation and model selection from coarsely sampled time series, without requiring analytical solutions. Crucially, we couple mechanistic stochastic models with their inference procedures to quantify how experimental design—specifically, sampling frequency and dataset size—governs estimator precision and model distinguishability. This analysis reveals optimal sampling regimes and sharp, resolution-dependent limits beyond which competing models become empirically indistinguishable. We validate the approach across four disparate systems—trajectories of optically trapped particles, human microbiome dynamics, social-media topic mentions, and forest population time series—recovering parameters and identifying when inference is fundamentally constrained by measurement resolution, thereby clarifying ongoing debates about dominant noise sources in these systems. Together, these results establish path-based Monte Carlo as a practical, general tool for inference and model discrimination in complex systems and provide principled guidelines for designing measurements that maximize information under real-world constraints.

In stochastic processes with absorbing states, the quasi-stationary distribution provides valuable insights into the long-term behavior prior to absorption. In this work, we revisit two well-established numerical methods for its computation. The first is an iterative algorithm for solving the nonlinear equation that defines the quasi-stationary distribution. We generalize this technique to accommodate general Markov stochastic processes, either with discrete or continuous state space, and with multiple absorbing states. The second is a Monte Carlo method with resetting, for which we propose a single-trajectory approach that uses the trajectory's own history to perform resets after absorption. In addition to these methodological contributions, we provide a detailed analysis of implementation aspects for both methods. We also compare their accuracy and efficiency across a range of examples. The results indicate that the iterative algorithm is generally the preferred choice for problems with simple boundaries, while the Monte Carlo approach is more suitable for problems with complex boundaries, where the implementation of the iterative algorithm is a challenging task.

Antimicrobial resistance is a threat to public health with millions of deaths linked to drug-resistant infections every year. To mitigate resistance, common strategies that are used are combination therapies and therapy switching. However, the stochastic nature of pathogenic mutation makes the optimization of these strategies challenging. Here, we propose a two-scale stochastic model that considers the effective evolution of therapies in a multidimensional efficacy space, where each dimension represents the efficacy of a specific drug in the therapy. The diffusion of therapies within this space is subject to stochastic resets, representing therapy switches. The boundaries of the space, inferred from coarser pathogen-host dynamics, can be either reflecting or absorbing. Reflecting boundaries impede full recovery of the host, while absorbing boundaries represent the development of antimicrobial resistance, leading to therapy failure. We derive analytical expressions for the average absorption times, accounting for both continuous and discrete genomic changes using the frameworks of Langevin and master equations, respectively. These expressions allow us to evaluate the relevance of times between drug switches and the number of simultaneous drugs in relation to typical timescales for drug resistance development. To study realistic therapy scenarios, we impose constraints on the number of administered therapies and/or their costs, which reveals nontrivial optimal drug-switching protocols that maximize the time before antimicrobial resistance develops while reducing therapy costs. Finally, we extend the model to consider single-cell heterogeneity to accurately capture the effects of individual mutations that result in drug resistance.

Species interactions—ranging from direct predator–prey relationships to indirect effects mediated by the environment—are central to ecosystem balance and biodiversity. While empirical methods for measuring these interactions exist, their interpretability and limitations remain unclear. Here we examine the empirical matrix of pairwise interactions, a widely used tool, and analyze its temporal variability. We show that apparent fluctuations in interaction strength—and even shifts in interaction signs, often interpreted as transitions between competition and facilitation—can arise intrinsically from population dynamics with fixed ecological roles. Experimental protocols further shape these estimates: the duration of observation and the type of setup in microbial growth studies (e.g., chemostats, batch cultures, or resource conditions) systematically affect measured interactions. Considering interactions across timescales enhances interpretability: short-term measurements primarily capture direct species couplings, whereas long-term observations increasingly reflect indirect community feedback. Taken together, these results establish short-duration inferences, obtained either directly or extrapolated, as a principled way to disentangle direct from indirect interactions. Building on this insight, we propose a model-inference approach that leverages multiple short time series rather than extended longitudinal datasets.

This paper introduces two techniques for computing the distribution of the absorption or first passage time of the drifted Wiener diffusion subject to Poisson resetting times to an upper hard wall barrier and to a lower absorbing barrier. The first method, which we call the Padé-partial fraction approximation, starts with the Padé approximation to the Laplace transform of the first passage time distribution, which is then exactly inverted by means of the partial fraction decomposition. The second method, which we call the multiresolution algorithm, is a Monte Carlo technique that exploits the properties of the Wiener process to generate Brownian bridges at increasing levels of resolution. Our numerical study reveals that the multiresolution algorithm has higher efficiency than standard Monte Carlo, whereas the faster Padé-partial fraction method is accurate in various circumstances and provides an analytical formula. Also, a closed-form exact expression for the expected first passage time is derived.

Ocean currents exhibit strong time dependence at all scales that influences physical and biochemical dynamics. Network approaches to fluid transport permit to address explicitly how connectivity across the seascape is affected by the spatiotemporal variability of currents. However, such temporal aspect is mostly neglected, relying on a static representation of the flow. We here investigate the role of current variability on networks describing physical transport across the Mediterranean basin. We first focus on degree distributions and community structure comparing ensembles of temporal networks that explicitly resolve time dependence and their aggregated, i.e., time-averaged, counterparts. Furthermore, we explore the implications of the two approaches in a simple reaction dispersal model for a generic tracer. Our analysis evidences that aggregation induces structural network changes that cannot be easily avoided, not even introducing a pruning of the aggregated adjacency matrix. We also highlight that, depending on the time scales considered, the importance of the temporal features of the networks can vary significantly. Finally, we find that the tracer evolution obtained from a temporal dispersal kernel cannot be always approximated by aggregated adjacency matrices, in particular during transients of the dynamics.

Approximate numerical methods are one of the most used strategies to extract information from many-interacting-agents systems. In particular, numerical approximations are of extended use to deal with epidemic, ecological and biological models, since unbiased methods like the Gillespie algorithm can become unpractical due to high CPU time usage required. However, the use of approximations has been debated and there is no clear consensus about whether unbiased methods or biased approach is the best option. In this work, we derive scaling relations for the errors in approximations based on binomial extractions. This finding allows us to build rules to compute the optimal values of both the discretization time and number of realizations needed to compute averages with the biased method with a target precision and minimum CPU-time usage. Furthermore, we also present another rule to discern whether the unbiased method or biased approach is more efficient. Ultimately, we will show that the choice of the method should depend on the desired precision for the estimation of averages.

This article analyzes and compares two general techniques of rare event simulation for generating paths of Markov processes over fixed time horizons: exponential tilting and stochastic bridge. These two methods allow us to accurately compute the probability that a Markov process ends within a rare region which is unlikely to be attained. Exponential tilting is a general technique for obtaining an alternative or tilted sampling probability measure, under which the Markov process becomes likely to hit the rare region at terminal time. The stochastic bridge technique involves conditioning paths towards two endpoints: the terminal point and the initial one. The terminal point is generated from some appropriately chosen probability distribution that covers well the rare region. We show that both methods belong to the class of importance sampling procedures by providing a common mathematical framework of these two conceptually different methods of sampling rare trajectories. We also conduct a numerical comparison of these two methods, revealing distinct areas of application for each Monte Carlo method, where they exhibit superior efficiency. Detailed simulation algorithms are provided.