Timescale separation and models of symbiosis: state space reduction, multiple attractors and initialization

Abstract Dynamic Energy Budget models relate whole organism processes such as growth, reproduction and mortality to suborganismal metabolic processes. Much of their potential derives from extensions of the formalism to describe the exchange of metabolic products between organisms or organs within a single organism, for example the mutualism between corals and their symbionts. Without model simplification, such models are at risk of becoming parameter-rich and hence impractical. One natural simplification is to assume that some metabolic processes act on ‘fast’ timescales relative to others. A common strategy for formulating such models is to assume that ‘fast’ processes equilibrate immediately, while ‘slow’ processes are described by ordinary differential equations. This strategy can bring a subtlety with it. What if there are multiple, interdependent fast processes that have multiple equilibria, so that additional information is needed to unambiguously specify the model dynamics? This situation can easily arise in contexts where an organism or community can persist in a ‘healthy’ or an ‘unhealthy’ state with abrupt transitions between states possible. To approach this issue, we offer the following: (a) a method to unambiguously complete implicitly defined models by adding hypothetical ‘fast’ state variables; (b) an approach for minimizing the number of additional state variables in such models, which can simplify the numerical analysis and give insights into the model dynamics; and (c) some implications of the new approach that are of practical importance for model dynamics, e.g. on the bistability of flux dynamics and the effect of different initialization choices on model outcomes. To demonstrate those principles, we use a simplified model for root-shoot dynamics of plants and a related model for the interactions between corals and endosymbiotic algae that describes coral bleaching and recovery.


Introduction
Metabolic models are used to describe physiological processes in living organisms. They range in mathematical and computational complexity from high-dimensional, parameter-rich, systems biology models used in biomedical applications to low-dimensional models with many fewer parameters exemplified by models based on Dynamic Energy Budget (DEB) theory (Kooijman, 2010). Some of Kooijman's pioneering work in this field was motivated by challenges in environmental management, notably in ecotoxicology (Kooijman & Metz, 1984), with his first book subtitled Theory and Applications in Ecotoxicology (Kooijman, 1993). This focus remains important, e.g. (Murphy et al., 2018), and subsequent research has greatly widened the scope of applications to many areas of conservation physiology, notably changes in a species' ecological niche in response to environmental change (Kearney & Porter, 2020). Lavaud et al. (2021) offer an overview of such applications.
DEB models have been used in innovative research for nearly two decades to describe networks for the exchange of metabolic products between interacting organisms (Kooi et al., 2004, Kooijman et al., 2003, but this ecologically important area is still underdeveloped. The conceptual foundation within DEB theory for such models is the work by Kooijman (2001), which sketches links involving excretion fluxes from animal host and algal symbiont in corals and interactions between root and shoot in a plant. Subsequent models for the coral symbiosis and the interaction of roots and shoots in plants have been developed in Muller et al. (2009), Cunning et al. (2017), Ledder et al. (2020) and Schouten et al. (2020). These models can be applied to describe responses to environmental stress, e.g. coral bleaching (Cunning et al., 2017, Eynaud et al., 2011. Mathematically, several of the cited models take the form of differential equations that involve a small number of state variables. For instance, the models in (Cunning et al., 2017) and Ledder et al. (2020) describe two 'players' who exchange metabolic products that are surplus to their own needs. This interaction is described with a set of algebraic equations that characterize the flows of elemental matter within an interaction network; see schematic Fig. 1. The algebraic equations arise from the assumption that reactions and translocations within the network occur on a much faster timescale than changes in the state variables for the biomass of the interacting organisms. Therefore, the reactions and translocations can be assumed to equilibrate virtually immediately. This timescale separation allows for simpler model formulation with fewer parameters and can greatly improve numerical speed for simulations (a concern when using computer-intensive methods for parameter estimation).
Even conceptually simple networks typically involve many interconnected fluxes. As a consequence, such networks can have multiple equilibrium states. In those cases, unambiguous tracking of only slow variables is not possible without at minimum additional information about the state of the fast dynamic network. An example for such a situation is the earlier mentioned model for the growth of the root and the shoot of a plant by Ledder et al. (2020). This model does not only need information on the biomass of the two organelles, but it also needs to track the current rates at which root and shoot share nitrogen and carbon, respectively. The initialization of these rates can determine whether the system starts with either the root growing or the shoot growing. A similar example is the model for the interactions between a coral host and its endosymbiotic algae by Cunning et al. (2017). In this model, the initial fluxes can not only influence the initial growth phase but also determine whether the system reaches a healthy state with a high symbiont density or an unhealthy state where the coral loses its symbiont and starves eventually (coral bleaching). The work by Cunning et al. (2017) demonstrates how the issue of multiple fast equilibria in DEB models has been approached previously with a discrete-time scheme. In this scheme, the ambiguity of the model definition is resolved by adding state variables for the fast variables and using the state of the system at the preceding time step in the update step.
In the present paper, we propose a new approach for the issue of tracking the fast network. The new scheme adds an explicit continuous-time representation of the fast dynamics by introducing additional state variables. Those additional state variables buffer some connections in the network and in this way complete the definition of the fast dynamics. Among the advantages of this new method is that it relies on continuous time and thus allows simulating models with standard solvers for ordinary differential equations.
While this new approach can be used generally as a numerical approximation of unknown fast dynamics, we show how under certain circumstances it can be interpreted mechanistically. We further connect this new approach to previous work by showing that the continuous-time approach is a generalization of the discrete-time scheme used before.
Our approach for completing implicitly defined models is based on additional state variables (buffers) that interrupt cycles in the metabolic network. Naturally, it is possible to track all auxiliary variables, but often it is possible to only track some of them and still fully define the model. This reduces the dimensionality of the system, which can simplify model analysis. We will show that graph theory can be used to decide which variables to track and to find the minimal number of additional state variables that are needed to fully describe such systems.
The paper is structured as follows: Section 2 describes the key steps in our approach to timescale separation. In Sections 3 and 4 we use the new methods in analyses of the simple root-shoot model from Ledder et al. (2020) and the more complex model for the dynamics of corals with their endosymbiotic algae partners from Cunning et al. (2017). In Section 5 we discuss our findings and offer an outlook on wider areas of application, including modeling complex metabolic networks that are relevant to pressing environmental challenges.

Timescale separation in metabolic models
Cyclic and acyclic metabolic networks.
The Introduction highlighted dynamic bioenergetic models in which the state X of interacting organisms (or of organs within a single organism) is described by a set of differential equations coupled to an auxiliary network of variables, Z. A general form for the differential equation for X is The auxiliary variables Z most commonly represent fluxes. To minimize abstraction, we refer to them as fluxes in the remainder of this section. Although the fluxes Z may depend on the state variables, they are typically assumed to change on a sufficiently fast timescale that they can be replaced at any time by their equilibrium values (given the current state X). We shall on occasions refer to these value as 'fast equilibrium'. A model is unambiguously defined in this way when the auxiliary variables can be expressed as simple functions of the state variables, An example for such a system is shown on the left side of Fig. 2. In this example, knowledge of the state variables permits direct evaluation of Z 1 , which in turn allows evaluation of Z 2 and Z 3 .
Such sequential evaluation is not necessarily possible if the fluxes in the network depend on each other in a cyclic fashion (right side of Fig. 2). In this case, the fast equilibrium states are defined implicitly as solutions of an equation (almost always non-linear) of the form With an implicitly defined network, it is possible that multiple values of Z satisfy the equilibrium condition. In this case, additional assumptions on fast dynamics are needed to fully specify the model. This is not a mere mathematical nicety: the two case studies in subsequent sections will show that the initial state of a network will frequently have lasting impact on the long-term (slow) dynamics of the full system. In short, to unambiguously define our models, we need to specify fast dynamics.

Differential equations for fast and slow dynamics
To complete models with implicitly defined auxiliary networks, we set up hypothetical fast dynamics that mimic the (generally unknown) actual fast dynamics. Those hypothetical fast dynamics should satisfy two conditions. First, at their equilibrium they should satisfy the model specification Z = g( X, Z) as stated in equation (3). Second, each flux Z i should approach its 'target' g i ( X, Z). Arguably, the simplest solution for these dynamics is to let each flux decay exponentially towards its target. We therefore propose to replace the implicit equation for Z with the dynamical system where the set of parameters {λ i } determines how quickly the fluxes approach the fast equilibrium. When the values for the λ i are high, the fluxes typically track their targets tightly. The anticipated properties of our scheme for model completion are supported by a mathematical result-Tikhonov's theorem (Klonowski, 1983, Tikhonov, 1952. Expressed informally, the theorem states that in systems with fast and slow variables the fast dynamics typically track an isolated fast equilibrium after some initial transient (i.e. the auxiliary variables are sticking close to their previous state). Appendix A gives a precise statement of the theorem together with three conditions for it to hold. Appendix B outlines in detail the reasoning that connects the theorem to our differential equation scheme and also shows that the current approach is asymptotically a generalized form of the discrete time method used in a previous work (Cunning et al., 2017).
Our choice of equations for fast dynamics was motivated by mathematical requirements and took no account of the physiological/biochemical nature of the networks. There is no reason to believe that they represent the 'real' underlying fast dynamics that are generally unknown. In some cases, however, the differential equations we introduced can plausibly represent original fast dynamics. The root-shoot model is such a case. We work out the details for this example in the section for this model. Generally, our numerical scheme emerges when the nodes of the auxiliary network are assumed to act as 'buffers' that delay inputs within the network. Those buffers then represent the state of the network. Assume the following fast network. The state of each part i of a network is described by a buffer E i that feeds into the network at a rate λ i E i . The input of the buffers is given by g i ( X, E) (where is a diagonal matrix with the λ i on its diagonal). This gives us the differential equation for the buffers: Letting Z = E we find This is equivalent to the implementation described earlier, equation (4).

Reducing the dimensionality of the fast network
The natural implementation of our approach is to track all auxiliary variables by additional state variables. This increases the dimensionality of the system by the number of auxiliary variables. However, often the auxiliary network can be already fully described with less state variables. This can be useful for analytical and numerical model analysis, particu-larly because auxiliary variables need to be initialized and this initialization then in turn can determine the long-term (slow) dynamics of the system. To reduce the dimensionality of the auxiliary fast dynamics, we find combinations of auxiliary variables through which the rest of the auxiliary network can be calculated.
Biologically, this state space reduction corresponds to implementing a time scaling argument within the fast dynamics-some processes are on a fast timescale and some on a very fast timescale. The fast processes are tracked with differential equations, and the very fast processes are assumed to equilibrate instantaneously and are calculated directly.
When there is no biological intuition on the relative speed at which the different fast processes act, graph theory can be used to find low-dimensional representations of the network. For this purpose, we represent the interactions of the model fluxes with a directed graph as shown in Fig. 2. The nodes of the graph represent the fluxes, and the arrows indicate the connections between the fluxes. An arrow from flux Z i to flux Z j indicates that Z i appears in the function for calculating Z j . When this graph is acyclic, the model specification readily defines the fluxes uniquely and no additional state variables are needed to simulate the system (left side of Fig. 2). On the other hand, when this graph is cyclic, the model equations may be satisfied for different combinations of the fluxes (right side of Fig. 2). Tracking a flux by an additional state variable removes this flux from the network of algebraic equations because this flux does not equilibrate instantaneously anymore. This interrupts all cycles in which the flux is involved. Thus, in order to fully define the system, we need to find combinations of fluxes that, when expressed by state variables, interrupt all cycles in the graph. In the example on the right side of Fig. 2, tracking any of the fluxes with a state variable is enough to interrupt all cycles and uniquely define the system. The remaining fluxes can be calculated one by one from the value of this flux (and the slow dynamics, which have not been specified in the figure). The process is formally set out in Appendix C.
The Mathematica code in the supplementary material demonstrates how to identify whether a combination of fast state variables leads to an acyclic graph and how to find the minimal number of buffers that renders the fast dynamic network acyclic and defines the model uniquely.
In the following sections, we will illustrate those concepts with examples.  shares surplus nitrogen not used with the shoot, and the shoot shares surplus carbon with the root.

Model formulation
A detailed diagram of the model is shown in Fig. 3b. The state variables, fluxes and parameters are summarized in Table 1.
The equations for the root-shoot model are as follows. The assimilation of nitrogen by the root and the assimilation of carbohydrates (short:carbon) by the shoot are where the parameters α N and α C describe the rates at which root and shoot produce nitrogen and carbon, respectively.
The version of the model we describe here assumes that growth of each organ is limited by the more limiting resource, i.e. growth is governed by a minimal synthesizing unit (SU) (Ledder et al., 2020). Growth of root and shoot is given by where the parameters η −1 R and η −1 S describe the N:C ratio in root and shoot biomass and ρ N and ρ C are the rates at which nitrogen and carbon are shared, respectively. In the original model formulation, the nutrient sharing rates are defined by the difference of assimilation and use, ρ N = g ρ N and ρ C = g ρ C , where

Multiple solutions of flux combinations
The model as defined above does not always define all fluxes uniquely. Given values for R and S, the model equations can be satisfied for different combinations of the surplus sharing fluxes ρ N and ρ C (the auxiliary network). Figure 4 shows solutions for the surplus sharing fluxes as functions of the root-shoot ratio R/S. An analogous plot is shown in Fig. 3 of Ledder et al. (2020).
The reason multiple flux combinations can satisfy the original model equations is that the model fluxes for surplus sharing and growth form a closed loop that feeds back to itself, as shown in Fig. 3b. The loop can be described as follows. When the root grows, it shares less nitrogen, thus the shoot grows less (or not at all as in our example) and shares more carbon, reinforcing continuous growth of the root. In the other way around, when the shoot grows, it shares less carbon, thus the root grows less (or not at all as in our example) and shares more nitrogen, and shoot growth is reinforced. In the simplest form, this loop can be seen by expressing the two fluxes for surplus sharing as functions of each other by plugging in the equations for the growth rates of root and shoot, This loop is schematized in Fig. 5  Adding state variables to the model to obtain a uniquely defined system and simulate the model As shown in Fig. 4, the network of auxiliary variables (the rates at which nitrogen and carbon are shared) can have up to three solutions for a given root-to-shoot ratio. Additional model assumptions are needed to simulate the system. Here we show how to unambiguously define the system using the approach outlined in Section 2. The method bases on adding state variables for the surplus sharing rates to the model equations so that the cycles within the auxiliary network are interrupted.
To derive the method mechanistically for the present example, we add buffers that gather the fluxes of nitrogen and carbon shared by the root and the shoot.
Nitrogen that is shared by the root is first captured in a buffer E ρ N . The input rate to the buffer is the rate at which the root produces surplus nitrogen, g ρ N . Additionally, we assume that the buffer empties at a rate λ ρ N E ρ N into the growth SU for the shoot. This means the buffer changes according to When the buffer is at equilibrium, dE ρ N /dt = 0, and therefore This rate ρ N at which nitrogen is shared can be expressed directly by the differential equation The rate nitrogen is shared ρ N obviously follows its equilibrium g ρ N (assuming the other state variables are constant). The speed at which the tracks its equilibrium is given by the parameter λ ρ N .
In the same way, we can introduce a differential equation for the rate at which carbon is shared, ρ C , The model is now fully defined as a system of four differential equations that track the states of the biomass of root and shoot R and S, and the rate at which nitrogen and carbon is shared, ρ N and ρ C . This choice of additional state variables is shown in the bottom row of Fig. 5.
A phase plane for the fast buffer dynamics is shown in Fig. 6. The simulations assume that the biomass of root and shoot are kept constant (assuming buffer dynamics are much faster than biomass dynamics). The plot shows the bistability of the (fast) system: when starting with a high rate of nitrogen shared ρ N , the system converges to an equilibrium with high ρ N ; when starting with a high rate of carbon shared ρ C , the system converges to an equilibrium with high ρ C . These two equilibria correspond to states where, respectively, only the root or only the shoot is growing.  Simulations of the long-term (slow) dynamics with two buffers are shown in Fig. 7. The initial values for root and shoot biomass correspond to the (fixed) values for those state variables in the phase plane for the fast dynamics in Fig. 6. As the slow dynamics progress, the phase plane of the fast dynamics would change. The long-term simulations are started with two different choices for the initial values for the fast variables. The simulations demonstrate that the initialization of the fast variables can not only influence a short initial transient of the fast variables, but also influence long-term behavior by determining whether the system starts either with only the root growing or with only the shoot growing.

Reducing the number of additional state variables
The bottom row of Fig. 5 shows how the root-shoot model can be uniquely defined with a single buffer for the nitrogen shared by the root or the carbon shared by the shoot. The loops in the model are interrupted by either adding a differential equation for ρ N or for ρ C .
We can create an adjacency matrix that shows which fluxes are directly connected (each row showing how the flux indicated on the left depends on the other fluxes). The matrix notation can be used to apply graph theory algorithms to identify whether a graph for the model fluxes is cyclic (and thus does not define the model uniquely) and which choices for adding buffers to the model result in an acyclic graph (and thus defines the model uniquely).
Without additional state variables, the graph for the fast dynamics is cyclic (top row of Fig. 5). Its adjacency matrix is Adding state variables for the nitrogen and/or carbon shared makes the graph acyclic. Adding state variables for both surplus sharing fluxes (middle row of

The coral symbiont-host model
The model by Cunning et al. (2017) describes the dynamics of corals and their symbiotic algae. Similarly to the rootshoot model, the coral host shares nitrogen and the symbiont shares fixed carbon with their partners. A sketch of the model is shown in Fig. 8. The model details are stated in Appendix D.

Fast and slow dynamics
Many of the auxiliary variables (fluxes) in this model depend on other fluxes in the model. As in the root-shoot model, this makes it possible that multiple values for the fluxes satisfy the model equations for given S and H values. Figure 10 shows solutions for the fluxes as functions of the symbiont-host ratio S/H, together with the corresponding growth rates of the symbiont and the host. The different solutions correspond to a 'functional' symbiosis (blue), a 'dysfunctional' symbiosis (orange) and an unstable equilibrium between those two states (green). In the functional symbiosis, photosynthesis, carbon shared by the symbiont and CO 2 available to the symbiont are high, while nitrogen shared by the host and photodamage of the symbiont are low. When the S/H ratio is at its corresponding equilibrium (when dS/dt/S − dH/dt/H = 0), this state has a positive growth rate of symbiont and hosts (dS/dt/S and dH/dt/H). The dysfunctional symbiosis describes the opposite case, including negative growth rates at the corresponding equilibrium S/H ratio.

Adding state variables to the model to obtain a uniquely defined system and simulate the model
The model can be completed by adding state variables for all 12 fluxes involved in loops (the variables in the network shown in Fig. 9). Supplementary Fig. 1 shows the minimal combinations of buffers needed to complete the system. As all code, the implementation for finding those combinations is available online. The plots in the supplementary material show that the system needs at least two buffers. All twobuffer combinations include j CP as one of the additional state variables, and either j HG , j SG , ρ C or ρ N as the second additional state variable. For the coral model, the fast state variables do not have an obvious mechanistic interpretation because of the complex interactions between the auxiliary variables. One could change the model formulation slightly to accommodate for buffers with clear mechanistic interpretations, as worked out for the root-shoot model. This could be done by choosing fluxes to buffer that represent transfer of substances (such as with our current choices the photosynthesis rate j CP , but not the symbiont growth rate j SG ). However, working out the details is out of scope for this paper. For the present analysis, we introduced the fast state variables just as a numerical method to track solutions of the auxiliary network.  as the initial values of the chosen buffers can affect the model dynamics and the final state to which the system converges. The left side of Fig. 11 represents an initialization tool that shows how the initial values of the fast state variables can influence the dynamics of those state variables. The figure demonstrates the bistability of the fast dynamics, whereby the photosynthesis rate j CP and the symbiont growth rate j SG have been chosen as additional state variables. The fast dynamics can converge to either a state with high or low photosynthesis rate. The right side of Fig. 11 shows how this bistability of the fast dynamics in turn can influence the slow dynamics of the system (the change of biomass of symbiont and host). Depending on the initial photosynthesis rate, the system converges either to a functional symbiosis with a high S/H ratio and positive growth rates or a dysfunctional symbiosis with a low S/H ratio and negative growth rates. Figure 12 Table 2 in Cunning et al. (2017), with additionally L = 30, N = 1.5 × 10 −6 and X = 0. of the fast dynamics and the consequences on the slow system.

Discussion
DEB theory has been underutilized by ecological and evolutionary theorists modeling species interactions. It offers models of intermediate complexity that take account of thermodynamic and evolutionary constraints on the flows of energy and elemental matter at both ecological and evolutionary timescales (Jusup et al., 2017). One possible reason for this neglect is that much 'mainstream' theory focuses on changes in qualitative properties of simplified dynamical systems. Full DEB models of species interactions are typically too equationand parameter-rich to permit insightful qualitative analysis. To bridge this gap, many studies reduce the number of state variables, e.g. Chapter 9 of Kooijman (2010). A few studies archive further simplification by appealing to timescale separation, e.g. Poggiale et al. (2020). In similar spirit, we have here analyzed the dynamics arising from timescale separation in symbiotic interactions that involve shared metabolic products.
DEB models with timescale separation are of particular interest for modeling species interactions because they maintain the detailed physiological dynamics of single-organism DEB models while keeping model complexity low by assuming that physiological processes proceed faster than ecological dynamics. In published DEB-inspired models, timescale separation has been implemented by setting up differential equations for state variables that are solved simultaneously with algebraic equations. These algebraic equations can, in principle, be derived from underlying fast biochemical pro- cesses that equilibrate virtually instantaneously, but the details for those fast processes are often not explicitly provided in the model description. If these algebraic equations have a unique solution, simulating and analyzing such models is straightforward because the algebraic equations, together with the differential equations, unambiguously define the model dynamics. However, when the algebraic equations have multiple (or no) solutions, this is not true and the model dynamics are not defined uniquely anymore. In the present work, we showed how such ambiguity can be resolved by hypothesizing explicit fast dynamics that approach the solutions of the algebraic equations in the model specification. We emphasized that this type of ambiguity in the model definition is not only a technical mathematical issue since the initialization of the fast variables can impact long-term dynamics of the model system.
When setting up bioenergetic models that involve networks of fast processes, one often realizes quickly that such models can easily become complex and difficult to analyze. To address this issue, we offer a recipe for reducing the complexity of the underlying fast networks. The recipe is based on choosing a few auxiliary variables to be tracked directly through differential equations and evaluating the remaining quantities explicitly from these state variables. Biologically, this corresponds to assuming that the fast dynamic processes themselves act on different timescales: fast dynamics that are described through state variables and very fast dynamics that are assumed to equilibrate immediately and can be calculated using the rest of the system. Reducing the dimensionality of the system in this way can simplify numerical and analytical analysis. This simplification is especially useful for more complex models, such as our example for the symbiosis between corals and  Figure 12: Simulation of the coral model, implemented with fast state variables for photosynthesis j CP and symbiont growth j SG . From left to right and top to bottom: symbiont biomass (log scale), host biomass (log scale), ratio of symbiont to host, photosynthesis and symbiont growth rate. Blue curves: starting with high photosynthesis, j CP (0) = j CPm , and low symbiont growth, j SG (0) = 0, leads to a functional symbiosis with a high S/H ratio and positive growth rates. Orange curves: starting with low photosynthesis, j CP (0) = 0, and high symbiont growth, j SG (0) = j SGm , leads to a dysfunctional symbiosis with a low S/H ratio and negative growth rates. Parameter values as in Fig. 11. Initial symbiont to host ratio is S(0)/H(0) = 0.15. endosymbiotic algae (Cunning et al., 2017), or even more for the extension of this model with multiple symbionts (Brown et al. 2022).
We demonstrate these approaches for two published models: one for the growth dynamics of the root and the shoot of a plant (Ledder et al., 2020) and one for the growth dynamics of corals and their endosymbiotic algae (Cunning et al., 2017). We show how these models can be simulated and how their analysis can be simplified by reducing their dimensionality. Both examples illustrate how the choice of initial conditions can affect a system's state through cyclic feedbacks. In the example for the root and the shoot of a plant, an initially high carbon transfer triggers the system to start with the root growing, while an initially high nitrogen transfer triggers the system to start with the shoot growing. Similarly, in the coral example, the rate at which the symbiont shares carbon with its host turns out to be a particularly influential node in the fast network. Starting the system with a low photosynthesis rate can lead to a negative feedback loop, involving reduced activity of the host's carbon concentration mechanism that provides CO 2 to the symbiont. The decreased CO 2 transfer in turn additionally reduces the rate at which the symbiont shares carbon with the host. Eventually, this cycle leads to a breakdown of the symbiosis. In this example, the bistability of the fast dynamics can affect the final state of the slow dynamics, i.e. the initialization of the fast network determines whether the system settles at a functional symbiosis characterized by a high symbiont to host ratio and positive growth rates, or whether it settles at a dysfunctional symbiosis characterized by a low symbiont to host ratio (bleached corals) and negative growth rates.
The bleaching of coral reefs is just one example of how understanding the breakdown of a symbiosis is critical for conservation. DEB models of ecological interactions, analyzed with our methods for simulating and simplifying the fast dynamics, could be useful in similar contexts. Such other examples where our approach could be used to study symbiosis include situations where abiotic changes driven by human activity (such as rising temperatures) can decouple the growth rates of hosts and their symbionts. This can impact animal growth (Greenspan et al., 2020), the digestive system of humans (Carding et al., 2015), the cycling of organic matter by the microbiome of marine macroalgae (Minich et al., 2018) and the drought resistance of plants that harbor ectomycorrhizal fungi (Sapes et al., 2021).
The scenarios described above can be all interpreted as bioenergetic systems. However, these are not the only type of systems where complex fast dynamics can have alternative equilibrium states. For example, alternative states of fast state variables have been reported in models that couple evolutionary and ecological dynamics (Brown & Akçay, 2019, Geritz et al., 2002, Lehtinen & Geritz, 2019, in strictly ecological models (Rinaldi & Muratori, 1992) and in neurobiological models (Kurikawa & Kaneko, 2021, Wernecke et al., 2018. In these examples, fast dynamics were either simulated directly, or the models were analyzed without simulating the full coupled system of fast and slow dynamics. These examples reflect that models from different areas could warrant approaches similar to those described in the present work, which allow reduction of the dimension of the state space and approximation of fast dynamics efficiently for simulations. We end as we started-by emphasizing the unrealized potential for new DEB theory involving the exchange and interaction of metabolic products. Challenging areas include the ecology and evolution of the microbiome and its interaction with animal hosts, the evolution of mixotrophy, the emergence of microscale microbial consortia in aquatic systemsand of course there are many more. These and other fast processes may often drive population and community responses to changing environments, the core issue for conservation biology. Systematizing methodology for timescale separation in the absence of empirical data at all timescales will be essential in such work; we see this paper as a contribution to this effort. All code for the models has been implemented in Wolfram Mathematica 12.2 (Wolfram Research, 2021). The code is freely available at https://github.com/ferdi-p/time-scaleseparation/tree/main/code.

Data availability
No new data were generated or analysed in support of this research.

Funding
This work was supported by the National Science Foundation through research grants EF-1921356 (to H.V.M. and R.M.N.) and EF-1921425 (to R.C.). Support was also provided by a Simons Foundation Early Career Award in Marine Microbial Ecology and Evolution (to H.V.M.). A.R.D. was supported by an NSF Graduate Research Fellowship.

A Tikhonov's theorem
Tikhonov's theorem (Klonowski, 1983, Tikhonov, 1952 states that under certain conditions the fast-slow system as long as the point ( X, t) does not leave a certain region D of the state space.
The conditions are as follows: 1. Z = φ( X, t) is an isolated solution of h( X, Z, t) = 0 whose points are locally (asymptotically) stable equilibrium points of the adjoined system ( X, t being considered here as parameters, the point ( X, t) belonging to a bounded open region D). 2. The initial conditions ( X 0 , Z 0 , t 0 ) are such that the solution Z(τ ) of the adjoined system 3. X, f are n-vectors, Z, h m-vectors and f , h satisfy suitable regularity assumptions.

B Details to differential equations for fast and slow dynamics
To apply this theorem to our type of system, we assume that the algebraic equation for Z can be defined through the equilibrium states of a differential equation on a fast timescale (even when those equations are not explicitly stated in the model description). The form of the differential equation is generally Here γ is a parameter that determines the speed of the system. Tikhonov's theorem loosely tells us that the slow variables converge for γ → ∞ to the form is an isolated solution of h( X, Z) = 0 whose points are locally (asymptotically) stable equilibrium points In the absence of more information, we have to add assumptions to specify the model dynamics. It seems generally reasonable to assume that each component Z i of the network is attracted to its 'target' g i ( X, Z), i.e. the sign of the time derivative h i is the same as the sign of the difference between target and state, Qualitatively, this behavior can be captured by linearization, i.e.
where c i > 0 are constants that determine the (relative) speed with which each part of the network approaches its target.
Choosing λ i = γ c i , this brings us to equation (4) Note that the different λ i are generally not specified by the model (and they can in principle also depend on the state of the system, i.e. which equilibrium is in the vicinity of the auxiliary network). When applying the method to a new model, it can be useful to try different large values for the λ i and check convergence.
It can be useful to vectorize the model. Let In absence of any other information, a simple choice can be to choose all λ i equal and reasonably big.
The scheme we propose is asymptotically a generalization of the discrete time scheme implemented in earlier work (Cunning et al., 2017), when all λ i = λ and the time step t = 1/λ is small. This can be seen by expressing equation (4) as and solving for Z(t + t) so that Z(t + t) = g( X(t), Z(t)). (B.10)

C Dimensionality reduction
Divide the complete set of auxiliary variables into two sets, Z A and Z B , where ∪ refers to 'combining the two vectors and rearranging them in the original order') and then it is enough to set up dynamics for Z A and derive Z B from it, which completes the network. That is, Z A is described by a differential equation where A is a diagonal matrix that contains the numerical parameters λ i for how fast the different components of Z A approach their targets. The remainder of the network, Z B , can then be calculated directly: When the fast dynamics are at their equilibrium, any choice for which auxiliary variables to track as state variables satisfies the model specification. While in our experience the choices make little difference, they can principally differ in their dynamics. When setting up a new model, it can be useful to try different combinations of fast state variables and compare simulations.

D Coral model formulation
This section describes the details of the coral model by Cunning et al. (2017). The model tracks the biomass of symbiont S and host H. Additionally, it tracks various fluxes within and between the organisms, as shown in Table D1. The model parameters are shown in Table D2  that for the coral model the fluxes are measured per biomass (either of the host or the symbiont, depending on where the flux is originating from) while in the root-shoot model the fluxes describe the total values (i.e. already accounting for the biomass of the organ where the flux is originating from). Both models could be formulated in either way, but for consistency with the original papers we stick with the formulations used there. We also fixed two obvious typos in the original formulation of the coral model, changing slightly the equations for j HG and ρ N .
For convenience, we define notations for the functional responses used. For the uptake of a single substrate x the Michaelis-Menten-type SU is defined by For details, see Cunning et al. (2017).