Impact of Memory on Complex Dynamics Case Studies in Ecology and Epidemiology Moein Khalighi TURUN YLIOPISTON JULKAISUJA – ANNALES UNIVERSITATIS TURKUENSIS SARJA – SER. F OSA – TOM. 64 | TECHNICA – INFORMATICA | TURKU 2025 – – – – IMPACT OF MEMORY ON COMPLEX DYNAMICS Case Studies in Ecology and Epidemiology Moein Khalighi TURUN YLIOPISTON JULKAISUJA ANNALES UNIVERSITATIS TURKUENSIS SARJA SER. F OSA TOM. 64 | TECHNICA INFORMATICA | TURKU 2025 University of Turku Faculty of Technology Department of Computing Computer Science Doctoral Programme in Technology Supervised by Professor, Leo Lahti University of Turku Finland Doctor, Guilhem Sommeria-Klein University of Bordeaux France Doctor, Ville Laitinen University of Turku Finland Reviewed by Professor, Ville Mustonen Assistant Professor, Denis Patterson University of Helsinki Durham University Finland United Kingdom Opponent Doctor, Katharine Coyte University of Manchester United Kingdom The originality of this publication has been checked in accordance with the University of Turku quality assurance system using the Turnitin OriginalityCheck service. ISBN 978-952-02-0344-3 (PRINT) ISBN 978-952-02-0345-0 (PDF) ISSN 2736-9390 (PRINT) ISSN 2736-9684 (ONLINE) Painosalama, Turku, Finland 2025 UNIVERSITY OF TURKU Faculty of Technology Department of Computing Computer Science MOEIN, KHALIGHI: Impact of Memory on Complex Dynamics Doctoral dissertation, 163 pp. Doctoral Programme in Technology October 2025 ABSTRACT Dynamical systems are tools for modeling many real-world phenomena. However, one crucial but often overlooked aspect of such models is memory effects, where past conditions infuence future dynamics. Traditional models that lack memory are inadequate to represent the historical infuences essential for a broad comprehension of these phenomena. In this thesis, we investigate the infuence of memory on the behavior and dynam- ics of commonly used models. The research involves identifying different types of memory and systems, selecting appropriate modeling frameworks, and integrating memory components to quantify their impact on system dynamics. Emphasizing the importance of memory incorporation, this work focuses on models applied in ecology and epidemiology. We used fractional calculus to introduce memory into differential equation models, allowing the infuence of past events to wane gradually. This approach mirrors many real-world scenarios where more recent events carry greater weight. The study explores the application of fractional derivatives, which are derivatives of non-integer orders, to model long-term memory effects, focusing on differential equation-based models. Additionally, the development of computational tools offers a more targeted approach to analyzing memory-augmented models. The dissertation offers new perspectives on ecological systems by examining the infuence of memory on nonlinear behaviors, system stability, and responses to disturbances. It demonstrates how memory acts as a key factor driving complex system dynamics. In epidemiological systems, the work highlights the impact of memory on disease transmission dynamics, revealing long-term effects and showing how memory can alter outbreak predictions. These fndings emphasize the need for more sophisticated models that incorporate memory. Key fndings, supported by publications, illustrate the practical applications of this approach across both simulated and real-world data. Incorporating memory effects into current dynamical systems modeling practices has the potential to enhance their accuracy and applicability across various disciplines. KEYWORDS: Complex systems, fractional differential equations, memory effects, population dynamics iii TURUN YLIOPISTO Teknillinen tiedekunta Tietotekniikan laitos Tietojenkäsittelytiede MOEIN, KHALIGHI: Impact of Memory on Complex Dynamics Väitöskirja, 163 s. Teknologian tohtoriohjelma lokakuuta 2025 TIIVISTELMÄ Dynaamiset järjestelmät ovat välineitä monien reaalimaailman ilmiöiden mallintamiseen. Yksi ratkaiseva mutta usein unohdettu näkökohta tällaisissa malleissa on kuitenkin muis- tivaikutus, jossa menneisyyden olosuhteet vaikuttavat tulevaan dynamiikkaan. Perinteiset mallit, joista puuttuu muisti, ovat riittämättömiä kuvaamaan historiallisia vaikutuksia, jotka ovat olennaisia näiden ilmiöiden laajan ymmärtämisen kannalta. Tässä tutkielmassa tutkimme muistin vaikutusta yleisesti käytettyjen mallien käyttäy- tymiseen ja dynamiikkaan. Tutkimus käsittää erityyppisten muistien ja järjestelmien tun- nistamisen, sopivien mallinnuskehysten valitsemisen ja muistikomponenttien integroin- nin niiden vaikutuksen kvantifoimiseksi järjestelmän d ynamiikkaan. Korostaen muistin sisällyttämisen tärkeyttä tämä työ keskittyy ekologiassa ja epidemiologiassa sovellettuihin malleihin. Käytimme murtolaskennan avulla muistin sisällyttämistä differentiaaliyhtälömalleihin, jolloin menneiden tapahtumien vaikutus voi vähitellen hiipua. Tämä lähestymistapa heijastaa monia reaalimaailman skenaarioita, joissa viimeaikaisilla tapahtumilla on suurempi pain- oarvo. Tutkimuksessa tutkitaan murtojohdannaisten eli ei-kokonaislukuisten derivaattojen soveltamista pitkän aikavälin muistivaikutusten mallintamiseen keskittyen differentiaaliy- htälöihin perustuviin malleihin. Lisäksi laskennallisten työkalujen kehittäminen tarjoaa kohdennetumman lähestymistavan muistin vaikutuksesta vahvistettujen mallien analysoin- tiin. Väitöskirja tarjoaa uusia näkökulmia ekologisiin järjestelmiin tarkastelemalla muistin vaikutusta epälineaariseen käyttäytymiseen, järjestelmän vakauteen ja vasteisiin häiriöille. Se osoittaa, miten muisti toimii avaintekijänä, joka ohjaa monimutkaisten järjestelmien dynamiikkaa. Epidemiologisissa järjestelmissä työ korostaa muistin vaikutusta tautien leviämisdynamiikkaan, paljastaa pitkän aikavälin vaikutuksia ja osoittaa, miten muisti voi muuttaa taudinpurkausten ennusteita. Nämä havainnot korostavat, että tarvitaan kehittyneem- piä malleja, jotka sisältävät muistin. Keskeiset havainnot, joita julkaisut tukevat, havainnollistavat tämän lähestymistavan käytännön sovelluksia sekä simuloidussa että todellisessa aineistossa. Muistivaikutusten sisällyttäminen nykyisiin dynaamisten järjestelmien mallinnuskäytäntöihin voi parantaa niiden tarkkuutta ja sovellettavuutta eri tieteenaloilla. AVAINSANAT: kompleksiset järjestelmät, murtolukudifferentiaaliyhtälöt, muistivaikutukset, populaatiodynamiikka iv Acknowledgements I express my deepest gratitude to my supervisor and research director, Prof. Leo Lahti, for his invaluable guidance and support during this challenging yet rewarding journey. His expertise in diverse research areas, along with his ongoing encouragement, has profoundly shaped my thesis and professional development. I also wish to thank my second supervisor, Dr. Guilhem Sommeria-Klein, my third supervisor, Dr. Ville Laitinen, and the thesis pre-examiners, Prof. Ville Mustonen and Dr. Denis Patterson, whose detailed feedback signifcantly improved my doctoral dissertation. Their knowledge and constructive critiques were instrumental in refning my research. I am very grateful to both past and present members of the Turku Data Science group. Their contributions created a supportive and collaborative atmosphere, and their camaraderie made the challenges of my doctoral studies more manageable. I extend my thanks to all collaborators—working alongside such experts has been both an honor and a tremendous learning experience. A special thank you goes to my friends, who helped me maintain a healthy balance between work and life; their friendship has been a great source of strength. To my cousins and lifelong friends, Mohammad, Kourosh, and Dariush, thank you for your constant support. I owe a deep debt of gratitude to my family—my parents, Nasrin and Davood, and my brothers, Babak, Siamak, Majid, and my twin, Mohammad. Their endless love, understanding, joy, and support have been the foundation of my journey. Lastly, I express my deepest love and thanks to my beloved Razi, whose artistry vividly paints the essence of my ideas and whose presence in my life is a wellspring of boundless empathy and profound meaning. October 15, 2025 Moein Khalighi v Table of contents Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . v Table of contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of original publications . . . . . . . . . . . . . . . . . . . . . . ix 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 Population dynamics modeling . . . . . . . . . . . . . . . . . . 4 2.1 Dynamical systems . . . . . . . . . . . . . . . . . . . . . . . 4 2.1.1 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1.2 Alternative stable states . . . . . . . . . . . . . . . . 5 2.1.3 Resilience and resistance . . . . . . . . . . . . . . . 7 2.1.4 Perturbation . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.5 Tipping points and critical transitions . . . . . . . . . 8 2.1.6 Hysteresis . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1.7 Potential landscapes . . . . . . . . . . . . . . . . . . 9 2.2 Ecological models . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.1 Univariate models . . . . . . . . . . . . . . . . . . . . 12 2.2.2 Multivariate models . . . . . . . . . . . . . . . . . . . 14 2.2.3 Stability analysis . . . . . . . . . . . . . . . . . . . . 19 2.3 Epidemiological models . . . . . . . . . . . . . . . . . . . . . 27 2.3.1 Key quantities . . . . . . . . . . . . . . . . . . . . . . 29 2.3.2 Extended models . . . . . . . . . . . . . . . . . . . . 32 3 Memory effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.1 Basic concepts . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.1.1 Memoryless and memory-integrated dynamics . . . 39 3.1.2 Embedding memory in models . . . . . . . . . . . . 41 3.2 Memory in interacting systems . . . . . . . . . . . . . . . . . 44 3.2.1 Individual component memory . . . . . . . . . . . . 44 3.2.2 Coherent and varied memory . . . . . . . . . . . . . 44 vi 3.3 Modeling memory with fractional calculus . . . . . . . . . . . 46 3.3.1 Fractional calculus . . . . . . . . . . . . . . . . . . . 46 3.3.2 Caputo derivative . . . . . . . . . . . . . . . . . . . . 47 3.3.3 Stability regions . . . . . . . . . . . . . . . . . . . . . 49 3.3.4 Implementation . . . . . . . . . . . . . . . . . . . . . 50 3.4 Alternative concepts . . . . . . . . . . . . . . . . . . . . . . . 52 3.5 Alternative approaches . . . . . . . . . . . . . . . . . . . . . 53 4 Impact of memory on ecological systems . . . . . . . . . . . 55 4.1 Univariate systems . . . . . . . . . . . . . . . . . . . . . . . . 55 4.1.1 Memory and stability landscapes . . . . . . . . . . . 55 4.1.2 Stability landscape transformations . . . . . . . . . . 57 4.1.3 Impact of memory on resilience and resistance . . . 60 4.2 Multivariate systems . . . . . . . . . . . . . . . . . . . . . . . 62 4.2.1 Impact of memory on a model in the relaxed phase 63 4.2.2 Impact of varied memory on a model under pertur- bation . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.2.3 Exploring the frontiers of ecological memory . . . . 66 5 Impact of memory on epidemiological systems . . . . . . . 69 5.1 Impact of memory on a simple model . . . . . . . . . . . . . 69 5.1.1 Memory and fnal epidemic size . . . . . . . . . . . . 69 5.1.2 Memory and basic reproduction number . . . . . . . 70 5.1.3 Challenges in using varied memory . . . . . . . . . 71 5.2 Assessing memory in extended models . . . . . . . . . . . . 72 5.2.1 Ebola transmission dynamics . . . . . . . . . . . . . 73 5.2.2 COVID-19 transmission dynamics . . . . . . . . . . 74 6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 List of references . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 Original publications . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 vii Abbreviations ARFIMA ARIMA L-BFGS �� �0 EWS FDEs FFT � (�, �) (g)LV GPU 1−� ��0 HMC JIT LV ODEs REPL RMSD ℛ0 SEIR SINAR SIR Autoregressive fractionally integrated moving average Autoregressive integrated moving average Limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm Time-fractional Caputo derivative of order � ∈ R+, starting from the initial value �0 Early warning signals Fractional differential equations Fast Fourier transform A derivative function with variable � as a function of time � (generalized) Lotka-Volterra model Graphics processing unit Riemann-Liouville fractional integral of order 1 − � from initial value �0 Hamiltonian Monte Carlo Just-in-time Lotka-Volterra Ordinary differential equations Read-evaluation-print loop Root mean square deviation Basic reproduction number Susceptible-exposed-infected-recovered Sparse identifcation of nonlinear autoregressive Susceptible-infected-recovered viii List of original publications This dissertation is based on the following original publications, which are referred to in the text by their Roman numerals: I Moein Khalighi, Giulio Benedetti, Leo Lahti. Algorithm 1047: FdeSolver, A Julia package for solving fractional differential equations. ACM Transactions on Mathematical Software, 2024; 50(3): 22. II Moein Khalighi*, Guilhem Sommeria-Klein*, Didier Gonze, Karoline Faust, Leo Lahti. Quantifying the impact of ecological memory on the dynamics of interacting communities. PLoS Computational Biology, 2022; 18(6): e1009396. III Faïçal Ndaïrou, Moein Khalighi, Leo Lahti. Ebola epidemic model with dynamic population and memory. Chaos, Solitons & Fractals, 2023; 170: 113361. IV Moein Khalighi, Faïçal Ndaïrou, Leo Lahti. Epidemic transmission model- ing with fractional derivatives and environmental pathogens. International Journal of Biomathematics, 2024; 2450085. The original publications have been reproduced with the permission of the copyright holders. *Equal contribution ix List of co-authored publications not included in the thesis The following peer-reviewed publications were published during doctoral studies. Although they are not part of this thesis, they bear relevance to it. - Moein Khalighi, Leo Lahti, Faïçal Ndaïrou, Peter Rashkov, Delfm FM Tor- res. Fractional modelling of COVID-19 transmission incorporating asymp- tomatic and super-spreader individuals. Mathematical Biosciences, 2025; 380: 109373. - Hamid Safdari, Majid Rajabzadeh, Moein Khalighi. LDG approximation of a nonlinear fractional convection-diffusion equation using B-spline basis functions. Applied Numerical Mathematics, 2022; 171: 45-57. - Hamid Safdari, Majid Rajabzadeh, Moein Khalighi. Solving a non-linear fractional convection-diffusion equation using local discontinuous Galerkin method. Applied Numerical Mathematics, 2021; 165: 22-34. - Moein Khalighi, Leila Eftekhari, Soleiman Hosseinpour, Leo Lahti. Three- species Lotka-Volterra model with respect to Caputo and Caputo-Fabrizio fractional operators. Symmetry, 2021; 13 (3): 368. - Moein Khalighi, Mohammad Amirian Matlob, Alaeddin Malek. A new approach to solving multiorder time-fractional advection–diffusion–reaction equations using BEM and Chebyshev matrix. Mathematical Methods in the Applied Sciences. 2020; 44: 2964–2984. x 1 Introduction Dynamical systems modeling facilitates understanding the natural world’s complex dynamics. Real-world examples, such as the growth of the ecosystem population and the spread of diseases, demonstrate the demand for quantitative tools or models to describe and predict these changes. Existing models can effectively represent complex phenomena, especially when data is limited. However, the increasing availability of data calls for exploring more sophisticated approaches that account for additional complexities. As modeling techniques evolve, simplifed representations often overlook crucial factors such as historical infuences. Many mathematical models rely on memoryless frameworks, assuming that future states depend solely on present conditions. However, research increasingly shows that a system’s history affects how it responds to new situations [1; 2; 3]. This history forms a lasting imprint that gradually diminishes, infuencing future behavior and observations, known as memory effects. Therefore, the primary goal of this dissertation is to examine the impact of memory on the dynamics of complex systems. Ecological systems exemplify complex systems, where memory plays a key role in shaping how organisms and ecosystems respond to disturbances. This is evident in a variety of ecological contexts, from the adaptability of the human gut microbiome to larger ecosystems [4; 5; 6; 7; 8]. By integrating these historical infuences into models, we not only broaden our knowledge of biological and ecological systems but also improve strategies for public health and ecosystem management [5; 9; 10]. This approach elucidates the patterns and mechanisms that drive species adaptation and survival, highlighting the importance of historical legacies in ecological contexts at both micro and macro scales [8; 11]. A major challenge with models that incorporate memory effects is that they require more computational effort and are more complex than models without memory. Thus, this topic has been relatively unexplored in many areas such as social dynamics and microbiome studies, where there is a distinct lack of methodologies for incorporating these effects into standard models. Therefore, I quantifed the impact of memory on complex dynamics by systematically incorporating memory into standard models. Advancing research in memory effects and dynamical models requires a focused approach that identifes specifc types of memory infuences and the most suitable mod- eling techniques. Various felds have developed diverse methodological approaches to incorporating memory into models [12; 13; 14; 15], each offering unique advantages 1 Moein Khalighi and challenges. To address the complexities inherent in these systems, I have aimed to adopt a fexible mathematical framework capable of handling deterministic models. Such a framework should accurately model complex systems and effectively capture the impact of memory. Therefore, I refned the existing models by using fractional calculus, thus enhanc- ing their fexibility. Fractional calculus, which extends traditional calculus to allow differentiation and integration of arbitrary orders [16; 17], is a promising tool for modeling non-standard behaviors and incorporating memory effects into computa- tional modeling. In this regard, fractional derivatives offer a robust mathematical foundation for modeling memory, with a well-established theoretical framework and a growing body of analytical and numerical methods. This provides a strong basis for the methodology and allows for rigorous analysis and interpretation of memory-dependent phenomena [18]. The applicability of fractional derivatives spans various felds and real-world applications [19], as recent studies highlight its role in studying the dynamics of infectious diseases [20; 21]. Fractional-order models capture memory effects and long-range dependencies more accurately, enhancing epidemic simulations compared to classical models [22]. Consequently, I focused on integrating memory processes into ecological and epidemiological models, utilizing fractional derivatives to enhance the simulation of these complex systems. Hence, the research objectives can be summarized to address the following questions: 1. How can memory be incorporated into models of complex dynamics? This question addresses the conceptual and methodological integration of memory into existing frameworks. It involves defning memory, identifying suitable case studies, and systematically developing approaches for memory inclusion. 2. How do memory effects alter the dynamics of ecological models? Here, the focus is on elucidating the impact of historical dependencies on system dynam- ics, with particular attention to aspects such as stability, resilience, resistance, and transitional behaviors in response to disturbances. 3. How do memory effects alter the dynamics of epidemiological models? This question examines how historical dependencies might infuence the dynamics of epidemiological systems, which typically involve nonlinear behaviors, het- erogeneous connectivity patterns, and dynamic network structures [23]. By systematically incorporating memory effects into epidemiological models, I aimed to investigate their impact on the ability of models to simulate and rep- resent real-world epidemic data. Specifcally, I explored whether considering memory leads to meaningful differences in model performance, potentially revealing new patterns, behaviors, or avenues for future research. 2 Introduction 4. How can the incorporation of memory improve classical models? The emphasis here is on the theoretical contributions of memory-based modeling, elucidating discrepancies between classical and memory-integrated models. In addressing these questions, I refned existing population dynamics models, de- scribed by ordinary differential equations (ODEs), using fractional derivatives that span several publications. I began by developing a software package for incorporating memory effects into typical models in Publication I. Publication II delves into the impact of memory on ecological and biological models. Publications III and IV focus on simulating the spread of pandemic diseases with memory effects and ftting these models to real data. The structure of the thesis is designed to frst elucidate the char- acteristics of dynamical systems in Chapter 2, explore the concept of memory effects and the implementations in Chapter 3, and then discuss the impact of memory on various complex dynamics in Chapters 4 and 5. The thesis is concluded in Chapter 6. 3 2 Population dynamics modeling In this chapter, I introduce phenomena and models that are the focus of this thesis. I explore several felds, each focused on understanding how populations, abundance, or density change and interact over time. This initial exploration connects standard models to their practical applications, and I delve into the complexities of these phenomena to lay the groundwork for the thesis. 2.1 Dynamical systems In this section, I explore behaviors in population dynamics, emphasizing properties such as stability (equilibrium states, where the system rests in the absence of external forcing), resistance (the capacity to withstand perturbations), resilience (the speed of recovery following disruption), and nonlinear phenomena like tipping points (abrupt shifts in system states) and hysteresis (dependence of current states on historical trajectories) [24]. These characteristics determine how ecological systems respond to environmental disturbances and stressors. By examining these properties, I set the stage for exploring how memory effects may alter ecological properties. 2.1.1 Stability Deciphering emergent behaviors in complex systems depends on knowing their stabil- ity, especially with regards to interdependent and nonlinear interactions among their components. Stability is crucial for making a system more resilient to disturbances. Thus, understanding ecosystem stability is essential for creating effective interventions and management policies. Nevertheless, the concept of stability in ecosystems remains controversial, with interpretations varying widely among researchers in the ecological literature. This diversity in defnitions can lead to misunderstandings and ambiguities in ecological dis- cussions. A notable example is the study by Grimm and Wissel [25], which elucidates the many conficting defnitions of stability. Consequently, they proposed system- atically categorizing these defnitions and advocated for precise, context-specifc application of stability-related terms in ecological studies. This approach aims to reduce conceptual confusion and clarify the notion of ecological stability. Exploring specifc ecological contexts, like microbial communities, illustrates how the concept of stability varies depending on the situation. For instance, Gonze 4 Population dynamics modeling Time S ta te S ta te S ta te Time Time (a) (b) (c) Equilibrium point, asymptotic stability, or steady state Perturbation Figure 1. Asymptomatic stability in dynamical systems: (a) Illustrates a dynamical system converging towards a fxed or equilibrium point, indicating stability as there are no changes over time. (b) Displays the system converging towards an equilibrium point with minimal fuctuations, demonstrating near-steady state behavior. (c) Depicts the system’s stable response to minor disturbances, highlighting trajectories that return to equilibrium. et al. discuss this understanding in the complex world of microbial ecology [26]. They defne stability in terms of linear asymptotic stability, permanence, persistence, temporal stability, robustness, and structural stability, each refecting different re- sponses to disturbances and environmental changes. These defnitions highlight the context-dependent nature of stability in ecological studies. Here, I explore the concept of stability in dynamical systems, primarily by study- ing asymptotic stability, or steady states. These states manifest at fxed, or equilibrium, points in the system’s behavior, where there is either no change (Figure 1a) or only mi- nor fuctuations (Figure 1b) over time. This form of stability is characterized not just by the system’s trajectories remaining in proximity after a minor disturbance but also by their eventual convergence back to the equilibrium point, underlining a return to a state of stability, Figure 1c. I employ linear stability analysis and other mathematical tools to assess these conditions, elaborated in Section 2.2.3. This approach enables us to minimize complexity and avoid potential confusion in our stability assessments. 2.1.2 Alternative stable states The state of a system can change in response to variations in certain control parameters, exhibiting different behaviors depending on the nature and magnitude of these changes. In some cases, the system’s equilibrium shifts smoothly, while in others, abrupt transitions can occur. Understanding how a system responds to these changes aims to predict and manage its behavior under varying conditions. One way to analyze these responses is through a bifurcation diagram, which illustrates how a system’s equilibrium states evolve as a function of control parameters or external infuences. These diagrams provide valuable insights into the stability and potential transitions within the system. Such conditions are typically explored by gradually changing key infuences to refect how the system might behave in different real-world situations [24; 26]. Figure 2a presents an example of a bifurcation diagram, demonstrating how equilibrium states smoothly change with varying conditions. 5 Moein Khalighi Some systems exhibit general insensitivity to changes across a broad range of conditions, yet respond dramatically near specifc thresholds. This behavior is illus- trated in Figure 2b, where the system maintains a stable state until it crosses a critical point, beyond which there is a sudden and signifcant shift in state. Such transitions have been observed in various natural systems, including shallow lakes, which can shift from a clear to a turbid state due to nutrient overload [27]. When a critical nutrient threshold is exceeded, the lake rapidly transitions to a turbid state, leading to signifcant ecological impacts such as excessive algal growth and biodiversity loss [24]. However, some dynamical systems go beyond having just one stable state. Such systems with alternative stable states present more complex scenarios, which can be illustrated using a bifurcation diagram (Figure 2c). Here, the equilibrium response curve, infuenced by environmental conditions [28; 27], folds to show two distinct stable states (solid lines) separated by unstable equilibrium points (dashed line). This scenario, where a system can exist in multiple stable states under identical conditions, is called multistability. Figure 2. Bifurcation diagrams of three distinct systems (modifed from [24]): (a) Depicts a dynamical system gradually adapting its equilibrium to changing conditions. (b) Illustrates a system’s relative insensitivity within certain external condition ranges, followed by sudden responses near a threshold. (c) Demonstrates a catastrophe fold in the equilibrium curve, emphasizing the system’s ability to switch between two stable states (solid lines), separated by unstable states (dashed line). Small changes can quickly and signifcantly affect the system’s state, especially at tipping points �1 and �2. Blue arrows in these panels indicate the trajectories of system dynamics, starting from an initial value and converging towards the equilibrium point marked by a blue circle. (d) Depicts mono-stable dynamics, where the system converges to a single equilibrium point. (e) Presents a bistable system, emphasizing how initial conditions determine the convergence to distinct equilibrium points. The light blue curve starts below the critical density (dashed line), or the unstable point, and consequently declines toward the lower equilibrium branch. However, the dark blue curve converges to the upper one due to its initial condition. Multistability refers to a system’s ability to settle into different stable states depending on its initial conditions and external disturbances—or even solely on its initial conditions when the surrounding environment remains unchanged. For example, 6 Population dynamics modeling the climate could alternate between a glaciated state and a predominantly ice-free state under identical external conditions [29]. This characteristic can trigger oscillations or lead to the disappearance of certain components upon reaching critical thresholds. Focusing on population dynamics, a low-density population may either recover its carrying capacity or decline toward extinction. These different outcomes show the importance of initial conditions and external pressures in determining a population’s future. Changing conditions can drive a population toward extinction if its density falls below a critical level. This creates two alternative stable states: one at the population’s carrying capacity and the other at a lower (or even zero) density. Unlike single-stable systems (Figure 2d), in this situation, a population cannot recover if it drops below the critical density (dashed line in Figure 2e), especially in harsh conditions where a certain population size is needed for a favorable environment. Multistability can demonstrate the concepts of critical transitions, hysteresis, and tipping points—key phenomena thoroughly explained within the framework of catastrophic bifurcations [24]. These concepts are further detailed in the subsequent subsections, demonstrating how small changes in system parameters can trigger signifcant, and often irreversible, shifts in the system’s state. However, before delving into these topics, I frst introduce the basic concepts of resilience and resistance. 2.1.3 Resilience and resistance Resilience, a multifaceted term widely used across different disciplines, mainly denotes the capacity of a system to absorb or adapt to disruptions and alterations [24]. Different areas of study use this concept to understand how systems respond to disturbances, each adding their unique perspective [25]. This concept is evident in ecosystems, where alterations may favor certain species, and in societies, where social structures and relationships undergo evolution [24]. Drawing on Hodgson et al. [30], I treat resilience and resistance as two sides of the same coin, respectively, corresponding to the return-time and state-change dimensions of system responses. Here, resilience refers to the speed with which a system returns to its pre-disturbance state after a small shock, shorter return times (or higher recovery rates) imply greater resilience, while resistance denotes the system’s ability to endure a disturbance without shifting to an alternative stable state. Together, these properties clarify both a system’s capacity to fend off change and, when change does occur, to recover from it [31]. In Section 2.2.1, I formalize these defnitions for univariate models, mirroring the bivariate “resistance + recovery” framework advocated in [30]. Such defnitions emphasize the importance of identifying relevant perturbations, enhancing their applicability and signifcance, and aiding their interpretation. This approach allows for the study of how memory affects system resilience and resistance in response to these perturbations. 7 Moein Khalighi 2.1.4 Perturbation In mathematical modeling, understanding ecosystem resilience and resistance involves the analysis of the system’s response to various perturbations. In an ecological context, perturbations refer to temporary changes in environmental conditions that cause a notable deviation in an ecosystem’s structure or function. In this thesis, I impose disturbances by altering one or more parameter values in the studied systems, which allows us to simulate perturbations in response to various environmental factors such as changes in resource availability. For instance, changing parameters, such as birth or growth rates in a microbial ecosystem model [26] or fre ignition and tree regrowth rates in an ecological landscape model [6] can represent these perturbations. This approach is closely related to the theoretical sections of this thesis, as discussed in section 2.2. I categorize perturbations into three main types (detailed in Publication II), each with unique patterns and impacts on ecosystems. These types offer distinct per- spectives for analyzing and understanding the resilience and resistance of ecological models. Types of perturbations Pulse perturbations: These are acute, short-term disturbances that abruptly impact the ecosystem but are relatively infrequent. The ecosystem is allowed time to recover or adjust between these occasional perturbations, enabling us to study the system’s immediate reactions and recovery processes (as shown in Figure 1c). Periodic perturbations: In contrast, periodic perturbations occur at regular intervals, bringing a rhythmic disturbance to the ecosystem. This allows the examination of the ecosystem’s resilience through its responses to consistent and predictable disturbances, which might be seasonal or follow another set pattern. Stochastic perturbations: These perturbations are characterized by their randomness and unpredictability. They introduce uncertainty in the analysis, which is essential to understanding how ecosystems respond to irregular, un- foreseeable disturbances, offering insights into the adaptability and robustness of the ecological systems. 2.1.5 Tipping points and critical transitions In catastrophic bifurcations, gradual changes in conditions or parameters can cause sudden shifts in a system’s stability and behavior, leading to transitions between stable states. These transitions, often sudden and potentially irreversible, occur in many real-world systems [27; 7] and reveal critical thresholds that can fundamentally alter 8 Population dynamics modeling how systems evolve. For example, the catastrophic bifurcation diagram 2c visually maps a system’s response to gradual parameter adjustments, identifying critical points �1 and �2 where small changes can trigger signifcant transformations. At these points, the system cannot smoothly move from the upper to the lower part of the curve. These points represent tipping points—critical thresholds where slight variations in conditions can lead to dramatic shifts in a system’s dynamics. Recognizing early warning signals (EWS) that precede such critical transitions is vital across various disciplines, enabling the prediction and potential mitigation of abrupt changes [29]. EWS are statistical indicators, such as rising autocorrelation and variance, that can be detected before critical transitions, offering insights into a system’s approach toward instability [32]. Identifying and responding to these signals can help prevent irreversible consequences in complex systems, such as those found in ecosystem management and climate regulation [32; 29]. 2.1.6 Hysteresis Hysteresis highlights the structure-dependent nature of system responses, showing that the paths to recovery can be quite different from those that led to a change, due to the system’s internal dynamics and multiple stable states [33; 27]. In critical transitions, hysteresis means that returning to the original state is not straightforward once a system has shifted. The system may need conditions that are different from those that caused the transition in order to return to its initial state. For example, in the scenario shown in Figure 2c, if the system starts in a stable high state near point �1 and then shifts to a lower stable state due to small changes, returning to the original high state might require a substantial effort, possibly pushing the system all the way to the tipping point �2. This complexity makes reversing the effects of critical transitions challenging, emphasizing the importance of anticipating hysteresis when managing natural systems. Moreover, hysteresis is a key indicator of multi-stability, as it prevents the system from returning to its original state even if the initial conditions are restored. 2.1.7 Potential landscapes In classical mechanics, a system’s potential energy landscape directly determines how a particle moves under conservative forces [24]. For instance, a ball placed inside a potential well (such as a cup or valley) will settle at the lowest energy state and oscillate around it if perturbed. This visualization corresponds precisely to the mathematical formulation of potential energy in a given system [34]. In many felds, particularly ecology and dynamical systems, this potential land- scape is also termed a stability landscape to emphasize its role in depicting equilibria and system resilience. 9 Moein Khalighi This concept has inspired similar representations in other felds, often as a heuris- tic way to describe stability, transitions, and system dynamics [28; 35; 36]. However, these applications can range from rigorous derivations to more loosely defned analo- gies. In simple one-dimensional dynamical systems, an exact potential function can often be constructed, as long as the dynamics can be expressed in gradient form. In contrast, for more complex systems, especially those involving multiple dimensions or non-conservative forces, an exact potential function may not exist, and the landscape representation becomes more of a conceptual tool rather than a strict mathematical construct [37]. Figure 3. Stability landscapes illustrating mono- and bistable states in dynamical systems (adapted from [28; 37]). The top section depicts a system’s stability under varying conditions, with valleys representing stable equilibria and hilltops as unstable states. Below, a bifurcation diagram (see Figure 2c) highlights equilibria between tipping points �1 and �2, where flled markers indicate stable states and unflled markers denote unstable states. Changes in conditions can shift the system between states, demonstrating critical transitions and hysteresis. In this analogy (Figure 3), a ball represents the system’s state, with valleys corre- sponding to stable equilibria (attractors) and hilltops to unstable states (repellors). The ball’s position refects the system’s current state, and its motion—driven by external forces such as disturbances or gradual shifts in control parameters (e.g., environmental stressors)—determines transitions between states. A mono-stable system has a single valley, while multi-stable systems feature multiple valleys separated by hilltops. Resilience corresponds to the ball’s tendency to return to a valley after perturba- tions, while resistance quantifes the magnitude of disturbance the system can absorb 10 Population dynamics modeling without transitioning to a new basin of attraction. The metaphor also illustrates critical transitions and hysteresis. Consider a bistable system (purple landscapes in Figure 3) where the ball initially resides in the right valley. As control parameters shift (e.g., increasing stress), the landscape deforms: valleys merge or vanish, reducing stability. Crossing a tipping point (�1) triggers an irreversible transition to a new basin (left valley). Reversing this transition requires surpassing a different tipping point (�2), demonstrating hysteresis—the system resists returning to its original state due to intervening hilltops. Key quantitative features include basin size (width/depth of valleys), attractor curvature (rate of return to equilibrium), and potential depth (resistance to pertur- bations) [38]. These metrics, computationally analyzed in Section 2.2.3, determine dynamical behavior. For example, deeper basins imply greater resilience, while steeper slopes correlate with faster recovery. This framework bridges intuitive visualization and dynamical theory: Landscapes directly map to equations governing system behavior, while basins and potentials quantify resilience [38]. Its utility extends to ecology, climate science, and other felds where multi-stability and tipping points are critical. 2.2 Ecological models The mathematical models in this thesis track the changes in the population over time. Such models are applied widely in ecology, epidemiology, economics, and sociol- ogy to analyze drivers of population change, including birth/death rates, migration, resource availability, competition, predation, and disease dynamics, as well as so- cioeconomic factors like technological advancements and market forces [39]. These systems range in complexity from simple growth equations to nonlinear differential frameworks such as: �� = � (�, �), (1) �� where � is the population changing over time � by equation function � (�, �), denot- ing the forces driving the system’s dynamics. In ecological contexts, � (�, �) encodes growth drivers (e.g., resource abundance) and constraints (e.g., density dependence). While Equation (1) represents a deterministic frst-order ordinary differential equa- tion, other approaches—including partial differential equations, stochastic models, and higher-order systems—are used to address spatial heterogeneity, randomness, or delayed feedbacks [39]. Population dynamics in ecology exemplify the concepts discussed in previous sections. Typically, populations often exhibit rapid initial growth under abundant resources, followed by stabilization as they approach environmental carrying capac- ities [39]. Interdependent cycles, such as predator-prey interactions, highlight how population trajectories are shaped by interspecifc relationships. Ecosystem stability further depends on resilience and resistance to structural shifts caused by gradual 11 Moein Khalighi changes reaching critical thresholds. For instance, the Allee effect demonstrates extinction risks for low-density populations struggling to locate mates [40]. The following sections will explore the typical ecological deterministic models and their properties, depending on the dimensions of these systems. 2.2.1 Univariate models One of the simplest forms of a population dynamics model is the Exponential growth model, generally expressed by the equation function ��(�), where � is the growth rate. Nevertheless, the limitation is that its growth does not stop and eventually goes to infnity. Thus, considering an environment’s carrying capacity (�), limiting how large the population can grow, we have the logistic growth model [41]: �� � = ��(1 − ). (2) �� � Unlike the exponential model, this logistic model stabilizes at a fnite value, aligning with the carrying capacity �. Further modifcations to the logistic growth model introduce additional terms to account for new factors relevant to specifc modeling purposes. For instance, in the study of desertifcation, � represents the vegetation biomass, with infuences such as herbivory and land degradation shaping its dynamics. Desertifcation, primarily driven by overgrazing, deforestation, and unsustainable land use, impacts the carrying capacity �, making it a dynamic variable [42]. Hence, to model the impact of herbivory, an additional term like − �� can be included to represent the reduction of �+� vegetation due to grazing, where � and � represent the herbivore population density and half-saturation density, respectively [43]: �� � �� = ��(1 − ) − . (3) �� � � + � The equilibrium points of the model (3) are found by setting the equation function �� to zero, i.e., ��(1− � )− = 0. One trivial equilibrium is � = 0, corresponding � �+� to the absence of vegetation. Other equilibrium points depend on the values of the parameters. In biologically relevant scenarios, the model can have at least one positive stable equilibrium and one unstable point between these stable equilibria [43]. There- fore, this model could be bistable, highlighting the importance of initial conditions in determining the eventual state of vegetation density. The details for assessing the stability of systems’ equilibrium points are provided in section 2.2.3. Quantifying potential landscapes To study the stability of nonlinear systems, we can use the concept of potential landscapes discussed in Section 2.1.7. The potential energy of the general model (1), 12 Population dynamics modeling serving as its stability landscape, is calculated by the equation: �� �� = � (�, �) = − , (4) �� �� where � (�) is the corresponding potential function. In this framework, the negative derivative of � (�) forms the right-hand side of the differential equation. To determine � (�), we perform the integration of � (�, �) over the state variable �: ∫ � (�) = − � (�, �) ��. (5) This formula is used to analyze the stability of certain systems in Section 2.2.3. I will then theoretically demonstrate how memory infuences this stability landscape in Sections 4.1.1 and 4.1.2. Resilience I defne resilience (or recovery rate) as follows [44]: ( ) |�0−�� |−|�0−�� | |�0−�� |+|�0−�� |Resilience = , �� − �� where �0 is the original stable state of the system, �� is the state at the end of the perturbation (time ��), �� denotes the state at which the system is considered to have recovered suffciently close to �0 (at time ��), see Figure 4. Time x0 xn xp S ta te s Perturbation t0 tp tn Figure 4. Illustration of resilience using a quantitative approach in a one-dimensional system. The system is initially at the stable state �0. A perturbation from time �0 to �� moves the state to ��. After the perturbation ceases, the system recovers, reaching a state �� (within an acceptable proximity of �0) by time ��. Intuitively, this metric quantifes how rapidly the system returns toward its orig- inal state after being perturbed. A higher resilience value indicates that the system 13 Moein Khalighi recovers more quickly relative to the magnitude of the disturbance, as measured by the difference |�0 − ��|. It is important to note that in many ecological systems, the precise stable state �0 and the exact recovery time �� are diffcult to observe because these systems are under continuous perturbation [31]. In practical studies, these quantities are often estimated or substituted with measurable proxies that refect the system’s dynamic response. Resistance I defne resistance as the maximum level of pulse disturbance an ecosystem can absorb without transitioning to an alternative stable state [38]. In other words, resistance measures the inherent ability of an ecosystem to maintain its structure and function when subjected to external perturbations. Intuitively, a system with higher resistance can endure larger disturbances while remaining within the same stability basin. Factors such as biodiversity, disturbance history, and resource availability contribute to this capacity. This concept is particu- larly relevant in deterministic models where a threshold disturbance level can be more clearly identifed [38]. However, as with resilience, it is important to recognize that in real-world ecolog- ical systems, the stable state and its relaxation time are rarely observed directly due to the continuous nature of environmental disturbances [24]. Consequently, resistance is often estimated using surrogate measures or approximated thresholds that refect the ecosystem’s adaptive capacity. 2.2.2 Multivariate models In ecology, systems often consist of numerous interacting species or components, resulting in complex, high-dimensional dynamics that simple univariate or low- dimensional models cannot adequately represent. Multivariate models are essential for understanding these intricate systems and fnd wide application in areas such as conservation, climate change prediction, biodiversity assessment [45], habitat restora- tion, policy formulation [46], invasive species management, resource management, land-use planning, and evaluating ecosystem services [47]. These diverse applications show the role of multivariate models in ecological research. Beyond macroscopic ecology, the feld of microbial ecology has undergone ad- vancements in recent years, driven largely by high-throughput measurement technolo- gies such as metagenomic sequencing [48]. However, practical time series analyses in this domain frequently face challenges, including measurement uncertainties, limited observations, and the complexity of high-dimensional data. Microbial communities offer an ideal experimental platform for tackling these issues [49; 50], as they can be manipulated and observed under controlled conditions, enabling researchers to test theoretical models, such as ODEs [26; 33]. 14 Population dynamics modeling Among these theoretical tools, the Lotka-Volterra (LV) equations stand out for their ability to represent fundamental ecological interactions, such as predation and competition. For instance, research has utilized LV models to investigate the ecolog- ical roles of key species in the human gut microbiome [51], shedding light on the principles that govern these microbial communities. Furthermore, extended mathemat- ical simulations based on these models have demonstrated the existence of alternative stable states within human gut microbial communities [33]. In this section, I will explore the key principles and methodologies for modeling multivariate systems, with a particular focus on LV models and their extensions. These approaches provide frameworks for understanding interactions in ecological, biological, and other complex networks, illuminating both natural and engineered phenomena. Lotka-Volterra : R� For modeling multivariate systems, when � ∈ R� and � → R�, � > 1, LV models are a popular starting point [52]. These are sets of two or more differential equations used to describe predator-prey relationships. In its most basic form, the model can be represented by a pair of equations: �� = �� − ���, �� (6) �� = ��� − ��, �� where � and � represent the prey and predator populations, respectively, while �, �, �, and � are parameters denoting growth rates and interaction effects, such as predation rates. Oscillatory dynamics characterize this model. The versatility of the LV model extends beyond predator-prey dynamics, covering other ecological interactions [53; 26]. A simple modifcation introduces additional factors that could create more equilib- rium points, such as a more complex interaction term or carrying capacities for both variables. The latter is achieved by expanding the terms �� and �� to ��(1−�/�1) and �� (1 − �/�2), respectively, where �1 and �2 represent the carrying capacities of the prey and predator populations: �� � = ��(1 − ) − ���, �� �1 (7) �� � = ��� − �� (1 − ). �� �2 This adaptability enables us to describe species abundance dynamics in the human gut microbiome, parameterized using in vitro interaction experiments [51]. I use it 15 Moein Khalighi to demonstrate the impact of memory on their dynamics in Section 4.2 and explore it in Publication II. Over the years, numerous extensions have been developed to capture a wide range of complexities and scenarios. In Publication II, I adopted one such extension, which is described in further detail in Section 2.2.2. Moreover, these extensions are not confned to biology; they have also been applied in diverse felds—such as economics and social science—to analyze competitive and cooperative dynamics [54; 55]. Interactions Components engage in multifaceted interactions in ecology and sociological structures, constituting complex networks within dynamic communities. These interactions drive community function, offering a framework for analyzing and comprehending the intricate relationships within a network [56]. Through network inference and matrix interaction [57], we can identify important connections between different parts of a community, whether they are organisms in an ecosystem or entities in a social network. Consider, for instance, a generalized Lotka–Volterra (gLV) model that accounts for interactions in a multivariate community. In its simple linear form, the model is expressed as �∑ ��� = ��(�� + ��� �� ), �� �=1 where �� represents the growth rate of the �th component, � > 1 is the number of components, and ��� is the interaction coeffcient between components � and �. This model extends (6) to accommodate models with a large number of variants. The interaction coeffcient ��� covers both intra-specifc interaction coeffcients (� = �), and inter-specifc interaction coeffcients (� ̸= �), respectively. Such a model can be represented in matrix form, compressing the interdepen- dencies and interactions within the community, as Figure 5 demonstrates for a four- components community. The interaction matrix, symbolized by � in Figure 5(b), constructs dynamical models. This matrix can be derived from static network data or direct experimental observations [58], and it forms the foundation for simulating how communities of interacting members behave under different scenarios. Negative and positive values within the matrix interactions signify the nature of pairwise interactions—antagonistic or synergistic, respectively. Moreover, the magnitude of these elements refects the strength or intensity of the interactions, quantifying the impact of one component on another. The matrix’s diagonal elements represent self-interaction, denoting the infuence of components on their own state. These diagonal elements are sometimes omitted for simplicity or illustrated with curved arrows looping back to their origin in graphical representations of the model, 16 Population dynamics modeling Interaction strength Growth rate levels Mutualism Predation Amenalism Commensalism Competition Neutral + + + + _ _ __ 0 0 0 0 (a) Interacting Community (b) Mathematical Model Figure 5. Illustration of a four-component community model within a generalized Lotka-Volterra (gLV) framework. (a) A graphical diagram of the model, highlighting the interplay among components through color-coded arrows: blue for positive and orange for negative interactions. These arrows symbolize the nature of each pairwise interaction. (b) The gLV model’s depiction of community population dynamics integrates all components within the vector �. This dynamics is governed by their growth rates � and the interaction matrix �. The circle sizes in vector � visually represent the growth magnitude of each component. Matrix interaction � quantitatively represents interdependencies among the community members, in which darker hues signify stronger relationships. as shown in Figure 5(a). To understand a model, we need to understand the nature of interactions. Utilizing biological terminologies, the pair-wise interactions are categorized into six fundamen- tal classifcations [58]: Mutualism (where both beneft), Predation (one preys upon the other), Competition (they compete for resources), Amensalism (one is harmed, the other is unaffected), Commensalism (one benefts, the other is unaffected), or Neutral (neither affects the other). By collecting data, we can predict interactions within a community, recognizing that non-random patterns typically arise from specifc factors. Identifying these interactions involves determining positive connections, inferred when two elements co-occur or show similar patterns across multiple observations. Conversely, negative connections are deduced when elements consistently avoid each other. However, interpreting the true signifcance of these interactions presents challenges. Positive connections might result from cooperation or shared preferences, while negative ones could be due to confict, predator-prey dynamics, or competition. Furthermore, these interactions may not be immediate; for example, an increase in one element might only later cause a decline in another. Ongoing research increasingly integrates time-series data, multi-omics approaches, and refned statistical methods to discern direct effects from indirect ones, thereby advancing our understanding of how complex interaction patterns shape community dynamics. 17 Moein Khalighi In computer science, deriving interaction networks from binary or quantitative data is called network inference, a feld with a signifcant history in method develop- ment [57]. These techniques, widely used in genetics, are now increasingly applied in environmental studies [56]. Network inference methods can be categorized into identifying binary interactions and revealing more complex connections. This thesis primarily examines how memory infuences system dynamics using a constant parameter set and interaction matrix. While evaluating the interaction matrix is beyond the scope of this thesis, its role in transmitting memory effects across interacting communities is undeniable. Nonlinear interaction models While the generalized Lotka-Volterra model offers a foundational framework for un- derstanding linear interactions, many ecological systems display nonlinear dynamics that demand more advanced modeling approaches. A key inspiration for such models is a biochemical rate equation formulation of gene expression, such as genetic toggle switch [34], introduced by Gardner et al. [59]. In their work, two repressors mutually inhibit each other, described by the equations �� �1 = − �, �� 1 + � � �� �2 = − �, �� 1 + �� where � and � are repressor concentrations, �1 and �2 are synthesis rates, and � and � denote repression cooperativity. This system exhibits bistability, where one repressor dominates while suppressing the other, a principle extended by Gonze et al. [33] to microbial communities with three or more species. Building on this concept, the nonlinear inhibition model captures inhibitory interactions among species. The rate of change of species abundances �� is given by [33]: � �� ∏ ��� �� = �� (����({��}) − ����) , ��({��}) = , (8)�� �� + �� �� � �=1 � ̸=� where �� is the intrinsic growth rate of species �, �� governs intraspecifc competition (or self-limitation), and ��({��}) is an inhibition function aggregating the effects of other species � ̸= �. The product of Hill functions [60], parameterized by interaction constants ��� and Hill coeffcient �, mirrors the cooperative repression in Gardner’s model, introducing a saturating effect: as a competitor’s abundance �� rises, its inhibition of �� increases but eventually levels off. This formulation enables the model to represent diverse competitive dynamics, with inhibition strength varying based on competitor abundance and species sensitivity 18 Population dynamics modeling (via ���). Unlike the linear terms of the Lotka-Volterra model, the Hill functions provide a realistic depiction of ecological interactions—such as resource competition or toxin production—where effects saturate at high abundances [34; 60]. Gonze et al. [33] adapted this framework to explore microbial communities, showing that mutual inhibition among three species can lead to multi-stability, where one species dominates while maintaining others at low levels, resulting in alternative stable states under identical conditions. Their fndings highlight the model’s ability to replicate observed microbial ecosystem behaviors. In Publication II, I apply this model to investigate how memory infuences the temporal dynamics of microbial abundances under environmental perturbations. 2.2.3 Stability analysis Various methods and metrics are employed to analyze the stability of systems de- scribed by ODEs [39]. These analyses depend on examining equilibrium points, where the system’s rate of change is zero, i.e., ��/�� = 0. The stability of these points, whether attracting or repelling, is critical for assessing the system’s resilience. These equilibrium points are asymptotically stable if the system returns to them after small perturbations, meaning that trajectories converge to the equilibrium over time. Using the function � (�) of a differential equation (1), we derive characteristic equations whose roots help ascertain these stability conditions. The subsequent sections are divided based on system dimensions, starting with the simpler one- dimensional systems and then advancing to the more complex multi-dimensional systems. One-dimensional systems One-dimensional models are an ideal starting point for illustrating fundamental stability concepts as they allow clear analytical and graphical interpretations. To demonstrate basic ideas like stable and unstable equilibria, I begin with a linear exponential growth model, one of the simplest cases that can still exhibit these properties. I then move on to a cubic polynomial model, which introduces multi- stability and catastrophic bifurcations—key phenomena that cannot arise in the purely linear case. Along the way, I will connect these models to the stability landscapes introduced in Section 2.1.7 to show how potential functions can visually express the system’s behavior. Linear exponential growth model We start with one of the simplest one-dimensional systems that can exhibit both stable and unstable equilibria: �� = �(� − �), (9) �� 19 Moein Khalighi �� where is the rate of change of the system, proportional to the deviation of �� the current state � from an equilibrium �. Its linearity allows for straightforward analytical solutions, making it an ideal frst step in demonstrating fundamental stability concepts. For instance, in a population context, if the number of individuals exceeds �, the population grows faster; if it is below �, the population declines toward this baseline. The sign and magnitude of � determine whether the system converges to or diverges from �. This setup can represent scenarios where external conditions enforce a natural equilibrium, leading to stabilization or oscillation around the point �. The solution to this equation is: �(�) = � + (�0 − �)��� , (10) where �0 is the initial value. The top row of Figure 6 illustrates the solutions for different values of �: negative (stable equilibrium), zero (neutral equilibrium), and positive (unstable equilibrium). x dx/dt dx/dt dx/dt x x x t a a a x t a x r = 0r < 0 r > 0 t a Figure 6. System dynamics and phase space behavior. (Top row) Dynamics of �(�) for equation (9) across different initial conditions (gray lines). Panels depict scenarios for the parameter � being negative (left), zero (center), and positive (right), illustrating the system’s response under varying stability conditions. (Bottom row) Corresponding phase space diagrams showing the rate of change �� as a function of � for each parameter value �. �� �� Additionally, we can analyze as a function of � through phase space plots �� �� (shown in the bottom row of Figure 6). These plots depict how the slope of �� corresponds to the system’s stability [61]: 20 Population dynamics modeling �� � < 0: The phase space plot shows a negative slope. Here, if � < �, is positive, �� �� indicating movement towards �; if � > �, is negative, leading the system �� back towards �. This convergence characterizes � as a stable equilibrium point or attractor. � = 0: The plot displays a horizontal line, indicating no change regardless of �’s initial position, categorizing it as a neutrally stable state. � > 0: The plot reveals a positive slope. If � deviates from �, the change in � accelerates away from �, marking the equilibrium as unstable. Phase space plots, such as those in Figure 6, provide visual insights into the dy- namics of the system. They demonstrate how trajectories converge towards or diverge from equilibrium points, marked by solid circles for stable points and empty circles for unstable points. While these visualizations may appear simplistic, they illustrate the behavior of nonlinear dynamical systems in more complex cases, providing initial insights for understanding system behavior [61]. Cubic polynomial model and potential landscape While the linear model is insightful for basic stability analysis, it cannot describe more complex phenomena like catastrophic bifurcations or multiple stable equilibria. To illustrate such behavior as simply as possible, we move to a one-dimensional polynomial model of degree three: �� = � + � − �3 . �� This model is among the simplest systems that can exhibit multiple stable equilibria, bifurcations, and hysteresis, making it an ideal choice for demonstrating nonlinear stability concepts in a clear and tractable way. Since population variables � often cannot be negative, I shift the variable (e.g., � → � − �) to ensure all equilibrium points are non-negative, giving: �� = � + (� − �) − (� − �)3 , (11) �� and using (5), the corresponding potential (or stability) function becomes: � (�) = −�� − 1(� − �)2 + 1(� − �)4 . 2 4 Rather than merely solving these equations, I prioritize comprehending the system behaviors visible in phase space graphs and potential landscape plots (introduced in Section 2.1.7). Therefore, in the following, I assess the equilibrium points in the �� phase space plot of as a function of � . �� 21 Moein Khalighi For � = 0, the system simplifes to: �� = (� − �) − (� − �)3 . �� Phase space plots for this scenario are illustrated in the central diagram of the upper row in Figure 7. Adjusting � from zero alters the graph’s position: decreasing � shifts it downwards while increasing � shifts it upwards. Figure 7. Stability analysis of equation (11) across different parameter values. (Upper plots) Phase space plots showing solution shifts for positive and negative �, with transitions at bifurcation points �1 and �2. (Middle plot) Bifurcation diagram illustrating hysteresis with saddle-node bifurcations at � = �1 and � = �2. (Lower plots) Potential functions (stability landscapes) demonstrating how equilibria change as � varies. Initial states in the left minimum require � > �2 to switch to the right minimum, and states in the right minimum require � < �1 to switch back. Potential depth and curvature at � = 0 for the right equilibrium are highlighted in green. The key consideration is the number of intersections of the graph with the hori- zontal axis, indicating equilibrium points. With � = 0, three equilibrium points exist, and slight variations in � maintain this number. However, specifc values of � reduce the number of equilibrium points to two, corresponding to bifurcation points. These are denoted �1 and �2 for negative and positive critical values. Beyond these critical values, the system has only one equilibrium point. A qualitative change occurs at the parameter values �1 and �2, marking a transition from three equilibrium points to one. 22 Population dynamics modeling Stable and unstable equilibrium points can be identifed from the phase space plot (upper row in Figure 7). Notably, in nonlinear systems, when � = �1 or � = �2, the slope becomes zero—a condition previously termed neutrally stable for Equation (9). However, in Equation (11), examination of the state and changes near the curve’s tangent to the horizontal axis reveals that the fow converges from one side and diverges from the other. These are known as half-stable or saddle points, indicated by half-flled circles [61]. The relationship between the state space, bifurcation diagram, and potential function is demonstrated in Figure 7. This fgure includes a phase plot where the parameter � shifts from less than �1 to greater than �2. Initially, for � < �1, the plot shows a single intersection with the x-axis, indicating a stable equilibrium point. As � increases to �1, a saddle point emerges, leading to a bifurcation into a stable and an unstable equilibrium point—these points correspond to a local minimum and maximum on the potential function respectively. In the range �1 < � < �2, the system exhibits two stable equilibria (minima of the potential function) and one unstable equilibrium (a maximum of the potential function) between them. When � exceeds �2, the left minimum disappears, enabling the system to shift to the only remaining equilibrium point on the right. If we start with � > �2 and decrease the control parameter, the system remains at higher values of � until it crosses �1, demonstrating hysteresis. This phenomenon, where the system’s state depends on its historical parameter changes, is introduced in Section 2.1.6. From these phase space plots, the structure of the potential function as a stability landscape becomes clearer. We can now theoretically explain why some valleys are deeper or shallower, why slopes vary, and why the basins of attraction differ in size. Dakos and Kef [38] thoroughly investigated components of the stability landscape, such as potential depth and curvature. Potential depth measures the distance between the valley bottom and the hilltop. Potential curvature indicates the slope of the potential function at stable points and is quantifed by the derivative of the equation function at these points [38], � ′′ (��) = −� ′ (��), where �� indicates a stable point of the system. This derivative corresponds to the eigenvalue in a one-dimensional Jacobian. These key elements are also illustrated in Figure 7, where � = 0. Multi-dimensional systems This section lays the groundwork for the applications and analyses presented in subsequent chapters, as the models developed in my thesis publications are multi- dimensional. In later sections 3.2 and 3.3, I will examine how memory effects are integrated into these models. Understanding state-space diagrams and eigenvalue- based stability analysis enables the examination of deterministic, memory-integrated models of multi-entity interactions, facilitating the analysis of their properties and the interpretation of memory effects (Publications II-IV and 4.2). 23 Moein Khalighi State space diagrams Similar to the phase space plots introduced for one-dimensional systems, we can use state space diagrams for two-dimensional systems to illustrate the evolution of states. A state space plot of a two-dimensional system, where each axis represents a state variable, is a top-down projection of the system’s movement on a three-dimensional stability landscape [37]. In this landscape, the x and y axes correspond to the state variables, while the z-axis represents potential energy or stability. The trajectory in the state space plot mirrors the path of a ball rolling on the stability landscape, moving towards valleys (stable states) and away from hills (unstable states). Figure 8 shows how state space diagrams visualize dynamical systems and reveal behaviors under various initial conditions. Each axis represents a state variable, illus- trating how state attributes change over time. They highlight stability, multistability, and oscillation patterns, especially in cyclical or seasonal systems. Focusing on the classic predator-prey model (6), its biological implications can be observed. If a neutrally stable point is present in this model, continuous oscillations may occur [39]. This is evident in state space diagrams, where axes represent the predator and prey populations, showing cyclical dynamics as closed loops, Figure 8a. While ideal models suggest perpetual oscillations, real-world complexity (like envi- ronmental limits, time delays, or external factors) can alter these idealized neutral stability scenarios [61]. For example, adding ecological factors can change loop shapes and sizes, refecting the carrying capacity and birth and death rate variations. In predator-prey dynamics, oscillations converging on an equilibrium point in- dicate stabilization after fuctuations. Two primary outcomes are convergence to extinction or coexistence [34]. In the frst scenario, oscillations dwindle until one or both species become extinct, Figure 8b. For instance, if prey numbers fall criti- cally low, they might not recover, leading to their extinction and, subsequently, the predators’, due to food shortage. State space depicts this as trajectories spiraling towards the X-axis (prey extinction) or the origin [42] (both species extinct). In the second scenario, both species stabilize at certain population levels, Figure 8c. Fac- tors like increased carrying capacity, altered birth/death rates, or adaptive behaviors (like predators fnding alternative food sources), can lead to coexistence [39]. State space shows this as trajectories spiraling to a non-axial point, indicating stabilized population levels of both species. Nullclines, representing the loci in a system’s state space where the rate of change of variables zeroes out, are essential in determining dynamical models’ stability regions and trajectories [61]. By marking the equilibrium points where system variables cease to change, nullclines help predict the system’s local behavior near these points, potentially indicating the type of equilibrium (e.g., node, saddle, spiral, center) present. In the context of the LV models, which simulate species competition, including those without oscillatory dynamics, nullclines offer insights into various scenarios [61]. The system’s response, whether converging towards or repelling from equilibrium, 24 Population dynamics modeling Figure 8. Illustrations of various two-dimensional population dynamics in state space: (a) Idealized oscillatory behavior, forming a closed loop indicative of continuous cycles. (b) Oscillatory patterns leading to extinction, with trajectories spiraling towards axes, symbolize population decline. (c) Oscillatory trends stabilize at coexistence, where trajectories spiral to a non-axial equilibrium, signifying balanced population levels. (d) Non-oscillatory behavior with a single nullcline intersection, depicting a stable coexistence point as trajectories converge. (e) Non-oscillatory scenarios without nullcline intersections, forecasting the extinction of one population as trajectories converge towards an axis. (f) Bistable conditions with nonlinear nullclines, marked by three intersection points, denoting one unstable and two stable equilibrium states, refecting complex population interactions. depends on its parameter values and the confguration of the nullclines. For example, Figure 8d depicts nullclines intersecting at a steady state, suggesting a stable equilib- rium where the system’s states converge. Conversely, Figure 8e presents a scenario without nullcline intersection in the biologically relevant quadrant, predicting the extinction of one species and a single steady state along one axis. This highlights the role of nullclines and system parameters in shaping the dynamics and outcomes in models of species competition. In more complex scenarios, the phenomenon of bistability can emerge, as il- lustrated in Figure 8f, through the portrayal of nonlinear nullclines. Three distinct points of nullcline intersection mark this confguration, denoting the convergence of trajectories. One intersection signifes an unstable equilibrium, while the others 25 Moein Khalighi represent stable states. The trajectory and eventual state of the system are heavily infuenced by the initial abundance of each species. In such a bistable framework, one species can gain a competitive edge, not leading to the extinction of the other but rather establishing a state of coexistence where both species persist, albeit with one displaying a dominant presence [26]. This manifestation of multistability, akin to model (8), explored in Publication II, demonstrates the rich dynamics inherent in ecological systems. Furthermore, in Chapter 4.2, I delve into another facet of bistability by studying the impact of memory on model (7) characterized by linear nullclines intersecting at a single point [51]. Interestingly, this model results in the extinction of one species, a scenario potentially attributable to intense competition and the interplay of “loss-loss” interactions between the two variables. Linearization and eigenvalues State space diagrams provide a simplifed view compared to the stability landscape, as they do not convey the velocity of state changes, a feature distinctly shown in the stability landscape. Therefore, a more theoretical approach is necessary to compre- hend the types of equilibrium points and the scenarios described, through analysis rather than solving the models. Furthermore, for high-dimensional systems, the be- havior becomes complex and less intuitive. In this case, a theoretical method involves using linearization and eigenvalue analysis [61]. I apply local linearization around equilibrium points to achieve a linear approximation of the systems. Suppose the function � : R� → R� represents our system, and �* is an equi- librium point, where � (�*) = 0. To examine local stability around this point, we linearize the system using a Taylor expansion [61]: � (�) ≈ � (� * ) + �(� − � * ), where � is the Jacobian matrix: [ ] ��� � = . ��� �=�* Since � (�*) = 0, we have � (�) ≈ �(� − �*). Thus, �� ≈ �(� − � * ). �� An equilibrium point �* is called hyperbolic if all eigenvalues of the Jacobian � have nonzero real parts, i.e., Re(��) ≠ 0 for all �. Hyperbolicity is essential because, by the Hartman–Grobman theorem [62; 63], the local behavior of a nonlinear system near a hyperbolic equilibrium is topologically conjugate to the behavior of its linearized approximation. In practical terms, this means the linearization fully refects the qualitative stability properties, including attraction, repulsion, and saddle-type structures around the equilibrium point. 26 Population dynamics modeling To analyze the stability, we need to study the eigenvalues of the Jacobian matrix � . The solutions to the equation �� ≈ �� are found as �(�) = ����(0), where ��� �� is the matrix exponential, a generalization of the scalar exponential function. For the diagonalizable matrix � , the matrix exponential ��� can be computed using its eigenvalues and eigenvectors. If � can be expressed as � = � �� −1, where � is a diagonal matrix of eigenvalues �� of � , and � is the matrix of corresponding eigenvectors, then ��� = � ���� −1 , where ��� is a diagonal matrix with entries ��� � . Each term ���� indicates the growth or decay of the system’s response along the direction of the corresponding eigenvector, infuenced by the real part of ��. The behavior of �(�) = ����(0) is determined by the nature of the eigenvalues��, as detailed in the box below. Stability analysis in deterministic models using eigenvalues Negative real part (Re(��) < 0): Components ���� decay exponentially to zero as � → ∞. Therefore, any initial perturbation eventually dies out, indicating stability. Positive real part (Re(��) > 0): Components ���� grow exponentially as � → ∞. If any �� has a positive real part, the system will respond to even small perturbations with growing deviations from equilibrium, indicating instability. Purely imaginary eigenvalues (Re(��) = 0): Solution components oscillate, resulting in periodic or quasi-periodic behavior around the equilibrium. Stabil- ity is not guaranteed and requires further investigation, possibly with nonlinear analysis or methods like Lyapunov functions [64]. Zero or repeated eigenvalues: Special attention is required as the presence of zero or repeated eigenvalues with insuffcient geometric multiplicity might necessitate additional investigations into the system’s response to perturbations. These methods, applicable individually or in combination, offer diverse perspec- tives [64]. Some are more apt for linear systems, others for nonlinear ones, and they might provide insights into local or global stability. 2.3 Epidemiological models Epidemiological models aim to describe the spread of infectious diseases within populations and are the standard quantitative tool for understanding transmission dynamics. As specialized forms of population dynamics models, they are similar to ecological models in that they display alternative stable states—such as disease-free and endemic equilibria—but differ in their incorporation of human behavior and intervention strategies. These models capture the roles of key components, including the susceptible, infected, and recovered states, and describe how interactions among individuals, social behavior, population density, and immunity drive the dynamics of disease spread [23]. 27 Moein Khalighi Specifcally, the dynamics of epidemics are infuenced by individual interac- tions, movement patterns, and behavioral adaptations. For example, “super-spreader” compartments can drive the epidemic in a manner akin to keystone species in an ecosystem [65]. In contrast, public health interventions—such as quarantine and social distancing—alter the course of an epidemic by reducing transmission rates and creating conditions that weaken the disease’s progression. Occasionally, an epi- demic may exhibit a threshold effect; at low infection levels, insuffcient host contacts can lead to the elimination of the disease [23]. These observations highlight the importance of incorporating non-medical factors and human behavior when designing epidemic management strategies. Modeling disease spread within populations has garnered signifcant attention in recent years [23]. Many epidemiological studies utilize compartmental models, which segment a population into distinct groups based on their disease status, such as sus- ceptible (�), infected (�), or recovered (�). These models describe the transitions be- tween compartments through differential equations. Notably, the susceptible-infected- recovered (SIR) model, initially introduced by Kermack and McKendrick [66], and its extension, the susceptible-exposed-infected-recovered (SEIR) model, are among the most widely employed compartmental frameworks. In general, a compartmental model can be represented mathematically as: ��� = inputs − outputs, �� where �� is the size of compartment �, and the inputs and outputs refect processes like infection or recovery. Specifcally, the dynamics of the SIR model are governed by the following differ- ential equations: �� �� = −� , �� � �� �� = � − ��, (12) �� � �� = ��, �� where � is the transmission rate and � is the recovery rate, and � is the total popula- tion, satisfying the relation � + � + � = � . Figure 9(a) illustrates the evolution of these compartments over time for a specifc case. The SEIR model refnes this approach by adding an “Exposed” compartment, which accounts for individuals who have been exposed to the disease but are not yet infectious [67; 68]. These models have been extensively modifed to account for the complexities of disease transmission and control [68] and have found applications beyond epidemiology, including modeling the spread of rumors, the dynamics of opinions, the propagation of religion, and the spread of information on social net- works [69]. In Publications III and IV, I utilize two of these extended models, which are described in further detail in Sections 2.3.2 and 2.3.2. 28 Population dynamics modeling 2.3.1 Key quantities Basic reproduction number The study of epidemiology often centers around the basic reproduction number, denoted as ℛ0. This metric measures the contagiousness of an infectious disease by refecting the average number of new infections one infected individual can cause within a susceptible population [70]. A value of ℛ0 > 1 suggests that the disease will likely spread, as seen with several COVID variants, while ℛ0 < 1 indicates a probable decline in transmission. Understanding ℛ0 is vital for assessing the potential scope of an epidemic and devising effective containment strategies. Factors infuencing ℛ0 include the duration of infectiousness, the likelihood of transmission per contact, and the frequency of contacts [70]. These elements can vary, infuenced by the epidemiological triad of agent, host, and environment, which encompasses public health resources, policies, and physical surroundings. For instance, measles has exhibited a wide range of ℛ0 values, affected by local social and environmental conditions, demonstrating how ℛ0 can fuctuate with changes in human behavior and social dynamics [71]. The basic reproduction number ℛ0 is calculated as follows [70]: ℛ0 = Infectious contacts per unit time × Duration of infectiousness Thus, for the dimensionless type of SIR model (12), it yields 1 � ℛ0 = � × = . � � The threshold parameter ℛ0 indicates whether an infectious disease can spread through a population. If ℛ0 < 1, it suggests that the recovery or mortality rate (�) exceeds the transmission rate (�), leading to a decline in new cases until the disease eventually dies out. Conversely, if ℛ0 > 1, the transmission rate surpasses the recovery rate, meaning each infected person spreads the disease to more than one other individual on average, resulting in an outbreak or epidemic as the disease spreads faster than it can be contained. Calculating the basic reproduction number, ℛ0, is not always straightforward due to the complexity of real-world transmission dynamics. Although often treated as a single number, ℛ0 is not unique. Its value depends on the choice of model structure (for example, whether age classes or spatial heterogeneity are included), the assumed distribution of generation times, and the defnition of what constitutes a secondary case. Two models that both ft observed data can yield different ℛ0 estimates, and the same disease can have multiple context-specifc ℛ0 values that refect differences in mixing patterns or intervention coverage [71]. To accurately model these dynamics, it is essential to incorporate additional considerations, parameters, and compartments into the mathematical models, which complicates the estimation of ℛ0. Nevertheless, the next generation matrix method is available to estimate ℛ0 [70; 72]. 29 Moein Khalighi This method involves delineating the compartments that represent the infected states within the model and defning their role in generating new infections. Subse- quently, a transition matrix is constructed that refects the average number of infections each individual causes before they either recover or are removed from the pool of infectious agents. In this method, states are divided into infected and non-infected categories, and two primary matrices are formulated: � , which describes the force of new infections, and � , which describes transitions of new infections. ℛ0 is then determined as the dominant eigenvalue of the matrix product � � −1, which quantifes the average number of secondary infections generated by an infected individual. Final epidemic size The outbreak size of an epidemic (or fnal epidemic size) can be determined by calculating the proportion of the population that remains uninfected once the epi- demic ceases. This calculation is typically encapsulated by an implicit equation that links the epidemic’s fnal size to the basic reproduction number, ℛ0, and the initial conditions [73]. By considering a dimensionless system where � + � + � = 1, we can exclude � � from the analysis. We then derive �� = −1 + , leading to �(�) = −�(�) +�� �� � ln �(�) + �. The constant � can be determined using initial conditions. Setting� �(∞) = 0 allows us to solve for �(∞). Assuming that �(0) is nearly 1, we establish the relationship [73]: �(∞) = �−ℛ0(1−�(∞)), where �(∞) represents the fraction of the population that remains susceptible when the epidemic ends. The proportion of the population that becomes infected is given by 1 − �(∞). If �(0) is not close to 1, slight modifcations are needed, but the relationship remains essentially the same. This relationship illustrates how the disease state shifts from disease-free to endemic as ℛ0 surpasses the critical threshold of 1, see Figure 9(b). These principles are further explored in the subsequent sections. Disease-free and endemic equilibria The disease-free and endemic equilibria are the two common stable states of ODE models describing epidemics. The stability of these equilibria is assessed by evaluating the bifurcation diagram over the parameter ℛ0. The disease-free equilibrium describes a situation where a population has com- pletely eradicated a disease, with no new infections occurring. This state is reached when ℛ0 < 1, indicating that the disease can no longer maintain itself among the population [71]. Achieving this state is a key objective of public health measures, which seek to lower ℛ0 through methods like vaccinations, quarantines, and social 30 Population dynamics modeling (a) SIR Dynamics Time F ra c ti o n o f P o p u la ti o n 1 0 Susceptible Infected Recovered 0 0 (b) Final Epidemic Size R0 O u tb re a k S iz e 1 - S ( ) 1 1 Figure 9. Overview of simple epidemic models and the impact of ℛ0 on fnal epidemic size. (a) SIR model dynamics show disease progression over time by categorizing the population into susceptible (�), infected (�), and recovered (�) groups. The graph represents changes in each compartment over time (e.g., days) based on certain parameters, illustrating a decrease in susceptibles, a peak and decline in infections, and an increase in recoveries. (b) The relationship between the fnal epidemic size, 1 − �(∞), and the basic reproduction number (ℛ0) is shown. This highlights how ℛ0 infuences the fnal infected proportion, with a critical bifurcation at ℛ0 = 1 indicating either disease fade-out or outbreak. distancing practices. Examining the conditions and stability of a disease-free equi- librium is essential for epidemiologists and public health authorities to determine how and when a disease might be effectively eradicated or managed. Additionally, this state provides a framework for assessing the impact of disruptions, such as the introduction of infected individuals or reductions in vaccination coverage, on the population’s health [23]. On the other hand, the endemic equilibrium signifes a state where a disease consistently persists within a population at a stable rate over time [72]. From a phenomenological perspective, this equilibrium is not just about the continuous presence of the disease but also how individuals and communities experience and adapt to it. This involves cultural, social, and economic adjustments that people make to live with the disease as a constant part of their environment. For instance, in areas where a disease is endemic, daily routines, health practices, and community interactions are often shaped around minimizing the risk and managing the impact of the disease, refecting a deeply integrated response that transcends mere biological factors [23]. Backward bifurcation of the disease-free equilibrium The forward bifurcation at ℛ0 = 1, presented above, is a hallmark of the basic SIR framework; however, numerous realistic transmission mechanisms disrupt this sim- plicity. When factors such as partial immunity, vaccination with waning protection, density-dependent contact saturation, exogenous reinfection, or nonlinear treatment effects are introduced, the disease-free equilibrium may undergo a backward (subcriti- 31 Moein Khalighi cal) bifurcation (see, for instance, the tuberculosis model in [74]). In such scenarios, a stable endemic equilibrium can coexist with a stable disease-free equilibrium even when ℛ0 < 1, resulting in bistability. This phenomenon is analogous to the catas- trophic bifurcation previously discussed in Section 2.1.2; thus, Figure 9(b) would resemble the catastrophic fold depicted in Figure 2(c). A small adverse perturbation, such as a temporary lapse in vaccination, could push the system beyond a saddle node and trigger a sudden shift toward high disease prevalence. Consequently, merely reduc- ing ℛ0 below one might not guarantee disease elimination; additional interventions may be required to move the system out of the endemic basin of attraction. 2.3.2 Extended models Modern compartmental models have evolved to incorporate varying infection rates and changing societal dynamics, making them applicable to a wide range of diseases and scenarios. Refning epidemiological models now requires considering numerous factors. For example, in Publication III, I explore the impact of incorporating memory effects alongside dynamic population variables on an Ebola transmission model. In Publication IV, I present a COVID-19 transmission model that includes both environmental pathogens and memory effects, thereby highlighting the importance of sanitation and targeted interventions in infectious disease control. Ebola transmission modeling The management of outbreaks such as Ebola necessitates advanced and reliable mathematical models for disease transmission. Traditional models often struggle to represent dynamic population sizes and long-term epidemiological trends. In Publication III, I address these limitations by developing a model that incorporates variable population dynamics. In our extended epidemic framework, the total population is divided into eight distinct compartments, each representing a specifc stage in the disease progression. This structure builds on traditional SEIR models by incorporating additional compart- ments to consider more complex transmission dynamics and control measures. The total population � (�) is given by: � (�) = �(�) + �(�) + �(�) + �(�) + �(�) + �(�) + �(�) + �(�), where each compartment is defned as follows: Susceptible (�): Individuals at risk of contracting the disease begin in the sus- ceptible compartment �(�). The birth or augmented population rate consists of two components, �1, density-independent birth/augmentation rate, and �2, density- dependent birth/augmentation rate, contributing an increase of (�1 − �1� )� . These individuals become exposed through contact with various infectious groups: contact 32 Population dynamics modeling with actively infected individuals contributes the term �� �� , with hospitalized indi-� viduals the term �ℎ �� , with individuals dead but not yet buried the term �� ��, and � � with asymptomatic infectious individuals the term �� ��. Additionally, natural deaths � occur at a rate of (�1 + �2� )�, where �1 is the density-independent death rate, and �2 is the density-dependent death rate. Thus, the evolution of �(�) is described by: �� �� �ℎ �� �� = (�1 − �2�)� − �� − �� − �� − �� − (�1 + �2� )�. �� � � � � Notice that the coeffcients are scaled by 1/� to account for frequency-dependent contact rates. Exposed (�): Once a susceptible individual is infected, they enter the exposed compartment �(�), where they are infected but not yet infectious. The progression to the infectious stage occurs at rate � (i.e., ��), while natural deaths occur at rate (�1 + �2� )�. Hence, �� �� �ℎ �� �� = �� + �� + �� + �� − �� − (�1 + �2� )�. �� � � � � Infected (�): Exposed individuals become symptomatic and capable of transmitting the disease when they move into the infected compartment �(�) at rate ��. Within this compartment, individuals may transition to an asymptomatic infectious state at rate �1� , be hospitalized at rate �� , or die due to the disease at rate �� . Natural deaths also occur at rate (�1 + �2� )�: �� = �� − (�1 + � + � + �1 + �2� ) �. �� Asymptomatic infectious (�): The compartment �(�) comprises individuals who, although asymptomatic, remain capable of transmitting the disease. They enter �(�) from the infected compartment at rate �1� and from the hospitalized compartment at rate �2� . These individuals eventually recover completely, moving into the com- pletely recovered compartment �(�) at rate �3�, while also experiencing natural deaths at rate (�1 + �2�)�: �� = �1� + �2� − (�3 + �1 + �2� ) �. �� Hospitalized (�): Infected individuals who are isolated and receiving medical care are placed in the hospitalized compartment �(�) at rate �� . Within this compartment, some patients improve and move to the asymptomatic state at rate �2� , while others die and contribute to the burial process at rate �2� . Additionally, natural deaths occur at rate (�1 + �2� )�: �� = �� − (�2 + �2 + �1 + �2� ) �. �� 33 Moein Khalighi Dead but not buried (�): This compartment �(�) represents individuals who have died from the disease but have not yet been buried. They enter �(�) at rate �� and are later processed for burial at rate �1� or removed via incineration at rate ��: �� = �� − (�1 + �)�. �� Buried (�): The buried compartment �(�) includes individuals who have been interred. The infow into �(�) is the sum of those buried from �(�) at rate �1� and those from the hospitalized compartment at rate �2� . These individuals are then further removed from the transmission chain by incineration at rate ��: �� = �1� + �2� − ��. �� Completely recovered (�): Individuals who fully recover from the disease and are assumed to acquire immunity enter the completely recovered compartment �(�) at rate �3�. Natural deaths occur in this group at rate (�1 + �2� )�: �� = �3� − (�1 + �2� )�. �� The following system of differential equations thus governs the overall dynamics of the model: �� �� �ℎ �� �� = (�1 − �2� )� − �� − �� − �� − �� − (�1 + �2� )�, �� � � � � �� �� �ℎ �� �� = �� + �� + �� + �� − �� − (�1 + �2� )�, �� � � � � �� = �� − (�1 + � + � + �1 + �2� ) �, �� �� = �1� + �2� − (�3 + �1 + �2� ) �, (13) �� �� = �� − (�1 + �)�, �� �� = � � − (�2 + �2 + �1 + �2� ) �, �� �� = �1� + �2� − ��, �� �� = �3� − (�1 + �2� )�. �� The model introduces several innovative features compared to the original model [75]. The birth term (�1 − �2� ) represents a density-independent infux of susceptibles, 34 Population dynamics modeling moderated by a density-dependent reduction as the population grows. Additionally, natural deaths are modeled by the term (�1 +�2� ), which departs from the traditional use of fxed per capita death rates. Finally, to include the impact of social behaviors on epidemic dynamics, a new parameter � is introduced to represent the incineration rate. I also assume that �1 > �1 (births exceed natural deaths) and �1 > � (the natural death rate is lower than the incineration rate). Both birth and death rates vary smoothly and monotonically with the population size � . Summing the eight compartmental equations yields the net rate of change of the total population: �� = (�1 − �1)� − (�2 + �2)� 2 − (� − �1 − �2� ) [� + �] . �� The frst two terms describe density-dependent births and natural deaths, while the fnal term accounts for removals via burial and incineration. In this model, deriving the basic reproduction number, ℛ0, requires additional effort compared to the simple SIR framework introduced in Section 2.3.1. It can be obtained using the next generation matrix method [70; 72] as follows: ( ) ( ) � ���1�2�3 + �3 �1�1 + ��2 �� + �1�2��� + �2�3��ℎ ℛ0 = ( )( ) , (14) �1�2�3 �1 + � + � + �5 � + �5 with �1 = �2 + �2 + �1 + �2� * , �2 = �3 + �1 + �2� * , �3 = �1 + �, and �1 − �1 �5 = �1 + �2� * , where � * = . A detailed stability analysis and further �2 + �2 mathematical insights are provided in Publication III. COVID-19 transmission modeling In Publication IV, I introduce an enhanced compartmental model for COVID-19 transmission that builds upon the framework proposed by Mwalili et al. [76]. My model not only incorporates environmental pathogens but also extends to include the dynamics associated with deceased individuals. This model presents the multifaceted nature of COVID-19 spread, including both direct and environmental transmission pathways, while accounting for behavioral adaptations and mortality. The population is divided into seven compartments, each representing a distinct stage of disease progression or environmental viral presence: Susceptible (�): Individuals at risk of infection, who may contract the disease through environmental or direct contact. New susceptibles are continuously added to the population, such as through births (Λ). The infection dynamics involve trans- mission from environmental pathogens (�1Φ1) and direct contact with infectious individuals (�2Φ2). Additionally, if exposed individuals do not develop the disease, 35 Moein Khalighi they may return to the susceptible pool (��). Natural mortality within this group is also accounted for (��). Thus, the dynamics is governed by: �� = Λ − �1Φ1 − �2Φ2 + �� − ��. �� Exposed (�): Individuals infected but not yet contagious, in the latent phase. The transition into this group occurs as susceptible individuals contract the virus, indicated by the terms �1Φ1 and �2Φ2, representing infection via environmental pathogens and direct contact, respectively. The group diminishes as individuals revert to being susceptible (��), succumb to natural mortality (��), or progress to infectious states (��). Thus, its equation is: �� = �1Φ1 + �2Φ2 − �� − �� − ��. �� Asymptomatic infectious (��): Infected individuals capable of transmission with- out exhibiting symptoms. A fraction of the exposed individuals transition to this asymptomatic infectious state, represented by (1 − �)��, where � is the proportion that becomes symptomatic. The population of asymptomatic infectious individu- als decreases due to natural mortality (���) and disease-induced deaths (���), and individuals recover at the rate ����. Thus, the equation is: ��� = (1 − �)�� − (� + �)�� − ����. �� Symptomatic infectious (�� ): Infected individuals displaying symptoms and capa- ble of transmission. The infow of symptomatic infectious individuals is represented by ���, where � indicates the fraction of exposed individuals who develop symp- toms. The decline in this group occurs due to natural mortality (��� ), disease-specifc mortality (��� ), and recovery (�� �� ). Thus, this is modeled as: ��� = ��� − (� + �)�� − �� �� . �� Recovered (�): Individuals who have recovered and are presumed immune. The population of recovered individuals increases through the recovery of both symp- tomatic (�� �� ) and asymptomatic (����) infectious individuals. The only decrease in this group occurs due to natural mortality, denoted by ��. Thus, the dynamics is: �� = �� �� + ���� − ��. �� 36 Population dynamics modeling Deceased (�): Individuals who have succumbed to the disease. The population of deceased individuals increases due to disease-specifc mortality among both symp- tomatic and asymptomatic infectious groups, represented by the term �(�� +��). The outfow from this compartment, modeled by ��, represents the removal processes, such as burial. Thus, its equation is: �� = �(�� + ��) − ��. �� Environmental pathogens (� ): The concentration of virus in the environment that contributes to transmission. Virus particles are shed into the environment by infectious individuals at different rates, depending on whether they are asymptomatic (����) or symptomatic (�� �� ). The natural decay of the virus in the environment is denoted by �� � . Thus, this is described by: �� = ���� + �� �� − �� �. �� The following system of differential equations thus governs the overall dynamics of the model: �� = Λ − �1Φ1 − �2Φ2 + ��(�) − ��(�), �� �� = �1Φ1 + �2Φ2 − ��(�) − ��(�) − ��(�), �� ��� = (1 − �)��(�) − (� + �)��(�) − ����(�), �� ��� = ���(�) − (� + �)�� (�) − �� �� (�), (15) �� �� = �� �� (�) + ����(�) − ��(�), �� �� = �(�� (�) + ��(�)) − ��(�), �� �� = ����(�) + �� �� (�) − �� � (�). �� �(�)� (�) �(�) (��(�) + �� (�))In which Φ1 = and Φ2 = model the infec- 1 + �1� (�) 1 + �2 (��(�) + �� (�)) tion incidence functions with saturation effects, refecting adaptive behaviors like social distancing. Similar to the previous section, I employ the next generation matrix method [70; 72] to determine the basic reproduction number (ℛ0) of this model, as outlined below: ( ( ) ( )) Λ� � (1 − �) �� � ��(1 − �)ℛ0 = ( ) �2 + + �1 + , (16) � ����� ����� �� ����� �� ����� 37 Moein Khalighi where �� = � + � + �, ��� = � + � + ��, ��� = � + � + �� . Detailed stability analysis and further mathematical insights are provided in Publica- tion IV. 38 3 Memory effects In this chapter, I delve into the systematic analysis of memory effects in order to address the frst research question about incorporating memory. I introduce the concept of memory effects in dynamical systems. This concept shapes the basis of this dissertation and the publications. Finally, I discuss implementing memory effects and other relevant approaches and theories within the domain. 3.1 Basic concepts 3.1.1 Memoryless and memory-integrated dynamics Memoryless dynamics A system is memoryless if its future dynamics/evolution is determined solely by its current state—no additional information about how that state was reached is required. In other words, all the necessary information to predict the system’s future is contained in its present state, and the past does not infuence the outcomes or the destinations of dynamics (Figure 10a). This is analogous to the “Markov property” in stochastic processes [77]; however, since I focus on deterministic models, I will avoid stochastic terminology and refer to the memoryless property. A discrete example is a simple thermostat control system without hysteresis. In such a system, the heater turns on or off based solely on the current temperature reading. For instance, if the temperature drops below a threshold, the heater turns on; if it rises above another threshold, it turns off. This illustrates an extreme case of the memoryless property, where no history is ever taken into account. In continuous examples, many standard population dynamics systems exhibit an analogous but slightly different type of memoryless property. For instance, exponential population growth or decay depends solely on the current population size [61]; the rate of change at any given moment is determined by the present population, with no dependence on how the system arrived at that state. Similarly, early-stage bacterial [78] or yeast growth in a nutrient-rich environ- ment [79], or an introduced species thriving in a new habitat with abundant resources and few predators [80], can each approximate exponential behavior under ideal con- ditions. In these scenarios, the underlying differential equation depends solely on the present population (or substance amount) to determine the rate of change rather than tracking the entire past trajectory. Although real-world factors—such as resource 39 Moein Khalighi depletion or competition—eventually break pure exponential trends, such effects can often be incorporated using alternative memoryless models (e.g., logistic growth with carrying capacity). Under simplifed or short-term conditions, however, the exponential model suffces for the system’s memoryless nature. without memory with memory P o p u la ti o n d y n a m ic s P o p u la ti o n d y n a m ic s P o p u la ti o n d y n a m ic s P o p u la ti o n d y n a m ic s Past Present FuturePast Present Future Past Present FuturePast Present Future Time series (a) (b) Figure 10. Conceptual illustration of memory effects in dynamical systems. (a) In a memoryless system, past states do not infuence future dynamics. Therefore, identical current or recent conditions always lead to the same future outcomes, regardless of historical differences. (b) In a system with memory, historical states affect future outcomes. As a result, even if current or recent states are identical, variations in past conditions can lead to different future scenarios. Dynamics with memory In contrast to memoryless dynamics, the dependency on historical data makes dynam- ics with memory more complex and often leads to more realistic models of physical, biological, or social systems, where history plays a signifcant role in determining fu- ture outcomes [7; 81; 82]. Figure 10b conceptually illustrates cases in which historical factors shape future dynamics. For example, in a forest ecosystem recovery after a fre, the future state of the forest ecosystem after a wildfre depends signifcantly on its past states [7]. The speed and manner of recovery are infuenced by the current conditions (like weather or soil fertility) and the forest’s historical characteristics. These include the severity and frequency of past fres, the dominant vegetation types before the fre, and the historical interactions between different species in the forest [7]. The memory of the ecosystem, in terms of its resilience and response to disturbances, is a critical factor in its recovery process. As another example, microbial communities can illustrate systems with memory 40 Memory effects effects. For instance, antibiotic use impacts gut microbial communities, affecting not only illness-causing bacteria but also the diverse gut microbes [81]. These drugs induce long-term changes in microbial composition, resistance patterns, and functional abilities [83]. Repeated antibiotic treatments infuence the microbiome’s future response, fostering antibiotic-resistant bacteria and reducing benefcial microbes. This creates memory effects in the microbial ecosystem, where its current state refects past antibiotic exposures [84; 85]. In epidemiological studies, the dynamics of disease spread within a population previously exposed to the same or similar pathogens are infuenced by memory effects [82]. These effects arise from the collective immunity developed through past exposures, vaccinations, and recoveries from earlier outbreaks [86; 87; 88]. Such historical interactions with the disease shape the population’s immune response and impact future transmission patterns [89; 90]. Consequently, the current state of disease spread is determined not only by immediate factors, such as present interactions and transmission rates, but also by this accumulated immune memory. The presence of this memory introduces complexity into the system’s dynamics, transitioning it from a straightforward model to a more intricate one where past states infuence future outcomes. This characteristic accounts for the non-linear interplay between the disease’s current spread and the population’s historical relationship [22]. 3.1.2 Embedding memory in models Modeling memoryless dynamics Consider a memoryless time series �0, �1, . . . , �� representing population density, where �0 is the initial state. The differential equation �� = � (�, �) �� describes the system’s evolution, where � (�, �) represents the instantaneous rate of change of � with respect to time, depending on both time � and the current state �(�). This equation can be approximated using the Euler method [61], which updates each subsequent state ��+1 as ��+1 = �� + ℎ� (��, ��), where ℎ is a small time step. This approach illustrates that the model only depends on the current state �� to determine the next state ��+1, with no input from earlier states. As ℎ becomes very small, recursively applying Euler’s method over many small intervals makes this discrete update approach a continuous change. This transition leads to the integral form, ∫ �� �� = �0 + � (�, �(� )) ��, �0 41 Moein Khalighi where the sum of changes over these intervals closely approximates the integral over the entire period. The integral also emphasizes that the model’s output is solely dependent on the initial state �0 and the integral of the rate function � , without any memory of previous states. This characteristic defnes the memoryless nature of these models. Modeling dynamics with memory The aforementioned standard differential equations are insuffcient for capturing the dynamics of history-dependent systems. They lack a mechanism—specifcally, a memory kernel—a function that weights the infuence of past states over time. This limitation necessitates a revised equation that incorporates a formula for the current state based on past states: ∫ �� �� = �0 + �(�� − �)� (�, �(�)) ��, (17) �0 where �(�� − � ) is the memory kernel that accounts for past infuences, and �0 and �� are the starting and current times, respectively. In a scenario devoid of memory effects, one might conceptualize the kernel as �(�� − �) = 1, suggesting uniform infuence over time. More precisely, the integral form of a standard differential equation suggests an immediate, unweighted dependence on prior states. Yet, the memory-based equation could simplify to resemble the memoryless scenario if � were a Dirac delta function (�(�� − �) = �(�� − �) [16; 17]), which theoretically aligns with the differential approach but fails to represent the impact of historical states. Figure 11a illustrates the convolution of the equation function � and kernel �. Each case study requires a tailored determination of this kernel function to weigh the relevance of past events and better model the observed phenomena. Detecting memory kernels from data is a complex task; however, evaluating the impact of memory is comparatively more straightforward. As a result, researchers explore different types of memory kernels to understand how various memory effects infuence the behavior of dynamical systems [91]. To classify memory kernels, Ogle et al. [12] identifed three main components of memory: (1) memory length, which measures how far into the past conditions or states — whether short- or long-term — infuence current processes or states; (2) temporal pattern, which refects how the importance of past conditions varies over time, including the presence of signifcant time lags; and (3) memory strength, which quantifes the magnitude to which past conditions impact. Figure 11b shows examples of different memory kernels. These can range from those with short or long memory lengths to those with large or small time lags. For example, power-law kernels are used to model systems with long-range memory [92], where even distant past states continue to shape future dynamics. Exponential decay 42 Memory effects time time current time initial time Memory Kernel Memoryless Dynamics = = time Dynamics with Memory current time Convolution (a) (b) Power Law Exponential Small Lag Time-varying or Stochastics Big Lag Modeling Memory Eects Memory Kernels Figure 11. (a) This panel demonstrates the convolution of an equation function � with a sample kernel �, illustrating memory effects in a dynamical system �, modifed by the integral kernel �. (b) This panel presents a variety of memory kernels spanning a range of memory lengths from short to long and incorporating various time delays. Power-law kernels, ideal for long-range memory systems, demonstrate how past states substantially impact the dynamics of future systems. In contrast, exponential decay kernels are typical of short-term memory systems, where the effects of prior events quickly fade away. Two memory kernels with small and big lag types are useful where specifc past moments signifcantly infuence current states. The last is stochastic kernels that model systems with random or variable memory patterns. kernels often model systems with short-term memory where the infuence of past events exponentially diminishes over time. Moreover, memory effects in scientifc modeling can be represented by intro- ducing time delays into differential equation models. Different delay types–such as discrete delays (a fxed fnite lag) and distributed delays (a continuous spread of past infuences)–explicitly incorporate memory by allowing prior states to affect current dynamics [93]. This dependence on history can alter a system’s behavior; for instance, large delays are often destabilizing, whereas shorter delays can enhance stability [94]. Because memory-integrated models introduce additional parameters (such as the delay length or a memory kernel), parameter estimation techniques are required to infer these values from empirical data, thereby quantifying how strongly and over what timescale past states infuence the system. Stochastic memory kernels (i.e., time-varying memory functions) are employed to model systems that exhibit inherently random memory effects, thereby modeling unpredictable infuences as seen in fnancial markets, turbulent physical processes, or protein reaction kinetics [95]. In this thesis, I focus on the power-law kernel, valued for its slow decay rate that operates without lag effects, allowing early system states to maintain a persistent, though gradually weakening, infuence on the present. This property enables the kernel to capture long-range dependencies, mirroring the scale-invariant behaviors observed in many complex systems [96; 97; 98]. Unlike rapidly decaying kernels, 43 Moein Khalighi such as exponential ones, which may impose artifcial time scales and diminish the continuous impact of history [92], the power-law kernel preserves a realistic, unbroken connection to past states. Standardizing the memory pattern across studies assigns a consistent value to memory strength, facilitating comparisons across different systems. Furthermore, this alignment with natural scaling phenomena not only enhances the model’s realism [99] but also streamlines the mathematical framework [100] (see Section 3.3.2), establishing a robust foundation for deeper investigations in subsequent sections. 3.2 Memory in interacting systems For community-based scenarios, memory effects become especially complicated in multi-dimensional systems modeling (see Section 2.2.2). In these systems, interactions among components—beyond memory effects—add further layers of complexity. Consequently, this thesis focuses on a particular type of memory, and for alternative approaches, I refer readers to Sections 3.4 and 3.5. 3.2.1 Individual component memory In multivariate systems, I focus on specifc types of memory effects in which each component’s past states infuence not only its own current state but also the behavior of other components through their interactions, as illustrated in Figure 12a-b. In this case, the state of each component in a three-dimensional model, such as � , � , and �, can be systematically described as follows: ∫ �� �� = �0 + �1(�� − �)� (�, �, �), ��, ∫ �� 0 � �� = �0 + �2(�� − �)�(�, �, �), ��, (18) ∫ �0 �� �� = �0 + �3(�� − � )�(�, �, �), ��, �0 where �� (for � = 1, 2, 3) represents the memory kernels for each component. The functions � , �, and � in the equations model the interactions among the components � , � , and �, respectively. 3.2.2 Coherent and varied memory In a multi-dimensional system, I classify these memory effects into two categories based on the components’ memory strength levels, whether uniform or varied. When all components exhibit identical memory strengths (�� = � for all �), I term this a coherent memory, depicted in Figure 12c. On the other hand, a varied memory is 44 Memory effects Figure 12. Illustration of memory effects in a tri-variate system with interacting components. (a) Dynamical systems: Each component is represented by a distinct color. The colors fade into transparency towards the past, symbolizing the diminishing infuence of memory. This illustrates the principle that the past states of each component infuence the current state of the system. (b) Constant interaction: A toy model of three components linked by steady interactions (gray connections). Each component’s memory not only infuences its own current state but also affects other components through these interconnections. (c) Coherent memory: This scenario occurs when all components display uniform memory strength levels, characterized by an identical power-law memory kernel. (d) Varied memory: This is observed when at least one component exhibits a distinct memory strength, different from the others. observed when at least one component shows a different memory strength (�� ̸= �� for some � ̸= �), as seen in Figure 12d. Systems with coherent memory are generally easier to analyze. For instance, when all components share the same memory strength, their behavior tends to be uniform in response, simplifying the governing equations and making them easier to solve [17]. Analytical techniques, such as Laplace or Fourier transforms, can be more easily applied under such conditions because each term in the equation transforms similarly [16; 101]. This simplifcation benefts the algebraic manipulation in the transformed domain and the numerical computation, making the latter more effcient. In real-world applications, coherent memory systems are often easier to interpret and ft empirical data because they present fewer variables in the optimization problem. Stability analysis methods are also better established for systems with coherent memory, facilitating a more straightforward assessment of the system’s qualitative behavior through phase-plane analysis [102] (see Section 4.2.1). However, this ease of analysis comes at the cost of losing the rich behavioral diversity that varied memory systems can offer. While coherent memory systems are more straightforward, varied memory systems provide a broader spectrum of behaviors [103; 104; 105]. Varied memory adds fexibility, allowing for a better ft to empirical data and accommodating a wider range of multi-scale dynamics [106]. 45 Moein Khalighi Although this increases the complexity of analytical and numerical evaluations, it provides a more detailed and versatile framework for representing complex real-world systems. 3.3 Modeling memory with fractional calculus This section explores how fractional calculus incorporates memory effects into dy- namical systems through non-local operators, particularly fractional derivatives. Un- like classical integer-order derivatives—which depend solely on the instantaneous state—fractional derivatives involve a convolution integral that weights past states via a kernel. This framework provides a well-defned and systematic way to include mem- ory dynamics in differential equations [18], helping to analyze and measure memory effects in models, especially those of complex systems with interacting components. Here, I demonstrate fractional calculus in representing and implementing memory effects. 3.3.1 Fractional calculus Fractional calculus, distinguished by its focus on derivatives and integrals of non- integer orders, has a rich history in mathematics, dating back to the 17th century. It began with Leibniz’s idea of half-ordered differentiation in his correspondence with L’Hôpital [107]. Over the centuries, notable mathematicians like Abel, Riemann, and Liouville have advanced this feld, establishing its theoretical foundations [17; 16]. Nonetheless, for all its roots in pure mathematics, it is only in the contemporary era that fractional calculus found resonance in practical applications [19]. The value of fractional derivatives lies in their superior capacity to describe com- plex systems compared to traditional derivatives [100], as they account for memory and hereditary properties inherent in many real-world phenomena. This offers insights into systems’ overall behavior and history, extending their utility to diverse felds. For instance, they are employed to comprehend the characteristics of materials with both fuid and solid properties [108] and to analyze the complexities of fnancial markets [109]. A key feature of fractional derivatives is their ability to incorporate memory effects into traditional integer-order systems [22; 110]. To contextualize our study within the broader scholarly landscape, I reference the thorough exploration of fractional derivatives and integrals presented in [111; 17; 16] in various scientifc and engineering domains. These emphasize the signifcance of fractional calculus in practical situations and connect theoretical principles to their practical implementations. Furthermore, the perspectives offered in the review found in [19] shed light on the increasing applications of fractional calculus across various disciplines. This review provides a wealth of information for both novice and seasoned researchers. 46 Memory effects 3.3.2 Caputo derivative Unlike the unique defnition typically used for integer-order derivatives, fractional derivatives can be defned in multiple ways. This multiplicity stems from the richer mathematical framework of fractional calculus, which offers various formulations to express the memory effects observed in different contexts and systems [111]. Among these, the differential operator proposed by Caputo stands out for its widespread acceptance and usage within the scientifc community [16; 17], largely due to its ad- vantages in handling non-zero initial conditions and decaying memory effects—points I detail later in this section. Today, the Caputo derivative is recognized as one of the most widely used fractional differential operators. Consider the notation �� for the time-fractional Caputo derivative of order�0 � ∈ (0, 1], starting from the initial value �0. For a function � that is differentiable, its fractional derivative at time � according to Caputo is: ∫ � 1 � ′ (�) �� 1−� �� �(�) = � � ′ (�) = . (19)�0 �0 Γ(1 − �) �0 (� − �)� 1−� Here, ��0 denotes the Riemann-Liouville fractional integral of order 1 − �, expressed as: ∫ � 1 �(� ) �� �� � 0 �(�) = , � ∈ R+ , (20) Γ(�) �0 (� − �)1−� where Γ is the gamma function. In Equation (19), the integral accumulates contributions from the past derivative � ′ (�), weighted by (� − �)−� , where the gamma function normalizes the kernel. This integral emphasizes the non-local nature of the fractional operator and embodies the memory effects by summing over all previous times � . Multivariate models utilizing fractional derivatives, as follows: ��� �0 �� = ��(�, X), (21) � ∈ N>1 , �� ∈ (0, 1], X = {��}, implicitly introduce a memory correlation function, which creates a dependency between the current state of a system and its historical trajectory through a convolution integral. Alternatively, this relationship can be expressed using an equivalent Volterra integral equation as [17]: ∫ �� ��(��) = ��(�0) + ��� (� − �) ��(�, ��(� )) ��, (22) �0 where the memory kernel is defned by (� − �)��−1 ��� (� − �) = . Γ(��) 47 Moein Khalighi Unlike classical derivatives, the fractional derivative incorporates a history of the function. The convolution with the power-law kernel ��� (� − �) means that every past state � contributes to the derivative at time �. The parameter �� modulates the strength of this memory effect—smaller values correspond to a longer-lasting memory, whereas �� = 1 yields the memoryless (local) derivative. Therefore, in my thesis and associated publications, I adopt 1 − �� as a measure of memory strength for the �th component, quantifed on a scale from 0 to 1. Pub- lication II demonstrates the infuence of a range of memory strength levels on the convergence of the logistic growth model. According to Ogle’s classifcation [12] (see Section 3.1.2), these kernels exhibit long-term memory with no lag effects but similar decaying patterns, where the strength of memory and decay can be adjusted based on the order of derivatives. The Caputo derivative is characterized by its distinct advantages [16; 17; 112]. Firstly, it is adept at modeling systems with non-zero initial values, as the derivative of a constant is inherently zero. It also features a singular power-law kernel, effec- tively representing systems with decaying memory effects. A parameter � > 0 exists such that lim�→∞ �−��� (� − �) is fnite for fxed � , highlighting its applicability in real-world scenarios [113; 100]. Incorporating power-law distributions into mathe- matical models offers advantages such as self-similarity and straightforward Laplace transforms, which enhance the analytical tractability of complex systems [17]. Despite some theoretical criticisms [114], the empirical robustness of power-law distributions and the mathematical properties they offer, such as semi-group properties and simple mathematical transformations, justify their widespread use [17]. These properties make them suitable for theoretical applications, including those employing Caputo derivatives to model complex phenomena with memory effects. This robust empirical foundation and mathematical fexibility form the basis of their application in my thesis to model memory effects, asserting the suitability of Caputo fractional derivatives in systems governed by power-law functions. Commensurate and incommensurate fractional orders In the framework, I distinguish between general memory terms and technical fractional- order terms to describe memory effects in the system. I will use coherent memory and varied memory when discussing memory effects in general (see Section 3.2.2), in qualitative terms. Referring to the mathematical formulation involving fractional derivatives, I introduce the terms commensurate fractional order and incommensurate fractional order. A commensurate fractional order system corresponds to a coherent memory system, in which Caputo fractional derivatives of a uniform order are ap- plied across all variables (�� = � for all �). The uniformity of the derivative order means that all components of the system share the same strength of memory effect, which maintains consistency and aids in analytical tractability. Conversely, an in- commensurate fractional order system corresponds to a varied memory system and 48 Memory effects uses Caputo fractional derivatives of different orders for different variables (�� ̸= �� for some � ̸= �). In an incommensurate system, different parts of the model exhibit different memory strength, enabling the representation of more complex behaviors and responses. This distinction in terminology is important to clarify, as commensurate and in- commensurate fractional orders are technical concepts frequently used in the following sections. 3.3.3 Stability regions In this section, I explore one of the key characteristics of fractional order systems as a foundational analysis for understanding how memory integration enriches the theoretical framework of complex systems. Consider the linear autonomous system: ��� 0 �(�) = ��(�), (23) where �(�) ∈ R� is an �-dimensional state vector and � is an � × � matrix with real entries. The operator �� denotes the Caputo fractional derivative of order �, with�0 0 < � < 1. For the system in Equation (23) to be asymptotically stable (see Sections 2.1.1 and 2.2.3), every eigenvalue � of � (collectively referred to as the spectrum, spec(�)) must satisfy the following condition [115]: �� |arg(�)| > for all � ∈ spec(�). 2 In other words, the phase angle of each eigenvalue must lie outside the cone defned by ±�� in the complex plane. Notice that when � = 1 (the integer-order case), this 2 condition becomes � |arg(�)| > , 2 which is equivalent to requiring that Re(�) < 0, aligning with the stability criteria discussed in Section 2.2.3. An interesting implication of this criterion is that fractional-order systems (0 < � < 1) tend to have larger stability regions compared to their integer-order coun- terparts [116]. Figure 13 illustrates the stability region for a system governed by the Caputo fractional derivative: The unstable region is confned within the cone bounded by lines at angles ±�� in the complex plane. Hence, for � < 1, the stable 2 region exceeds the conventional half-plane (Re(�) < 0) associated with integer-order systems. However, it is essential to consider the practical applicability of these observations. The scenario described is somewhat idealized, primarily focusing on homogeneous systems with coherent and specifc memory. For example, various memory types correlate with distinct stability regions, as highlighted in [116]. Importantly, in 49 Moein Khalighi Figure 13. Stability region for a system with the Caputo fractional derivative (0 < � ≤ 1). The unstable region lies within the cone defned by ± �� in the complex plane. 2 reality, ecological and social systems seldom remain in a relaxed phase devoid of disturbances for extended periods. Thus, the conclusions drawn from such ideal settings about the impact of memory on the stability of complex systems may not be entirely generalizable. For other approaches to stability analysis in ecological systems, I refer the reader to Sections 4.1.1 and 4.1.2, as well as Publication II. Additional theoretical analysis for epidemiological systems can be found in Publications III and IV. 3.3.4 Implementation Fractional calculus presents challenges that often exceed the scope of traditional analytical methods. In particular, the numerical implementation of fractional differen- tial equations (FDEs) is considerably more complex than that of their integer-order counterparts—a fact that can pose signifcant hurdles for practitioners without an extensive background in advanced mathematics [19]. To bridge this gap and catalyze broader community engagement, I developed FdeSolver, an open-source, state-of-the- art package implemented in Julia. As detailed in Publication I, FdeSolver is designed not only to provide effcient numerical solutions for FDEs but also to be accessible and adaptable for real-world applications. At its core, FdeSolver transforms FDEs into solvable Volterra integral equations. For instance, an � size vector form of a fractional-order differential equation �� �(�) = � (�, �(�))�0 is reformulated as ∫ � �(�) = ��−1[�; �0] + 1 (� − � )�−1� (�, �(�)) ��, (24) Γ(�) �0 50 Memory effects where ��−1[�; �0] denotes the Taylor polynomial of degree � − 1 for �(�) centered at �0, and � is the smallest integer that is greater than or equal to the derivative orders �. This approach leverages foundational numerical algorithms from [101; 117; 118] and builds on methods previously implemented in MATLAB routines [119; 120]. To address the non-smooth behavior inherent in fractional differential equations, FdeSolver employs product-integration rules alongside predictor-corrector methods and the Adams multi-step process [121]. In the case of infexible problems, the package incorporates an implicit method by frst restructuring the discrete formulation as �� = Ψ�−1 + �0��, (25) followed by an iterative Newton-Raphson procedure: [ ]−1 ( ) � [�+1] = � [�] � − �0�� (� [0] � [�]− ) − Ψ�−1 − �0�� , (26)� � � � with �� (��) representing the Jacobian matrix of � , Ψ�−1 denoting the term in- cluding all the explicitly known information of the previous steps, and [�] indicating the number of iterations. Recognizing that direct matrix computations can be im- practical for large time steps, FdeSolver integrates Fast Fourier Transform (FFT) techniques [122] to reduce computational costs and enhance effciency in large-scale simulations. FdeSolver is engineered to handle one-dimensional equations with orders given by positive real numbers while also extending its capabilities to high-dimensional systems with orders up to one and incommensurate derivatives. Although the examples presented in Publication I predominantly employ zero initial conditions, the package fully supports nonzero initial conditions, making it versatile for a wide range of practical scenarios. Benchmark comparisons indicate that FdeSolver outperforms existing MATLAB routines and Julia’s FractionalDiffEq package (v0.3.1)—in terms of both accuracy and computational effciency. A key performance advantage stems from Julia’s just-in-time (JIT) compilation and its multiple dispatch paradigm [123], which enable the creation of generic, ef- fcient, and elegant functions. Furthermore, the modern ecosystem provided by Julia—with its intuitive package manager and interactive Read-Eval-Print Loop (REPL)—enhances usability. FdeSolver’s open-source nature, comprehensive doc- umentation, and rigorous testing encourage community contributions and continual refnement. The package also boasts seamless compatibility with languages such as C, Fortran, R, and Python, facilitating integration into diverse computational projects. Publication I further demonstrates FdeSolver’s capabilities through fve key exam- ples, including applications in modeling microbial dynamics and epidemic disease spread, thereby reinforcing its practical relevance. Looking ahead, planned develop- ments include the incorporation of fractional linear multistep methods [124], support for models with time-varying derivatives [95], automated cross-validation to prevent 51 Moein Khalighi overftting, advanced parameter optimization, Bayesian approximation methods via the Turing.jl package [125], and GPU support [126] to accelerate computations. More technical details and case studies demonstrating FdeSolver’s practical appli- cations are elaborated in Publication I. 3.4 Alternative concepts In this thesis, I focus on the infuence of past population states on system dynam- ics. However, it is also possible to consider the effects of past perturbations or past interactions. For example, in the context of antibiotic resistance (explained in Sec- tion 3.1.1), the evolution of species, specifcally, the transition from nonresistant to resistant strains, can be infuenced by the history of antibiotic interventions. In this sense, memory relates to the emergence of resistance, suggesting that memory and resistance dynamics are interconnected or even represent different facets of the same phenomenon. This raises another consideration: in what way past perturbations, rather than past population sizes alone, shape current population dynamics. One particularly relevant concept is interaction memory, studying the impact of past interactions. Traditional models often assume constant interactions between com- ponents, an assumption that may not hold true in evolving systems, particularly those where past interaction strengths continue to shape current dynamics, whether popula- tion dynamics or interaction dynamics. For example, in ecological contexts, species interactions may weaken or strengthen over time in response to resource fuctuations or past conficts [127; 48]. Social networks also reveal that repeated collaborations or rivalries can accumulate over time, shaping ongoing relationships [128]. For example, Williams et al. [129] show how patterns of direct and extended reciprocity in repeated interactions can lead to diverse social norms, illustrating how a history of cooperative or adversarial behavior molds current group dynamics. time S tr e n g th o f in te ra c ti o n s Figure 14. Interaction memory in a dynamically interactive model. This model represents the concept of memory effects as the infuence of each interaction’s past on the current states of the interactions. The left plot presents a sample of a tri-variate model’s dynamic interactions, illustrating the components’ interplay. The color coding of each element signifes the specifc components involved in these interactions. The interactions are treated as symmetric (right plot), leading to just three distinct interactions in this three-dimensional framework, simplifying the representation of the interconnected system. To account for such histories, one could consider dynamic interactions, as illus- 52 Memory effects trated in Figure 14. Enhancing models to include the effects of past interactions on current values can be insightful in contexts where path dependence or hysteresis are signifcant drivers of the system’s behavior. As with component memory, interaction memory can be formalized using a matrix representation—its diagonal elements indicate the self-memory of each interaction. While incorporating interaction memory does increase model complexity, the potential gain in realism and predictive power can justify the added detail in many applications. A basic version of this model for binary interactions is proposed in [129]. Other innovative models offer diverse aspects of memory, such as delay effects, time-varying memory, and various memory kernels [94; 91; 95]. These models potentially suggest exploring various approaches, such as varying memory over time or across components and considering constant or dynamic interaction between components. While these models fall outside the primary focus of this dissertation, they could contribute to our understanding of memory in complex systems, covering a broad spectrum of phenomena and stages of system evolution. 3.5 Alternative approaches The selection of a modeling approach is linked to the specifc requirements of the task, be it the need for empirical robustness, mathematical precision, or adaptability to the complexities and high-dimensional nature of the data or case study. Beyond fractional derivatives, a variety of complementary methodologies have been developed to incorporate memory and handle complex dynamics in different contexts. For instance, the Mori-Zwanzig formalism [130], rooted in statistical physics and stochastic differential equations, addresses the challenge of high-dimensional systems with limited observable variables. By positing that the system’s future behavior can be inferred from the past data of a smaller set of variables, it aligns with Takens’ Theorem [131], which asserts that the history of a single aspect can capture the overall system behavior. While this framework and Takens’ theorem enable analyzing complex systems from limited data, they require careful consideration of delay times, embedding dimensions, and the separation of relevant from irrelevant components [132]. Time-series approaches offer alternative ways to model memory. ARIMA mod- els [133] address short-term dependencies by differentiating data for stationarity and then applying autoregressive and moving average components. In contrast, ARFIMA [134] extends this idea with fractional differencing, thereby accommo- dating longer-range dependencies similar to those found in fractional derivatives. Another approach, the Sparse Identifcation of Nonlinear Autoregressive (SINAR) methodology [135], merges sparse identifcation of nonlinear dynamics [136] with autoregressive modeling to predict behavior from aggregated data (e.g., opinion dy- namics in agent networks). Although robust in handling nonlinearities, SINAR’s memory is encoded through discrete autoregressive terms rather than continuously. 53 Moein Khalighi Statistical methods like Random Forest models can manage high-dimensional data when inferring memory, but they may struggle with overftting and reduced interpretability [14]. In contrast, stochastic models, which emphasize probabilistic transitions, handle uncertainty effectively but can be complex to specify [12]. Markov chain refnements have likewise advanced our understanding of memory-inclusive network dynamics [137; 138], with applications in community detection, air traffc analysis, and scientifc collaboration networks [128; 139]. In particular, incorporating memory into network models can alter spreading patterns and transmission times, applicable to describing phenomena such as disease spread [140; 141]. Moreover, combining machine learning with stochastic models extends the reach of these techniques into areas as diverse as particle physics and fnancial markets [142], highlighting the versatility of hybrid approaches. Recent explorations of fractional- order artifcial neural networks further demonstrate the broad applicability of fractional calculus in various modeling frameworks [143]. In general, these developments demonstrate how different methods—ranging from statistical and stochastic models to network-based analyses—can complement fractional derivatives, offering a rich toolbox for inferring and describing memory and complexity across a spectrum of real-world systems. 54 4 Impact of memory on ecological systems In this chapter, I explore the impact of memory on the behavior of ecological systems under various conditions and perturbations, thus, addressing the second research question. This involves exploring the infuence of coherent and varied memory, using commensurate and incommensurate fractional derivatives (see Section 3.3.2), on the stability, resilience, and resistance of systems. While Section 3.3.3 introduced the role of fractional derivatives in defning the stability regions of idealized models, this chapter extends that discussion with analytical and numerical analyses of both univariate and multivariate systems. 4.1 Univariate systems In this section, I focus on examining the effects of memory on bi-stable, one- dimensional systems, typically characterized by a potential function. Simple de- sertifcation models, briefy introduced in Section 2.2.1, exemplify such systems. Here, I present how memory infuences the studied systems’ stability landscape, resilience, resistance, and other critical attributes. Consider the Herbivory model (3) introduced in Section 2.2.1, which features two stable states, ��1 and ��2, separated by one unstable state, �� . Figure 15a illustrates that memory slows down the system’s dynamics. Notably, the inclusion of memory does not alter the outcomes in univariate systems; the equilibrium points remain the same in the relaxed phase, as shown (Figure 15b). This occurs because the roots of the derivative function are identical in systems with and without memory. 4.1.1 Memory and stability landscapes To better understand the impact of memory on dynamics, I analyze its effect on the stability landscape. As discussed in Section 2.1.7, the stability landscape is closely linked to a system’s equilibrium points. Building on the formula (5), to evaluate the stability landscape of a system with memory, we need to solve the following equation: ∫ �� ˙� = − ���, (27) �� 55 Moein Khalighi Figure 15. The dynamics of the Herbivory model (3) and its bifurcation diagram versus parameter �, the herbivore population density. (a) The model’s dynamics for three different initial values converge on two equilibrium points. Blue lines represent dynamics without memory, while red lines indicate dynamics with memory. (b) The bifurcation diagram shows the equilibrium points, which remain unchanged irrespective of memory effects. ��1 and ��2 denote stable states, and �� indicates an unstable state. in which �˙ denotes frst derivative of �(�) with respect to time, and the state �(�) depends on the history of � via a memory kernel �, i.e., on past values � (�, �(�)) weighted by �(� − �): ∫ �� �� = �0 + �(�� − �)� (�, �(�)) ��, �� ≤ �� ≤ ��. (28) �0 Here, �� and �� represent the boundaries that defne the range of states within the stability landscape, specifcally chosen to include both the valleys and the peak in a bistable system. Since the roots of the derivative functions remain unchanged by memory, the positions of the peaks and valleys in the landscape stay the same in the absence of perturbations (see Section 2.2.3 and Figure 7). In other words, the distance to the basin threshold—the gap between a stable state and the edge of its basin of attraction—remains constant regardless of whether the system has memory. However, the results indicate that memory affects other aspects of the stability landscape, such as the shapes of basins of attraction, their potential depths, and curvatures (see Figure 16a). This occurs because the trajectories by which states converge to equilibrium differ when memory is present (Figure 15a). Figure 16b further shows that the stable points (��1 and ��2) and the unstable point (�� ) correspond to the roots of the potential slope curve (refer to Section 2.2.3). The curvature of the potential, as indicated by the slopes (blue and red lines), refects the rate of change at the equilibrium points, with shallower slopes indicating slower dynamics. This deceleration effect, induced by memory, is evident in the altered curvature of the potential landscape. 56 Impact of memory on ecological systems Figure 16. (a) Stability landscape (potential energy) corresponding to the dynamics in Figure 15a. (b) The slope of the potential energy indicates the rate of change in dynamics and, consequently, the speed of transitions between states. ��1 and ��2 denote stable states, and �� indicates an unstable state. The analysis of Equation (27) poses computational challenges because each state � is infuenced by an evolving integral that refects the impact of initial conditions. This results in a dynamic landscape where trajectories vary based on their starting points. Consequently, estimating the stability landscape in systems with memory requires either accounting for or fxing the initial conditions. Future studies could further investigate how memory infuences the depth and curvature of the stability landscape across various univariate models. Such research would provide deeper insight into how memory modulates resilience and resistance, thereby affecting the system’s response to perturbations. 4.1.2 Stability landscape transformations Understanding how a system’s stability landscape responds to perturbations is crucial, given that natural systems are rarely in a state of equilibrium and are continually subjected to environmental disturbances. Here, I extend the analysis to consider the dynamic effects of memory during perturbations. In Section 2.1.7, I demonstrated that stability landscapes change with varying conditions. Our focus is on abrupt changes—such as sudden shifts in light, tempera- ture, chemical concentrations, or biotic factors—rather than gradual transitions. To model these changes, I introduce a sudden stress into the simulation. In memoryless models (which depend solely on the most immediate past state), the response to this perturbation is abrupt; the transformation of the stability landscape is sudden when the conditions before and after the perturbation are assumed to be invariant. However, in systems with memory, the situation is different. Following a pertur- bation, the stability landscape of a system with memory does not shift instantaneously. Instead, it gradually reverts to a state resembling its original confguration. Sys- 57 Moein Khalighi tems with memory sustain the lingering effects of such disturbances, leading to slow recovery or irreversible alterations in the stability landscape. As a result, the post- perturbation landscape exists in a state that is qualitatively distinct from its original confguration, as Figure 17 conceptually shows. To elaborate, the landscape’s valleys and peaks transform after perturbations, necessitating time for the system to regain its original state. Additionally, the landscape’s potential depth and curvature assume new shapes. Figure 17. A conceptual comparison of stability landscapes in memory-integrated and memoryless systems, on univariate bistable systems. This illustrates the system’s behavior across three phases: pre-, during, and post-perturbation. In memoryless systems, the landscapes before and after perturbations are the same. Conversely, in systems with memory, the landscape evolves both during and after the perturbation, affecting the system’s perturbation response. ��1 and ��2 denote stable states, and �� indicates an unstable state. Here, �0 is the initial time, �� marks the onset of the perturbation, and �� is the fnal time, when the system either returns to the original stable state or transitions to an alternative one. To illustrate these ideas, consider a system initially at the bottom of the right valley in Figure 17, denoted as �0 = ��2 . Suppose a scenario where a system can recover itself after a disturbance. During this perturbation, the system temporarily loses stability at ��2 due to a perturbed function ��, causing a transition towards �� ′ , near the unstable state �� . The potential energy under this perturbation is given by: ∫ ��2 ˙�� = − ���, (29) �� ′ such that for each �� ∈ [�� ′ , ��2 ]: ∫ �� �� = �0 + �(�� − �)��(�, �(� ))��. (30) ��2 58 Impact of memory on ecological systems Here, ��2 represents the initial time when the system is in state ��2 , and �� represents the current time during the transition from the higher state ��2 to the lower state �� ′ . After perturbation, assume that the system will settle in state ��2 ′ , close to the original stable state ��2 . The potential energy in this new state can be determined by: ∫ �� ′ ˙� = − 2 ���, (31) �� ′ ′where each �� ∈ [�� ′ , ��2 ] achieved by: ∫ �� �� = �0 + �(�� − �)�(�, �(� ))��, (32) ��2 ′where �� is the current time as the system transitions from �� to ��2 ′ , and � is a piecewise function of � and ��: { ��(�, �(�)) if � ∈ [��2 , �� ′ ],�(�, �(�)) = (33) � (�, �(�)) if � ∈ [�� ′ , ��2 ],′ ′′ ′where �� and ��2 denote the times at which the system reaches the states �� and ′ , respectively. ��2 Equation (31) reveals that post-perturbation, the system’s potential valley shape cannot fully revert to its original form due to the kernel � and the piecewise function � in the integral equation (32). Notably, when there is no memory in the system, the potential energy (31) after the perturbation depends solely on the function � (�, �(�)): ∫ �� ′ ˙� = − 2 ���, (34) �� ′ ′where each �� ∈ [�� ′ , ��2 ] is given by: ∫ �� �� = �0 + � (�, �(�))��, (35) �� ′ In this case, the system’s states are independent of past states, and the integral starts at �� ′ , the time when the system reaches the state �� ′ , rather than the initial time ��2 . Consequently, only the immediate subsequent state depends on �� ′ , which is determined from the perturbed dynamics, �� = ��, while the following states are�� governed solely by the original function. Since the potential values derived from Equation (34) are identical to those before the perturbation (as they follow the same functional form), the potential landscape remains unchanged and is unaffected by the perturbation. These comparisons demonstrate that memory not only alters the transient dy- namics but can also lead to qualitatively different stability landscapes following perturbations. 59 Moein Khalighi 4.1.3 Impact of memory on resilience and resistance I now examine how memory infuences resilience and resistance using the Herbivory model (3). Focusing on a system residing in a single stable state, I analyze its response to pulse perturbations. The fndings indicate that memory decreases resilience to perturbations, as demonstrated in Figure 18a. Initially, memory may accelerate the recovery process; however, over time, it slows the return to the pre-perturbation state, thereby affecting resilience as defned in Section 2.2.1. Due to its prolonged infuence, long-term memory can extend recovery times—for example, at a 10−3 threshold, memory is associated with an 87% reduction in the recovery rate. Weak perturbation Strong perturbation Figure 18. A comparative view of a bistable herbivory model reaction to pulse perturbation. Panel (a) illustrates that both models with and without memory recover after the perturbation ceases, albeit at differing recovery rates. Specifcally, considering a return to the original state within a tolerance of 10−3 , the incorporation of memory leads to an 87% reduction in the recovery rate. Panel (b) demonstrates that the memoryless model fails to withstand a more intense perturbation and transitions to an alternative stable state. In contrast, the memory helps the system return to its original stable state. In contrast, memory can enhance resistance when the system is near the bottom of its basin of attraction. In such cases, memory acts as a stabilizing force by preventing state transitions following a disturbance. This effect is evident in Figure 18b, where the memory-integrated system withstands a stronger perturbation and returns to its original state, while the memoryless model transitions to an alternative stable state. This stabilization is particularly signifcant in multistable systems, where transient pulse perturbations may mimic more enduring press perturbations, potentially triggering regime shifts. However, this behavior is observed when the system is not near the unstable point. Figure 19 illustrates a scenario in which memory induces a system transition to an alternative stable state. In this case, the system follows a gradual trajectory from the unstable point to the high-population stable state, but the perturbation redirects the trajectory toward the alternative stable state. Because the system with memory recovers more slowly and remains near the unstable point for longer, even a slight 60 Impact of memory on ecological systems perturbation can shift it to an alternative stable state, thereby reducing resistance. Figure 19. Memory-induced transition to an alternative stable state. The slower dynamics in the memory-integrated system keep it near the unstable point, allowing a slight perturbation to shift it to a different stable state. Additionally, the deformation of the stability landscape driven by memory tends to fatten the peaks near unstable points, as indicated by the reduced slope of the potential energy in Figure 16b. This fattening results in slower dynamics and prolonged instability, as illustrated in Figure 20. In systems with memory, unstable points and the boundaries of attraction basins become less sharp, causing the system to linger near these points for an extended period. Consequently, memory can also lead to prolonged instability. Figure 20. Memory-induced prolonged instability. The fattened potential near unstable points slows the dynamics, causing the system to remain near unstable equilibria for longer periods. These contrasting effects raise important questions about how memory modulates system behavior under more complex types of perturbations, such as periodic or 61 Moein Khalighi stochastic disturbances. While memory generally slows system dynamics, it does not always buffer against perturbations; in some cases, it can make the system more responsive and lead to more pronounced shifts. As discussed in Section 4.1.2, the stability landscape can be dynamic when facing multiple perturbations. For instance, Figure 21 shows two opposing scenarios under stochastic perturbations. In one case, memory acts as a buffer, enhancing both resilience and resistance by ensuring a gradual, controlled response that prevents sudden shifts. In the other case, memory leads to more pronounced fuctuations and increases the propensity for the system to shift to an alternative stable state. (a) (b) Figure 21. Contrasting effects of memory under stochastic perturbations. (a) Memory buffers fuctuations and prevents transitions to alternative stable states. (b) Memory amplifes fuctuations, making the system more prone to shifts. 4.2 Multivariate systems This section explores the infuence of memory on multi-dimensional systems, with a specifc emphasis on ecological models as outlined in Publications II. In ecological systems, memory effects—commonly referred to as ecological memory [6]—play a crucial role. Despite the advancements in mathematical modeling presented in Section 2.2.2, particularly in community dynamics modeling, a comprehensive under- standing of ecological memory’s role in microbial ecology remains elusive. Microbial communities are a key focus because memory effects can be inte- grated into multi-species models, signifcantly impacting community composition, stability, and responses to environmental changes. In these systems, ecological mem- ory—which may exhibit power-law behaviors due to factors such as intracellular processes and phenotypic heterogeneity [144]—critically shapes how communities respond to perturbations and how alternative community states emerge. Exploring memory effects in high-dimensional microbial systems is challenging. For example, even in a simple two-species community model, the stability landscape becomes three-dimensional, necessitating increased computational efforts or even 62 Impact of memory on ecological systems rendering the landscape undefned [37]. Consequently, methods used in univariate model analyses are often impractical here. Instead, one must analyze system dynamics under various conditions and memory levels and then contrast the outcomes. Therefore, I focus on the specifc model (7), which can simulate both real and synthetic data, such as those used in studies of human gut microbiota [51]. In this section, I examine a bistable two-variable of this LV competition model (7), where the dominance of one species over the other depends on initial conditions. For simplicity, I denote the two species as � and � , abstracting from specifc species identities while retaining the core dynamics of gut microbiota community behavior. The next section will explore how memory affects the dynamics of this system. 4.2.1 Impact of memory on a model in the relaxed phase In this subsection, no perturbations are applied, allowing us to isolate the impact of memory on dynamics. As illustrated in Figure 22a, the model’s state space reveals a compelling structure under given conditions. The system has three equilibrium points apart from the stable state at the origin; one of these is unstable, while the other two are stable, indicating bistability that ultimately leads to the extinction of one species. Figure 22b demonstrates the impact of incorporating coherent memory on com- munity trajectories. Although memory does not alter the equilibrium points, it slows down system dynamics (Figure 23a and 23d). Consequently, with identical initial species abundances, the community reaches the same extinction outcome, albeit at different times. This behavior is observed under consistent, unperturbed conditions. The characteristics of a univariate system with memory (discussed in Section 4.1) can be generalized into how coherent memory infuences a bistable multivariate model in both the relaxed and perturbed phases. However, the interplay between varied memory and perturbations is complex, with outcomes that are unpredictable and context-dependent. Hence, I only focus on examining some behaviors of systems with varied memory. Impact of memory types on systems’ behavior Coherent memory behaves in a similar manner to that in univariate systems (see Section 4.1). Insights from Publication II indicate that coherent memory acts as a stabilizing force—it can buffer against fuctuations by ensuring a gradual, controlled response to changes, thereby promoting coexistence or prolonging the dominance of a species, depending on past states. Conversely, when acting as an amplifer, memory may drive a sudden shift toward an alternative stable community composition, possibly triggering species extinction or dominance shifts. In order to avoid repeating previous results, I focus here on the effects of varied memory within the context of ecological communities. Introducing varied memory into the system unveils complex scenarios that trans- 63 Moein Khalighi Figure 22. Illustration of the state space of a two-species competition LV model (7). This model features two stable states (indicated by flled circles) and one unstable state (indicated by an unflled circle). Dark and light gray curves show distinct trajectories converging from various initial species abundances to the stable states, while the blue and red curves specifcally show how adding memory to the species’ behavior changes their path, starting from the same starred points. Panel (a) shows the state space without memory, panel (b) with coherent memory, and panels (c) and (d) with varied memory: in (c) only species � has memory, and in (d) only species � has memory, illustrating how varied memory can alter trajectory convergence. form community behavior, as visualized in Figures 22c-d. Notably, even with identical initial abundances under unperturbed conditions, the community may reach different extinction outcomes. This highlights the importance of historical infuence on system behavior. For example, Figure 23 shows that memory type can change the speed of dynamics. While a multivariate system with coherent memory typically slows dynamics in the absence of perturbations (Figures 23a and 23d), varied memory can either decelerate (Figures 23c and 23e) or accelerate (Figures 23b and 23f) dynamics, sometimes even driving the community toward an alternative stable state. Similarly, varied memory can lead to prolonged transient states where the com- 64 Impact of memory on ecological systems Population without memory Population with memory (a) (b) (c) (d) (e) (f) Figure 23. Community dynamics under varied memory strength levels. The top row shows trajectories from initial condition (0.5, 0.6), and the bottom row corresponds to (0.7, 0.3). Each column represents a different memory setting, illustrating that varied memory can either slow down or speed up dynamics and, in some cases, lead to alternative stable states. munity remains unstable before eventually shifting abruptly to stability (Figure 24). Using the ball on a hill and valleys metaphor, memory fattens the hill near an un- stable point, allowing the system to linger in a transient state—a buffer period that might be crucial for adaptation in changing environments. Note that, as discussed in univariate models, coherent memory can also function as an amplifer. Given the greater fexibility of varied memory compared to coherent memory, we will see in the next section that varied memory can act as both a buffer and an amplifer. 4.2.2 Impact of varied memory on a model under perturbation Varied memory can disrupt predictable patterns and steer the system toward unex- pected outcomes. Its presence indicates a complex interplay between stability and change, highlighting the intricate nature of ecological dynamics. In this context, memory—in its varied forms—plays a signifcant role in determining how ecological communities respond to both internal and external pressures, thereby infuencing their evolution and adaptation over time. 65 Moein Khalighi Population without memory Population with memory Figure 24. Varied memory prolongs transient instability. Memory fattens the stability landscape near unstable points, resulting in extended periods of instability before the system eventually converges to a stable state. Varied memory can either enhance or reduce resilience and resistance. Figure 25a shows a case where memory increases resistance by preventing a shift to an alternative stable state following a pulse perturbation, whereas the memoryless system transitions to a new state. In contrast, Figure 25b illustrates a scenario where memory reduces resistance by causing the system to shift to a new basin of attraction while the memoryless system remains in its original basin of attraction. Thus, in response to a pulse perturbation, varied memory may either cause a shift in community composition or help resist such a shift. Moreover, varied memory can yield contradictory responses to stochastic noise. In some cases, it acts as a buffer (Figure 25c), leading to gradual responses and enhanced stability; in others, it amplifes fuctuations (Figure 25d), making the system more prone to abrupt shifts than in the absence of memory. 4.2.3 Exploring the frontiers of ecological memory This chapter partially explored the third and fourth research questions by examining how memory enhances our understanding of system dynamics in ecology. For exam- ple, memory-induced fattening of stability landscapes explains prolonged instability near tipping points, and the concept of a time-dependent stability landscape supports dynamic stability, aiding in gradual ecosystem recovery after disturbances. While my approach offers novel insights, many fundamental questions remain before the potential of ecological memory in predictive ecological modeling can be fully harnessed. Future research should prioritize determining the precise conditions under which memory shapes system behavior, as well as exploring its role in specifc ecological contexts—such as gut microbiota dynamics—which could substantially enhance model accuracy. 66 Impact of memory on ecological systems Population without memory Population with memory (a) (b) (d)(c) Figure 25. Impact of varied memory on community responses to perturbations. (a) Under a pulse perturbation, memory prevents a transition to an alternative stable state, thereby enhancing resistance. (b) Alternatively, memory can induce a shift to a new basin of attraction, reducing resistance. (c) In response to stochastic noise, memory may buffer fuctuations, and (d) in another scenario, memory can amplify fuctuations, resulting in exacerbated fuctuations. A central issue is how best to represent ecological memory. It is still unclear whether memory should be modeled as a static, long-term property or whether a time- varying formulation would be more appropriate. Moreover, identifying the optimal mathematical form for the memory kernel is critical; for instance, one must decide whether to favor singular power-law kernels or to explore alternative non-singular functional forms. Another open question concerns the content of ecological memory. Should models consider only past system states, or is it necessary to also incorporate the history of external perturbations and interspecies interactions? Investigating the relative importance of these memory components—including lag effects, where certain past events exert a disproportionately strong infuence—could yield refned insights into community dynamics. For example, while recent events may generally carry greater weight, some historical events might trigger amplifed or delayed responses. A deeper understanding of the underlying biological mechanisms is also essential. It remains to be determined whether ecological memory primarily refects reliance on historical events or whether it is driven by intracellular inertia, such as variability in species’ lag phases, which describe periods of slow or negligible growth before 67 Moein Khalighi exponential expansion due to physiological adaptations. Clarifying these mechanistic details will guide the development of more biologically realistic models. Practically, future studies should evaluate model accuracy across diverse ecological contexts to determine the most effective and parsimonious ways to incorporate memory, poten- tially by comparing alternative formulations that demonstrably improve predictive power. Finally, empirical investigations are crucial to characterizing ecological memory in real-world communities. Such studies will not only clarify memory dynamics but also pave the way for developing techniques to manipulate and manage ecosystems. By integrating theoretical analysis, simulations, and empirical validation, we can fully unlock the potential of ecological memory to deepen our understanding of community dynamics and enhance our ability to predict and manage ecological systems. Furthermore, while the initial focus has been on community dynamics, the methodologies developed here have broad applicability across diverse complex systems, as will be demonstrated in the following chapter on epidemiological models. 68 5 Impact of memory on epidemiological systems In this chapter, I examine how incorporating memory into epidemiological models can alter disease transmission dynamics, directly addressing the third research question. Traditional epidemic models assume a memoryless process, where the future state depends solely on the present; however, real-world epidemics are often infuenced by past events over extended periods [145]. In Publications III and IV, I extended ad- vanced compartmental models (see Section 2.3.2) by introducing memory kernels into differential equations. These kernels decay over time with long-range dependencies, distinct from the power-law event distributions seen in epidemics [146; 147]. Incorporating such memory kernels is a reasonable approach because epidemics are affected by factors like variable incubation periods, delayed effects of interventions, and heterogeneous contact patterns [145; 22]. Moreover, research has demonstrated that the effciency of control measures—such as soft quarantine strategies—is closely tied to these long-range temporal dependencies [148; 124]. By integrating memory kernels into epidemiological models, I aim to account for the time-dependent impacts of interventions and the multi-scale nature of epidemic spread. In the following sections, I illustrate how this modeling advancement alters key properties of epidemic dynamics, thus improving our predictive capabilities of disease transmission and addressing the fourth research question. 5.1 Impact of memory on a simple model I begin by discussing the impact of memory on the basic SIR model (introduced in Section 2.3) before considering more advanced cases. 5.1.1 Memory and fnal epidemic size Memory effects can infuence the dynamics and fnal epidemic (or outbreak) size in the SIR model (see Section 2.3.1) [22]. Incorporating coherent memory through 69 Moein Khalighi fractional calculus, the model is expressed as �� � = −���, �0 ��� 0 � = ��� − ��, (36) �� � = ��, �0 where lower values of � (i.e., longer memory lengths) increase the epidemic threshold and reduce the outbreak size, thereby making the system more robust against the spread of infections [22]. Over time, the infuence of memory decays, and the model gradually reverts toward the characteristics of a memoryless system. These results suggest that coherent memory can help control and mitigate disease spread by elevating resistance thresholds and reducing epidemic sizes. However, these conclusions primarily apply to scenarios with coherent memory in relatively simple compartment models. In models with varied memory or more advanced structures (see Section 5.2), such broad conclusions may not hold true. 5.1.2 Memory and basic reproduction number In epidemic modeling, the basic reproduction number, ℛ0, is a central concept (see Section 2.3.1). Although memory can impact epidemic dynamics in the SIR model, it may not necessarily alter ℛ0. In the classical case (12) without memory effects (see Section 2.3.1), ℛ0 = �/�. When fractional derivatives of the same order are employed—as in (36) or its incommensurate version—ℛ0 remains �/�. Thus, adding memory in this manner does not change the basic reproduction number in the SIR model or even in more advanced models (as in Publication III, see Section 5.2.1). Nonetheless, some modifcations introduced by fractional derivatives can affect ℛ0, as discussed further. Adapting parameter scaling When incorporating fractional-order derivatives into classical models, one may con- sider reevaluating parameter scaling as stressed in [149]. Aligning parameters with the fractional derivative’s order maintains consistency and ensures that the units remain meaningful. For example, in a modifed SIR model that incorporates varied memory, the equations become: ��� �0 � = −��� ��, ��� � = ��� �� − ��� �, (37)�0 ��� �0 � = ��� �. In this framework, the basic reproduction number becomes dependent on the derivative orders, highlighting its memory-dependent nature. To determine ℛ0, we can linearize 70 Impact of memory on epidemiological systems the system around the disease-free equilibrium (typically � ≈ 1) and focus on the � equation. In the early stages, the � equation reduces to ( ) ��� ��� − ��� � ≈ �. �0 By analogy with the classical case, where ℛ0 = �/�, it is natural to defne the basic reproduction number for the fractional model as ( )�� ��� � ℛ0 = = . ��� � This effective reproduction number refects the impact of memory in the infected compartment. Note that if no parameter rescaling is applied (i.e., �� = �� = �� = 1), the classical result ℛ0 = �/� is recovered. For an epidemic to fourish, the effective growth rate must be positive: ��� > ��� , or equivalently, � �� > 1. Thus, this fractional formulation of ℛ0 highlights how ��� memory effects can alter transmission dynamics. This approach not only preserves mathematical rigor and enhances the model’s ability to capture the complex dynamics of disease transmission and progression, but it also highlights that while adjusting parameter scaling can provide a more robust quantitative framework, many choose not to modify it. In numerous cases, the omission of this adjustment does not signifcantly compromise the model’s ef- fectiveness, especially when the focus is on qualitative dynamics rather than strict quantitative accuracy, though it remains an important consideration for achieving optimal precision. More advanced compartmental systems introduce further challenges [150; 106], such as evaluating a memory-dependent ℛ0 and analyzing stability, and ensuring the existence and uniqueness of solutions. These issues require specialized techniques and tools, especially because a memory-dependent ℛ0 affects both local and global stability analyses at the disease-free equilibrium, which I applied in Publication IV. 5.1.3 Challenges in using varied memory Using fractional order derivatives introduces specifc challenges, particularly when accounting for varied memory, or incommensurate derivative orders [106]. These challenges emerge when replicating real-world data, such as modeling epidemic diseases [124; 150]. In traditional epidemic models, like the SIR model, the total population is typically assumed to be constant, where �+� +� = � . This assumption simplifes the analysis and supports many studies. However, incorporating fractional-order derivatives introduces variable rates of change for each compartment. This change means that the sum of the compartmental populations may no longer equal the total population, which complicates both the interpretation and stability analysis of the model [151]. 71 Moein Khalighi Consequently, traditional methods for determining boundedness and stability might fail, prompting the development of novel techniques. Thus, recent advance- ments in the feld have illuminated ways to manage such complexities [152]. Under specifc conditions, the exponential boundedness of solutions for systems with in- commensurate fractional derivatives has been demonstrated [153]. In this context, exponential boundedness means that there exist constants � > 0 and � such that |�(�)| ≤ ���� for all �. This condition guarantees that even if the solution grows over time, it does so at a rate no faster than ����, ensuring the system remains controlled. Moreover, fractional-order derivatives add complexity by making the system’s dynamics harder to interpret than those of traditional integer-order models [154]. Fully understanding the effects of varied memory and the interactions between compartments requires a deep grasp of fractional calculus and its applications in epidemiological modeling. Despite these challenges, fractional-order modeling offers insights into the impact of memory on epidemic dynamics (see Section 5.2.1 and 5.2.2). My works in Publications III and IV address these issues through theoretical analysis, establishing a robust foundation regarding the existence, uniqueness, boundedness, and stability of advanced epidemic models with varied memory. 5.2 Assessing memory in extended models This section addresses the fourth research question for specifc cases in epidemi- ology to illustrate how integrating memory into models enhances our theoretical understanding of complex systems. I have investigated how varying memory levels impact the system’s dynamics while keeping model parameters constant. The next step involves fne-tuning both memory and model parameters simultaneously to better align with observed data. Additionally, I explored how memory in each compartment affects model parameters, system dynamics, the ftting process, and the basic reproduction number–all of which are achievable through my solver package (Section 3.3.4). Utilizing my solver package, which integrates seamlessly with various optimiza- tion tools, I adjusted both the fractional derivative orders and the model parameters concurrently. Numerical computations were conducted using FdeSolver.jl (outlined in Publication I) to numerically solve the fractional differential equations. Parameter estimation and Bayesian inference were carried out using the Turing.jl package [125], specifcally through Markov Chain Monte Carlo (MCMC) sampling. Hamiltonian Monte Carlo (HMC) sampling was employed, supported by informed prior distribu- tions. Variances, for instance, were assigned inverse gamma priors to ensure positivity and proper regularization. Other model parameters were represented using truncated normal distributions to impose biologically or physically meaningful constraints. 72 Impact of memory on epidemiological systems For optimizing fractional model parameters, I utilized the Optim.jl package [155], constraining fractional derivative orders initially within biologically plausible ranges to identify feasible parameter spaces. I chose the Limited-memory BFGS algo- rithm [155] for its effciency and robustness in managing constrained, nonlinear optimization problems. After exploratory analyses with Turing.jl to determine suit- able priors, these priors guided the parameter ranges used in Optim.jl. This enabled the fne-tuning of fractional derivative orders, leading to optimal model fts to the ob- served data. Model outcomes were evaluated based on the root mean square deviation (RMSD) from actual data, defned as: ⎯⎸ � ⎸ ∑ ⎷ 1RMSD(�, �̂) = (�� − �̂�)2 , (38) � �=1 where, � denotes the number of data points, � denotes the estimated values, and �̂ denotes the observed values. In Publication I, I presented the detailed methodology for ftting fractional-order models to real-world data, accommodating both commensurate and incommensurate orders. I demonstrated the effectiveness of this approach by applying it to real datasets. Specifcally, for COVID-19 data in Publications I and IV, I adjusted the order of derivatives as a measure of memory in the system. Similarly, in Publication III, I determined the derivative orders when modeling the spread of Ebola. The results highlight how incorporating varying levels of memory enhances model fexibility and accuracy. However, introducing additional fexibility can also increase the risk of overftting. Therefore, future research could explore techniques, such as cross- validation or other appropriate criteria , to further mitigate overftting and optimize the balance between memory capacity and model parameters. Such studies will enhance the reliability and interpretability of fractional order models in practical applications. 5.2.1 Ebola transmission dynamics In Publication III, I developed a fractional differential model, based on Section 2.3.2, to ft long-term Ebola epidemic data of cumulative confrmed cases from Guinea, Liberia, and Sierra Leone (World Health Organization, 438 days). This model incor- porates dynamic population adjustments and memory effects to simulate the disease’s complex transmission patterns, including a notable data shift between days 212–216. I evaluated four models of increasing complexity: Model 1: A simplifed version of Equation (13) with a constant population � , adjusted birth rate �1 − �2� , and death rate �1 + �2� with a single constant � [75], reducing model fexibility. Model 2: An ODE system with integer orders (13), providing dynamic adjustments but limited by fxed derivative orders. 73 Moein Khalighi Model 3: An incommensurate fractional differential equation system: �� �ℎ �� �� ��� �(�) = (�1 − �2� )� − �� − �� − �� − �� − (�1 + �2� )�, � � � � �� �ℎ �� �� ��� �(�) = �� + �� + �� + �� − �� − (�1 + �2�)�, � � � � ��� �(�) = �� − (�1 + � + � + �1 + �2� )�, ��� �(�) = �1� + �2� − (�3 + �1 + �2� )�, (39) ��� �(�) = �� − (�1 + �)�, ��� �(�) = � � − (�2 + �2 + �1 + �2� )�, ��� �(�) = �1� + �2� − ��, ��� �(�) = �3� − (�1 + �2� )�, allowing variable rates of change for improved data ftting. Model 4: A hybrid approach combining Models 2 and 3, transitioning from ODE to FDE at day 216 to incorporate memory effects post-data shift. Prior to that point, no memory effects are assumed, refecting a period with consistent social behavior or policies. After day 216, memory effects are introduced to deliver the impact of newly emerging or altered social behaviors. Unlike the more traditional Model 1, which assumes a constant population, Mod- els 2–4 dynamically account for demographic changes (e.g., births, deaths, and migration), thereby altering the number of susceptible individuals and infuencing disease transmission. The basic reproduction number ℛ0 for Models 2, 3, and 4 aligns with Section 2.3.2, while Model 1’s ℛ0 differs due to its simplifed parameters [75]. Fitting real-world data illustrates the utility of my solver package in epidemiological modeling. Model performance was assessed using RMSD 38, with Model 3 outperforming Model 2. This indicates that incommensurate fractional derivatives, that is, varied memory effects, yield a more precise representation of the Ebola transmission dynam- ics. Sensitivity analysis revealed that the dynamics of Model 3 were more sensitive to derivative order changes in compartments � (infected individuals), � (buried in- dividuals), and � (dead but not buried individuals) compartments, whereas Model 4 exhibits high sensitivity in � and � (exposed individuals) compartments. Of all the models, the hybrid Model 4, incorporating optimized fractional orders, achieved the highest accuracy. 5.2.2 COVID-19 transmission dynamics In Publication IV, I developed a fractional-order compartmental COVID-19 model, explained in Section 2.3.2, to address limitations inherent in traditional integer-order 74 Impact of memory on epidemiological systems models, which often overlook environmental transmission and memory effects [76; 124]. The proposed fractional model, expressed by equations (40), enhances predictive accuracy and offers deeper insights into COVID-19 transmission dynamics. As discussed in Section 5.1.3, all model parameters—except for �1, �2, and �, which are dimensionless—have dimensions of 1/��� , where � corresponds to the compartments (�, �, ��, �� , �, �, � ). Thus, the model is: ��� �(�) = Λ�� − ��� Φ1 − ��� Φ2 + ��� �(�) − ��� �(�),1 2 ��� �(�) = ��� Φ1 + ��� Φ2 − ��� �(�) − ��� �(�) − ��� �(�),1 2 ��� ���� ��(�) = (1 − �)���� �(�) − (���� + ���� )��(�) − � ��(�),� ��� ���� �� (�) = ����� �(�) − (���� + ���� )�� (�) − � (40)� �� (�), ��� �(�) = ��� �� (�) + ��� ��(�) − ��� �(�),� � ��� �(�) = ��� (�� (�) + ��(�)) − ��� �(�), ��� � (�) = ���� ��(�) + ���� �� (�) − ��� � (�).� I calibrated this model using COVID-19 data from South Africa, specifcally daily new confrmed cases � , recoveries �, and deaths �, sourced from the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University [156]. The fractional-order model demonstrated a superior ft to real data compared to its integer-order counterpart, evidenced by lower RMSD values (38). This improvement demonstrates the model’s fexibility and capacity to capture the complexities of COVID-19 transmission, especially those related to environmental factors. The mathematical analysis presented in Publication IV employs key aspects of fractional calculus, including the Caputo fractional derivative and the Mittag-Leffer function [157; 17], which are integral to understanding fractional-order dynamics. The existence and uniqueness of solutions are established via fxed-point theory and the Picard–Lindelöf theorem [158], reformulating the system of equations to prove these properties. Furthermore, it confrms the boundedness and positive invariance of solutions within a biologically meaningful domain, ensuring stability and relia- bility. The basic reproduction number ℛ0 is computed using the next generation matrix method, correcting oversights in previous modeling approaches. The resulting expression, [ ( ) ����� (1 − �)���� ℛ0 = (Λ/�)�� ��� +2 ����� ����� ( )] �� �� ����� ��� (1 − �)���� +��� � + , (41)1 ��� ��� � ����� ����� � 75 Moein Khalighi where ��� ��� = ��� +��� +��� = ���� +���� +� = ���� +���� +� �� , ��� � , ��� � , is memory-dependent and differs from the reproduction number in the memoryless model (16). To quantify the effect of small changes in both a parameter and its associated derivative order � , the forward normalized sensitivity index is introduced: �ℛ0 � �ℛ0 = .� �� ℛ0 A positive index signifes that increasing � leads to a higher ℛ0, while a negative index indicates the opposite effect. Sensitivity analysis of the basic reproduction number ℛ0 revealed that several parameters, including those related to environmental contamination (e.g., �1, �� , and ��), signifcantly impact ℛ0. This analysis demonstrated that reducing environmental transmission could substantially lower ℛ0, suggesting effective intervention strategies such as enhanced sanitation and hygiene practices. Additionally, the sensitivity indices indicated that while most parameters had a positive infuence on ℛ0, some, like the order derivatives, showed varied impacts, providing deeper insights into how different factors contribute to the disease’s spread. The fndings highlight the model’s utility in guiding public health policies to control the outbreak. Of particular note is how fractional orders affect ℛ0. While most parameters in the fractional model exhibit similar behavior to those in the integer-order model, certain parameters, such as � (the progression rate from exposed to infected), can negatively infuence ℛ0. Moreover, the sensitivity index of �� (the derivative order of the exposed compartment) is negative, whereas those of other orders are predominantly positive. These fndings suggest that fractional derivative orders infuence ℛ0 in nuanced ways. Overall, the basic reproduction number for fractional-order models is slightly higher than in integer-order models, and environmental transmission reduction remains a robust control strategy. Finally, local and global stability analyses verify that the disease-free equilibrium (see Section 2.3.1) is locally asymptotically stable under specifc conditions and glob- ally stable when ℛ0 ≤ 1. These proofs leverage Lyapunov functions and Ulam–Hyers stability principles [159]. In summary, these fndings can advance fractional epidemiological modeling by introducing theoretical tools to analyze solution existence, uniqueness, boundedness, and stability. The fractional-order compartmental model developed here provides a framework for examining the transmission dynamics of COVID-19, with particular emphasis on the role of environmental pathogens. Model robustness is confrmed through theoretical analysis and validation against real-world data, stressing the critical importance of environmental sanitation and targeted public health interventions for effective pandemic control. 76 6 Conclusion The primary aim of this thesis was to examine the role of memory effects in dynamical systems, addressing a gap in traditional modeling that often relies on memoryless approaches. The work investigated how memory infuences complex dynamics by tackling key questions: how to systematically incorporate memory into models, how memory alters ecological properties and epidemiological quantities, and how memory-based approaches can enhance our theoretical understanding of complex systems. This thesis conceptualized and systematically incorporated memory into both ecological and epidemiological models. By utilizing fractional derivatives, I modeled long-term memory effects that enhanced system calibration. This approach provided a robust foundation for integrating memory into both simple and complex scenarios, ensuring a balance between computational tractability and model comprehensiveness. Despite the inherent challenges, advanced numerical techniques and modern open- source computational tools were employed to bridge the methodological gap left by conventional approaches. My exploration revealed the impacts of memory on system behavior across various contexts. In ecological models, memory was shown to infuence stability, resilience, and resistance, providing a quantitative and theoretical depiction of how ecosystems may respond to disturbances and affect overall stability. Similarly, epidemiological models are improved by incorporating memory effects, as demonstrated in studies of Ebola and COVID-19 transmission. This approach explicitly quantifes the impact of memory on key epidemiological metrics, thereby representing the complexities of real-world transmission dynamics and informing the development of more effective outbreak control measures. While my approach in this thesis is broadly applicable, future research should focus on identifying memory mechanisms in specifc contexts, such as microbial communities, and on exploring various memory types. Although incorporating different memory types is mathematically and phenomenologically complex, doing so could further improve prediction accuracy. In summary, this dissertation advances the modeling of dynamical systems by incorporating memory effects and addressing specifc research questions of applied relevance. This work expands our understanding of system dynamics and prepares the way for new research questions. 77 Declaration on the use of AI I used ChatGPT and Grammarly to enhance the writing style, grammar, and spelling. I carefully reviewed and revised the text as needed, taking full responsibility for the fnal content. 78 List of references [1] Daniel A. Power, Richard A. Watson, Eörs Szathmáry, Rob Mills, Simon T. Powers, C. Patrick Doncaster, and Bła ˙Zej Czapp. What can ecosystems learn? Expanding evolutionary ecology with learning theory. Biology Direct, 10(1):69, 2015. [2] Chaitanya S. Gokhale, Stefano Giaimo, and Philippe Remigi. Memory shapes microbial popula- tions. PLOS Computational Biology, 17(10):1–16, 2021. [3] Richard A. Watson and Eörs Szathmáry. How can evolution learn? Trends in Ecology & Evolution, 31(2):147–157, 2016. [4] David G. Angeler, Hannah B. Fried-Petersen, Craig R. Allen, Ahjond Garmestani, Dirac Twidwell, Wen-Ching Chuang, Victoria M. Donovan, Tarsha Eason, Caleb P. Roberts, Shana M. Sundstrom, and Carissa L. Wonkka. Chapter one - adaptive capacity in ecosystems. In David A. Bohan and Alex J. Dumbrell, editors, Resilience in complex socio-ecological systems, volume 60 of Advances in Ecological Research, pages 1–24. Academic Press, 2019. [5] Janne Bengtsson, Per Angelstam, Thomas Elmqvist, Urban Emanuelsson, Carl Folke, Margareta Ihse, Fredrik Moberg, and Magnus Nyström. Reserves, resilience and dynamic landscapes. AMBIO: a Journal of the Human Environment, 32(6):389–396, 2003. [6] Garry D. Peterson. Contagious disturbance, ecological memory, and the emergence of landscape pattern. Ecosystems, 5:329–338, 2002. [7] Jill F. Johnstone, Craig D. Allen, Jerry F. Franklin, Lee E. Frelich, Brian J. Harvey, Philip E. Higuera, Michelle C. Mack, Ross K. Meentemeyer, Margaret R. Metz, George L. W. Perry, et al. Changing disturbance regimes, ecological memory, and forest resilience. Frontiers in Ecology and the Environment, 14(7):369–378, 2016. [8] Reena Debray, Robin A. Herbert, Alexander L. Jaffe, Alexander Crits-Christoph, Mary E. Power, and Britt Koskella. Priority effects in microbiome assembly. Nature Reviews Microbiology, 20(2): 109–121, 2022. [9] Irine Ronin, Naama Katsowich, Ilan Rosenshine, and Nathalie Q Balaban. A long-term epigenetic memory switch controls bacterial virulence bimodality. eLife, 6:e19599, 2017. [10] Pamela Lyon. The cognitive cell: bacterial behavior reconsidered. Frontiers in Microbiology, 6, 2015. [11] Andreas H. Schweiger, Isabelle Boulangeat, Timo Conradi, Matt Davis, and Jens-Christian Svenning. The importance of ecological memory for trophic rewilding as an ecosystem restoration approach. Biological Reviews, 94(1):1–15, 2019. [12] Kiona Ogle, Jarrett J. Barber, Greg A. Barron-Gafford, Lisa Patrick Bentley, Jessica M. Young, Travis E. Huxman, Michael E. Loik, and David T. Tissue. Quantifying ecological memory in plant and ecosystem processes. Ecology Letters, 18(3):221–235, 2015. [13] Malcolm S. Itter, Jarno Vanhatalo, and Andrew O. Finley. Ecomem: an R package for quantifying ecological memory. Environmental Modelling & Software, 119:305–308, 2019. [14] Blas M. Benito, Graciela Gil-Romera, and H. John B. Birks. Ecological memory at millennial time-scales: the importance of data constraints, species longevity and niche features. Ecography, 43(1):1–10, 2020. [15] Valentin Schaefer. Alien invasions, ecological restoration in cities and the loss of ecological memory. Restoration Ecology, 17(2):171–176, 2009. 79 Moein Khalighi [16] Igor Podlubny. Fractional differential equations: An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. Elsevier, 1998. [17] Anatoly Aleksandrovich Kilbas, Hari M. Srivastava, and Juan J. Trujillo. Theory and applications of fractional differential equations, volume 204. Elsevier, 2006. [18] Maolin Du, Zaihua Wang, and Haiyan Hu. Measuring memory with the order of fractional derivative. Scientifc Reports, 3(1):3431, 2013. [19] HongGuang Sun, Yong Zhang, Dumitru Baleanu, Wen Chen, and YangQuan Chen. A new collec- tion of real world applications of fractional calculus in science and engineering. Communications in Nonlinear Science and Numerical Simulation, 64:213–231, 2018. [20] Kottakkaran Sooppy Nisar, Muhammad Farman, Mahmoud Abdel-Aty, and Jinde Cao. A review on epidemic models in sight of fractional calculus. Alexandria Engineering Journal, 75:81–113, 2023. [21] Yuli Chen, Fawang Liu, Qiang Yu, and Tianzeng Li. Review of fractional epidemic models. Applied Mathematical Modelling, 97:281–307, 2021. [22] Meghdad Saeedian, Moein Khalighi, Nahid Azimi-Tafreshi, Gholamreza Jafari, and Marcel Ausloos. Memory effects on epidemic evolution: the susceptible-infected-recovered epidemic model. Physical Review E, 95:022409, 2017. [23] Romualdo Pastor-Satorras, Claudio Castellano, Piet Van Mieghem, and Alessandro Vespignani. Epidemic processes in complex networks. Reviews of Modern Physics, 87(3):925, 2015. [24] Marten Scheffer. Critical Transitions in Nature and Society. Princeton University Press, 2009. ISBN 9780691122045. [25] V. Grimm and Christian Wissel. Babel, or the ecological stability discussions: an inventory and analysis of terminology and a guide for avoiding confusion. Oecologia, 109(3):323–334, 1997. [26] Didier Gonze, Katharine Z. Coyte, Leo Lahti, and Karoline Faust. Microbial communities as dynamical systems. Current Opinion in Microbiology, 44:41 – 49, 2018. [27] Marten Scheffer and Stephen R. Carpenter. Catastrophic regime shifts in ecosystems: linking theory to observation. Trends in Ecology & Evolution, 18(12):648–656, 2003. [28] Marten Scheffer, Steve Carpenter, Jonathan A. Foley, Carl Folke, and Brian Walker. Catastrophic shifts in ecosystems. Nature, 413(6856):591–596, 2001. [29] Timothy M. Lenton. Early warning of climate tipping points. Nature Climate Change, 1(4): 201–209, 2011. [30] Dave Hodgson, Jenni L. McDonald, and David J. Hosken. What do you mean, ‘resilient’? Trends in Ecology & Evolution, 30(9):503–506, 2015. ISSN 0169-5347. [31] Laurent Philippot, Bryan S. Griffths, and Silke Langenheder. Microbial community resilience across ecosystems and multiple disturbances. Microbiology and Molecular Biology Reviews, 85 (2):10.1128/mmbr.00026–20, 2021. [32] Marten Scheffer, Jordi Bascompte, William A. Brock, Victor Brovkin, Stephen R. Carpenter, Vasilis Dakos, Hermann Held, Egbert H. Van Nes, Max Rietkerk, and George Sugihara. Early- warning signals for critical transitions. Nature, 461(7260):53–59, 2009. [33] Didier Gonze, Leo Lahti, Jeroen Raes, and Karoline Faust. Multi-stability and the origin of microbial community types. The ISME journal, 11(10):2159–2166, 2017. [34] Leah Edelstein-Keshet. Mathematical models in biology. SIAM, 2005. [35] M. Scheffer, S. H. Hosper, M. L. Meijer, B. Moss, and E. Jeppesen. Alternative equilibria in shallow lakes. Trends in Ecology & Evolution, 8(8):275–279, 1993. [36] Beatrix E Beisner, Daniel T Haydon, and Kim Cuddington. Alternative stable states in ecology. Frontiers in Ecology and the Environment, 1(7):376–382, 2003. [37] Pablo Rodríguez-Sánchez, Egbert H. van Nes, and Marten Scheffer. Climbing Escher’s stairs: A way to approximate stability landscapes in multidimensional systems. PLOS Computational Biology, 16(4):1–16, 2020. [38] Vasilis Dakos and Sonia Kéf. Ecological resilience: what to measure and how. Environmental Research Letters, 17(4):043003, 2022. 80 LIST OF REFERENCES [39] Peter Turchin. Complex population dynamics: a theoretical/empirical synthesis (MPB-35). Princeton University Press, 2013. [40] Jordi Bascompte. Extinction thresholds: insights from simple models. In Annales Zoologici Fennici, pages 99–114. JSTOR, 2003. [41] Alex Tsoularis and J. Bruce Wallace. Analysis of logistic growth models. Mathematical Bio- sciences, 179(1):21–55, 2002. [42] Imanuel Noy-Meir. Stability of grazing systems: an application of predator-prey graphs. The Journal of Ecology, pages 459–481, 1975. [43] Robert M. May. Thresholds and breakpoints in ecosystems with a multiplicity of stable states. Nature, 269(5628):471–477, 1977. [44] Ashley Shade, Hannes Peter, Steven Allison, Didier Baho, Mercé Berga, Helmut Buergmann, David Huber, Silke Langenheder, Jay Lennon, Jennifer Martiny, Kristin Matulich, Thomas Schmidt, and Jo Handelsman. Fundamentals of microbial community resistance and resilience. Frontiers in Microbiology, 3, 2012. [45] Xiao Feng, Daniel S. Park, Cassondra Walker, A. Townsend Peterson, Cory Merow, and Monica Papes¸. A checklist for maximizing reproducibility of ecological niche models. Nature Ecology & Evolution, 3(10):1382–1395, 2019. [46] Emily Nicholson, Kate E. Watermeyer, Jessica A. Rowland, Chloe F. Sato, Simone L. Stevenson, Angela Andrade, Thomas M. Brooks, Neil D. Burgess, Su-Ting Cheng, Hedley S. Grantham, Samantha L. Hill, David A. Keith, Martine Maron, Daniel Metzke, Nicholas J. Murray, Cara R. Nelson, David Obura, Andy Plumptre, Andrew L. Skowno, and James E. M. Watson. Scientifc foundations for an ecosystem goal, milestones and indicators for the post-2020 global biodiversity framework. Nature Ecology & Evolution, 5(10):1338–1349, 2021. [47] Ammad Waheed Qazi, Zafeer Saqib, and Muhammad Zaman-ul Haq. Trends in species distribu- tion modelling in context of rare and endemic plants: a systematic review. Ecological Processes, 11(1):1–11, 2022. [48] Karoline Faust, Leo Lahti, Didier Gonze, Willem M. de Vos, and Jeroen Raes. Metagenomics meets time series analysis: unraveling microbial community dynamics. Current Opinion in Microbiology, 25:56–66, 2015. [49] Karoline Faust, Franziska Bauchinger, Béatrice Laroche, Sophie De Buyl, Leo Lahti, Alex D Washburne, Didier Gonze, and Stefanie Widder. Signatures of ecological processes in microbial community time series. Microbiome, 6(1):1–13, 2018. [50] Leo Lahti, Jarkko Salojärvi, Anne Salonen, Marten Scheffer, and Willem M De Vos. Tipping elements in the human intestinal ecosystem. Nature Communications, 5(1):4344, 2014. [51] Ophelia S. Venturelli, Alex V. Carr, Garth Fisher, Ryan H. Hsu, Rebecca Lau, Benjamin P. Bowen, Susan Hromada, Trent Northen, and Adam P. Arkin. Deciphering microbial interactions in synthetic human gut microbiome communities. Molecular Systems Biology, 14(6):e8157, 2018. [52] Peter J. Wangersky. Lotka-Volterra population models. Annual Review of Ecology and Systematics, 9(1):189–218, 1978. [53] Katharine Z. Coyte, Jonas Schluter, and Kevin R. Foster. The ecology of the microbiome: networks, competition, and stability. Science, 350(6261):663–666, 2015. [54] Christos Michalakelis, Charisios Christodoulos, Dimitrios Varoutas, and Thomas Sphicopou- los. Dynamic estimation of markets exhibiting a prey–predator behavior. Expert Systems with Applications, 39(9):7690–7700, 2012. [55] A.W. Wijeratne, Fengqi Yi, and Junjie Wei. Bifurcation analysis in the diffusive Lotka–Volterra system: an application to market economy. Chaos, Solitons & Fractals, 40(2):902–911, 2009. [56] Riet De Smet and Kathleen Marchal. Advantages and limitations of current network inference methods. Nature Reviews Microbiology, 8(10):717–729, 2010. [57] Mark Newman. Networks. Oxford university press, 2018. [58] Karoline Faust and Jeroen Raes. Microbial interactions: from networks to models. Nature Reviews Microbiology, 10(8):538–550, 2012. 81 Moein Khalighi [59] Timothy S Gardner, Charles R Cantor, and James J Collins. Construction of a genetic toggle switch in Escherichia coli. Nature, 403(6767):339–342, 2000. [60] Crawford S. Holling. Some characteristics of simple types of predation and parasitism. The canadian entomologist, 91(7):385–398, 1959. [61] Steven H. Strogatz. Nonlinear dynamics and chaos with Student Solutions manual: with applica- tions to Physics, Biology, Chemistry, and Engineering. CRC press, 2018. [62] Philip Hartman. A lemma in the theory of structural stability of differential equations. Proceedings of the American Mathematical Society, 11(4):610–620, 1960. [63] D. M. Grobman. Homeomorphisms of systems of differential equations. Doklady Akademii Nauk SSSR, 128:880–881, 1959. [64] Jean-Jacques Slotine and Weiping Li. Applied nonlinear control. pages 1123–1131, 1991. [65] Gary Wong, Wenjun Liu, Yingxia Liu, Boping Zhou, Yuhai Bi, and George F. Gao. MERS, SARS, and Ebola: the role of super-spreaders in infectious disease. Cell Host & Microbe, 18(4):398–401, 2015. [66] William Ogilvy Kermack and Anderson G McKendrick. A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 115(772):700–721, 1927. [67] Norman T. J. Bailey. The mathematical theory of infectious diseases and its applications. Charles Griffn & Company, 1975. [68] Herbert W. Hethcote. The mathematics of infectious diseases. SIAM Review, 42(4):599–653, 2000. [69] Simone Raponi, Zeinab Khalifa, Gabriele Oligeri, and Roberto Di Pietro. Fake news propagation: a review of epidemic models, datasets, and insights. ACM Transactions on the Web (TWEB), 16 (3):1–34, 2022. [70] Klaus Dietz. The estimation of the basic reproduction number for infectious diseases. Statistical Methods in Medical Research, 2(1):23–41, 1993. [71] Paul L. Delamater, Erica J. Street, Timothy F. Leslie, Y. Tony Yang, and Kathryn H. Jacobsen. Complexity of the basic reproduction number (R0). Emerging Infectious Diseases, 25(1):1, 2019. [72] P. van den Driessche and James Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences, 180(1): 29–48, 2002. [73] Joel C. Miller. A note on the derivation of epidemic fnal sizes. Bulletin of Mathematical Biology, 74(9):2125–2141, 2012. [74] Carlos Castillo-Chavez and Baojun Song. Dynamical models of tuberculosis and their applications. Mathematical Biosciences and Engineering, 1(2):361–404, 2004. [75] Ivan Area, Jorge Losada, Faïçal Ndaïrou, Juan J. Nieto, and Daniel D Tcheutia. Mathematical modeling of 2014 ebola outbreak. Mathematical Methods in the Applied Sciences, 40(17): 6114–6122, 2017. [76] Samuel Mwalili, Mark Kimathi, Viona Ojiambo, Duncan Gathungu, and Rachel Mbogo. Seir model for COVID-19 dynamics incorporating the environment and social distancing. BMC Research Notes, 13(1):1–5, 2020. [77] Robert M. Gray. Probability, random processes, and ergodic properties. Springer Science & Business Media, 2009. [78] Rosalind J. Allen and Bartlomiej Waclaw. Bacterial growth: a statistical physicist’s guide. Reports on Progress in Physics, 82(1):016601, nov 2018. [79] Ivanna Karina Olivares-Marin, Juan Carlos González-Hernández, Carlos Regalado-Gonzalez, and Luis Alberto Madrigal-Perez. Saccharomyces cerevisiae exponential growth kinetics in batch culture to analyze respiratory and fermentative metabolism. Journal of visualized experiments: JoVE, (139):58192, 2018. [80] Catherine L. Kelly, Lin Schwarzkopf, Iain J. Gordon, and Ben Hirsch. Population growth lags in introduced species. Ecology and Evolution, 11(9):4577–4587, 2021. 82 LIST OF REFERENCES [81] Camila Fontoura Acosta Ribeiro, Gislaine Greice de Oliveira Silva Silveira, Elizabete de Souza Candido, Marlon Henrique Cardoso, Cristiano Marcelo Espinola Carvalho, and Octavio Luiz Franco. Effects of antibiotic treatment on gut microbiota and how to overcome its negative impacts on human health. ACS Infectious Diseases, 6(10):2544–2559, 2020. [82] Alexander Pimenov, Theresa C. Kelly, Andrei Korobeinikov, Michael J. A. O’Callaghan, Alexei V. Pokrovskii, and Dmitrii Rachinskii. Memory effects in population dynamics: spread of infectious disease as a case study. Mathematical Modelling of Natural Phenomena, 7(3):204–226, 2012. [83] Cecilia Jernberg, Sonja Löfmark, Charlotta Edlund, and Janet K. Jansson. Long-term impacts of antibiotic exposure on the human intestinal microbiota. Microbiology, 156(11):3216–3223, 2010. [84] Lieselotte Vermeersch, Lloyd Cool, Anton Gorkovskiy, Karin Voordeckers, Tom Wenseleers, and Kevin J. Verstrepen. Do microbes have a memory? History-dependent behavior in the adaptation to variable environments. Frontiers in Microbiology, 13, 2022. [85] Antun Skanata and Edo Kussell. Ecological memory preserves phage resistance mechanisms in bacteria. Nature Communications, 12(1):6817, 2021. [86] Martin Krsak, Brian L. Harry, Brent E. Palmer, and Carlos Franco-Paredes. Postinfectious immunity after COVID-19 and vaccination against SARS-CoV-2. Viral Immunology, 34(8): 504–509, 2021. PMID: 34227891. [87] Angel N. Desai and Maimuna S. Majumder. What is herd immunity? JAMA, 324(20):2113–2113, 2020. [88] Laura Pérez-Alós, Cecilie Bo Hansen, Jose Juan Almagro Armenteros, Johannes Roth Madsen, Line Dam Heftdal, Rasmus Bo Hasselbalch, Mia Marie Pries-Heje, Rafael Bayarri-Olmos, Ida Jarlhelt, Sebastian Rask Hamm, et al. Previous immunity shapes immune responses to SARS- CoV-2 booster vaccination and Omicron breakthrough infection risk. Nature Communications, 14 (1):5624, 2023. [89] Prasanna Jagannathan and Taia T. Wang. Immunity after SARS-CoV-2 infections. Nature Immunology, 22(5):539–540, 2021. [90] Chadi M. Saad-Roy, Caroline E. Wagner, Rachel E. Baker, Sinead E. Morris, Jeremy Farrar, Andrea L. Graham, Simon A. Levin, Michael J. Mina, C. Jessica E. Metcalf, and Bryan T. Grenfell. Immune life history, vaccination, and the dynamics of SARS-CoV-2 over the next 5 years. Science, 370(6518):811–818, 2020. [91] Jin-Liang Wang and Hui-Feng Li. Memory-dependent derivative versus fractional derivative (I): difference in temporal modeling. Journal of Computational and Applied Mathematics, 384: 112923, 2021. [92] Wei Min, Guobin Luo, Binny J. Cherayil, S. C. Kou, and X. Sunney Xie. Observation of a power-law memory kernel for fuctuations within a single protein molecule. Phys. Rev. Lett., 94: 198302, May 2005. [93] Kenneth L. Cooke and Zvi Grossman. Discrete delay, distributed delay and stability switches. Journal of Mathematical Analysis and Applications, 86(2):592–627, 1982. [94] Yuguang Yang, Kevin R. Foster, Katharine Z. Coyte, and Aming Li. Time delays modulate the stability of complex ecosystems. Nature Ecology & Evolution, 7(10):1610–1619, 2023. [95] Sansit Patnaik, John P. Hollkamp, and Fabio Semperlotti. Applications of variable-order fractional operators: a review. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 476(2234):20190498, 2020. [96] Aaron Clauset, Cosma Rohilla Shalizi, and M. E. J. Newman. Power-law distributions in empirical data. SIAM Review, 51(4):661–703, 2009. [97] MEJ Newman. Power laws, Pareto distributions and Zipf’s law. Contemporary Physics, 46(5): 323–351, 2005. [98] Didier Sornette. Critical phenomena in natural sciences: chaos, fractals, selforganization and disorder: concepts and tools. Springer Science & Business Media, 2006. [99] Stephen J. Hardiman, Nicolas Bercot, and Jean-Philippe Bouchaud. Critical refexivity in fnancial markets: a hawkes process analysis. The European Physical Journal B, 86:1–9, 2013. 83 Moein Khalighi [100] Ralf Metzler and Joseph Klafter. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Physics Reports, 339(1):1–77, 2000. [101] Kai Diethelm, Neville J. Ford, and Alan D. Freed. A predictor-corrector approach for the numerical solution of fractional differential equations. Nonlinear Dynamics, 29(1-4):3–22, 2002. [102] Concepción A. Monje, YangQuan Chen, Blas M. Vinagre, Dingyu Xue, and Vicente Feliu-Batlle. Fractional-order systems and controls: fundamentals and applications. Springer Science & Business Media, 2010. [103] Yan Zhou, Hongxing Wang, and Heng Liu. Generalized function projective synchronization of incommensurate fractional-order chaotic systems with inputs saturation. International Journal of Fuzzy Systems, 21:823–836, 2019. [104] Kotadai Zourmba, Alhadji Abba Oumate, Gambo Gambo, J. Effa Effa, and Aliduo Mohamadou. Chaos in the incommensurate fractional order system and circuit simulations. International Journal of Dynamics and Control, 7:94–111, 2019. [105] Milad Shahvali, Mohammad-Bagher Naghibi-Sistani, and Hamidreza Modares. Distributed consensus control for a network of incommensurate fractional-order systems. IEEE Control Systems Letters, 3(2):481–486, 2019. [106] Mohammad Saleh Tavazoei and Mohammad Haeri. Chaotic attractors in incommensurate frac- tional order systems. Physica D: Nonlinear Phenomena, 237(20):2628–2637, 2008. [107] Keith Oldham and Jerome Spanier. The fractional calculus theory and applications of differentia- tion and integration to arbitrary order. Elsevier, 1974. [108] Mohammad Amirian Matlob and Yousef Jamali. The concepts and applications of fractional order differential calculus in modeling of viscoelastic systems: A primer. Critical Reviews™ in Biomedical Engineering, 47(4), 2019. [109] Vasily E. Tarasov. On history of mathematical economics: Application of fractional calculus. Mathematics, 7(6):509, 2019. [110] Hadiseh Safdari, Milad Zare Kamali, Amirhossein Shirazi, Moein Khalighi, Gholamreza Ja- fari, and Marcel Ausloos. Fractional dynamics of network growth constrained by aging node interactions. PLOS ONE, 11(5):1–13, 2016. [111] Edmundo Capelas de Oliveira and José António Tenreiro Machado. A review of defnitions for fractional derivatives and integral. Mathematical Problems in Engineering, 2014:238459, 2014. [112] Kai Diethelm, Roberto Garrappa, Andrea Giusti, and Martin Stynes. Why fractional derivatives with nonsingular kernels should not be used. Fractional Calculus and Applied Analysis, 23(3): 610–634, 2020. [113] Valentina V. Tarasova and Vasily E. Tarasov. Concept of dynamic memory in economics. Com- munications in Nonlinear Science and Numerical Simulation, 55:127–145, 2018. [114] Dazhi Zhao and Maokang Luo. Representations of acting processes and memory effects: General fractional derivative and its application to theory of heat conduction with fnite wave speeds. Applied Mathematics and Computation, 346:531–544, 2019. [115] Denis Matignon. Stability results for fractional differential equations with applications to control processing. In Computational Engineering in Systems Applications, pages 963–968, 1996. [116] Moein Khalighi, Leila Eftekhari, Soleiman Hosseinpour, and Leo Lahti. Three-species Lotka- Volterra model with respect to Caputo and Caputo-Fabrizio fractional operators. Symmetry, 13(3): 368, 2021. [117] Kai Diethelm, Neville J. Ford, and Alan D. Freed. Detailed error analysis for a fractional Adams method. Numerical Algorithms, 36(1):31–52, 2004. [118] Roberto Garrappa. On linear stability of predictor–corrector algorithms for fractional differential equations. International Journal of Computer Mathematics, 87(10):2281–2290, 2010. [119] Roberto Garrappa. Numerical solution of fractional differential equations: A survey and a software tutorial. Mathematics, 6(2):16, 2018. ISSN 2227-7390. [120] Zhuo Li, Lu Liu, Sina Dehghan, YangQuan Chen, and Dingyü Xue. A review and evaluation of numerical tools for fractional calculus and fractional order controls. International Journal of Control, 90(6):1165–1181, 2017. 84 LIST OF REFERENCES [121] Kai Diethelm and Neville J. Ford. Analysis of fractional differential equations. Journal of Mathematical Analysis and Applications, 265(2):229–248, 2002. [122] Kai Diethelm, Roberto Garrappa, and Martin Stynes. Good (and not so good) practices in computational methods for fractional calculus. Mathematics, 8(3), 2020. ISSN 2227-7390. [123] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B. Shah. Julia: A fresh approach to numerical computing. SIAM Review, 59(1):65–98, 2017. [124] Hadi Jahanshahi, Jesus M. Munoz-Pacheco, Stelios Bekiros, and Naif D. Alotaibi. A fractional- order SIRD model with time-dependent memory indexes for encompassing the multi-fractional characteristics of the COVID-19. Chaos, Solitons & Fractals, 143:110632, 2021. [125] Hong Ge, Kai Xu, and Zoubin Ghahramani. Turing: a language for fexible probabilistic inference. In International Conference on Artifcial Intelligence and Statistics, 2018, 9-11 April 2018, Playa Blanca, Lanzarote, Canary Islands, Spain, pages 1682–1690, 2018. [126] Tim Besard, Christophe Foket, and Bjorn De Sutter. Effective extensible programming: unleashing Julia on GPUs. IEEE Transactions on Parallel and Distributed Systems, 30(4):827–841, 2019. [127] Lawrence A. David, Arne C. Materna, Jonathan Friedman, Maria I. Campos-Baptista, Matthew C. Blackburn, Allison Perrotta, Susan E. Erdman, and Eric J. Alm. Host lifestyle affects human microbiota on daily timescales. Genome biology, 15:1–15, 2014. [128] Martin Rosvall, Alcides V. Esquivel, Andrea Lancichinetti, Jevin D. West, and Renaud Lambiotte. Memory in network fows and its effects on spreading dynamics and community detection. Nature Communications, 5(1):4630, 2014. [129] Oliver E. Williams, Lucas Lacasa, Ana P. Millán, and Vito Latora. The shape of memory in temporal networks. Nature Communications, 13(1):499, 2022. [130] Robert Zwanzig. Nonequilibrium statistical mechanics. Oxford University Press, 2001. [131] Floris Takens. Detecting strange attractors in turbulence. Springer, 2006. [132] Shinya Maeyama and Tomo-Hiko Watanabe. Extracting and modeling the effects of small-scale fuctuations on large-scale fuctuations by Mori–Zwanzig projection operator method. Journal of the Physical Society of Japan, 89(2):024401, 2020. [133] George E. P. Box, Gwilym M. Jenkins, Gregory C. Reinsel, and Greta M. Ljung. Time series analysis: forecasting and control. John Wiley & Sons, 2015. [134] Yanlin Shi and Kin-Yip Ho. Long memory and regime switching: a simulation study on the Markov regime-switching ARFIMA model. Journal of Banking & Finance, 61:S189–S204, 2015. [135] Niklas Wulkow, Péter Koltai, and Christof Schütte. Memory-based reduced modelling and data-based estimation of opinion spreading. Journal of Nonlinear Science, 31:1–42, 2021. [136] Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identifcation of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 113(15):3932–3937, 2016. [137] Philipp Singer, Denis Helic, Behnam Taraghi, and Markus Strohmaier. Detecting memory and structure in human navigation patterns using Markov chain models of varying order. PLOS ONE, 9(7):1–21, 2014. [138] Renaud Lambiotte, Vsevolod Salnikov, and Martin Rosvall. Effect of memory on the dynamics of random walks on networks. Journal of Complex Networks, 3(2):177–188, 2014. [139] Tiago P Peixoto and Martin Rosvall. Modelling sequences and temporal networks with dynamic community structures. Nature Communications, 8(1):582, 2017. [140] Takayuki Hiraoka and Hang-Hyun Jo. Correlated bursts in temporal networks slow down spreading. Scientifc Reports, 8(1):15321, 2018. [141] Oliver E Williams, Fabrizio Lillo, and Vito Latora. Effects of memory on spreading processes in non-Markovian temporal networks. New Journal of Physics, 21(4):043028, 2019. [142] Antonio Russo, Miguel A. Durán-Olivencia, Ioannis G. Kevrekidis, and Serafm Kalliadasis. Machine learning memory kernels as closure for non-Markovian stochastic processes. IEEE Transactions on Neural Networks and Learning Systems, pages 1–13, 2022. [143] Manisha Joshi, Savita Bhosale, and Vishwesh A Vyawahare. A survey of fractional calculus applications in artifcial neural networks. Artifcial Intelligence Review, pages 1–54, 2023. 85 [144] Emrah ¸ Power-law tail in lag time distribution underlies bacterialSims¸ek and Minsu Kim. persistence. PNAS, 116(36):17635–17640, 2019. [145] Jonas Neipel, Jonathan Bauermann, Stefano Bo, Tyler Harmon, and Frank Jülicher. Power-law population heterogeneity governs epidemic waves. PLOS ONE, 15(10):1–27, 2020. [146] Géza Ódor. Nonuniversal power-law dynamics of susceptible infected recovered models on hierarchical modular networks. Physical Review E, 103:062112, 2021. [147] Carles Falcó and Álvaro Corral. Finite-time scaling for epidemic processes with power-law superspreading events. Physical Review E, 105(6):064122, 2022. [148] Cesar Manchein, Eduardo L. Brugnago, Rafael M. da Silva, Carlos F. O. Mendes, and Marcus W. Beims. Strong correlations between power-law growth of COVID-19 in four continents and the ineffciency of soft quarantine strategies. Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(4):041102, 2020. [149] Kai Diethelm. A fractional calculus based model for the simulation of an outbreak of dengue fever. Nonlinear Dynamics, 71:613–619, 2013. [150] Nadjette Debbouche, Adel Ouannas, Iqbal M. Batiha, and Giuseppe Grassi. Chaotic dynamics in a novel COVID-19 pandemic model described by commensurate and incommensurate fractional- order derivatives. Nonlinear Dynamics, pages 1–13, 2021. [151] Caterina Balzotti, Mirko D’Ovidio, Anna Chiara Lai, and Paola Loreti. Effects of fractional derivatives with different orders in SIS epidemic models. Computation, 9(8), 2021. ISSN 2079-3197. [152] Kai Diethelm, Stefan Siegmund, and HT Tuan. Asymptotic behavior of solutions of linear multi- order fractional differential systems. Fractional calculus and applied analysis, 20(5):1165–1195, 2017. [153] Kai Diethelm, Ha Duc Thai, and Hoang The Tuan. Asymptotic behaviour of solutions to non- commensurate fractional-order planar systems. Fractional Calculus and Applied Analysis, 25(4): 1324–1360, 2022. [154] Mehmet Ali Balci. Fractional virus epidemic model on fnancial networks. Open Mathematics, 14(1):1074–1086, 2016. [155] Patrick Kofod Mogensen, John Myles White, Asbjørn Nilsen Riseth, Tim Holy, Miles Lubin, Christof Stocker, Andreas Noack, Antoine Levitt, Christoph Ortner, Blake Johnson, Dahua Lin, Kristoffer Carlsson, Yichao Yu, Christopher Rackauckas, Josua Grawitter, Alex Williams, Alexey Stukalov, Ben Kuhn, Benoît Legat, Jeffrey Regier, cossio, Michael Creel, Ron Rock, Thomas R. Covert, Benoit Pasquier, Takafumi Arakaki, Andrew Clausen, and Arno Strouwen. Julianlsolvers/optim.jl: Release 1.7.3, 2022. [156] Ensheng Dong, Hongru Du, and Lauren Gardner. An interactive web-based dashboard to track COVID-19 in real time. The Lancet Infectious Diseases, 20(5):533–534, 2020. ISSN 1473-3099. [157] Hamed Al-Sulami, Moustafa El-Shahed, Juan J. Nieto, and Wafa Shammakh. On fractional order dengue epidemic model. Mathematical Problems in Engineering, 2014:1–6, 2014. [158] Yu-Ming Chu, Saima Rashid, Ahmet Ocak Akdemir, Aasma Khalid, Dumitru Baleanu, Bushra R. Al-Sinan, and O.A.I. Elzibar. Predictive dynamical modeling and stability of the equilibria in a discrete fractional difference COVID-19 epidemic model. Results in Physics, 49:106467, 2023. [159] Soon-Mo Jung. Hyers-Ulam-Rassias stability of functional equations in nonlinear analysis, volume 48. Springer Science & Business Media, 2011. ISBN 978-952-02-0344-3 (PRINT) ISBN 978-952-02-0345-0 (PDF) ISSN 2736-9390 (PRINT) ISSN 2736-9684 (ONLINE) Pa in os al am a, Tu rk u, F in la nd 2 02 5