Chaos in synthetic microbial communities

Predictability is a fundamental requirement in biological engineering. As we move to building coordinated multicellular systems, the potential for such systems to display chaotic behaviour becomes a concern. Therefore understanding which systems show chaos is an important design consideration. We developed a methodology to explore the potential for chaotic dynamics in small microbial communities governed by resource competition, intercellular communication and competitive bacteriocin interactions. Our model selection pipeline uses Approximate Bayesian Computation to first identify oscillatory behaviours as a route to finding chaotic behaviour. We have shown that we can expect to find chaotic states in relatively small synthetic microbial systems, understand the governing dynamics and provide insights into how to control such systems. This work is the first to query the existence of chaotic behaviour in synthetic microbial communities and has important ramifications for the fields of biotechnology, bioprocessing and synthetic biology.


Author summary
In chaotic systems, infinitesimally small differences in the initial conditions will become amplified over time, making forecasting and prediction of behaviour impossible. Although we know that chaos can be observed in the complex networks of natural ecosystems, the field of biotechnology is interested in designing and building new microbial communities and the presence of chaotic behaviour is unexplored. In this paper, we present a statistical pipeline that can tell us how, when and why chaos arises in small microbial communities. We apply this approach to study a set of communities involving quorum sensing systems and amensal interactions through antimicrobial peptides. Out of 4182 interaction networks in these three strain communities, we identify the networks that have the highest propensity to produce chaos. We then explore the levers we can pull to bring these networks in and out of chaotic regimes. Our work is the first to look at chaos in synthetic microbial communities and indicates that chaos is an important design consideration. Introduction chaotic attractors [20]. The possibility of chaos in synthetic microbial communities, to our knowledge, has not been previously considered. Standard competitive gLV models can produce chaotic behaviour in four species networks [11]. We use these previous findings to demonstrate and validate our methodology, before applying it to models that describe interpopulation interactions that are more specific to mechanisms found in synthetic microbial communities.

Development of a novel statistical approach to identifying chaotic regions in a multidimensional parameter space
Approximate Bayesian Computation Sequential Monte Carlo (ABC SMC) is a method that can be used for model selection and parameter inference in dynamical systems [21] (Algorithm 1). This flexible algorithm can also be used to tackle the design question, namely what model topologies and parameters are able to give rise to some specified target qualitative behaviour [22]. ABC SMC requires a distance function, describing how far away a simulation is from the objective behaviour. When searching for chaotic beahviour, we use the maximal Lyapunov exponent (λ 1 ) to create a distance function. We calculate λ 1 by initialising two nearby orbits and measuring their divergence or convergence over the course of a simulation (see Methods and Algorithm 1). λ 1 < 0 corresponds to linear stability, λ 1 = 0 corresponds to periodic oscillations, and λ 1 > 0 corresponds to chaos. While these rules hold true for infinite time, our simulations run for a finite time, meaning these boundary rules can be noisy. To identify chaos, we therefore define a threshold above 0 where we can be sure simulations have chaotic behaviour.
First, we demonstrate the use of ABC SMC in resolving a chaotic parameter distribution in a competitive gLV system. Competitive gLV equations are commonly used in ecological population modelling, and have similarly been used to model microbial communities [23]. They describe generic negative interactions between species that could represent competition for nutrients or amensal interactions. Competitive gLV systems take the form where N i is the size of a species population, r i is the growth rate, n is the number of species in the population and α the interaction matrix, describes the amensal interactions between pairs of species in the system. To simulate the chemostat environment, we set the diagonal as a dilution rate, D, which is the same for all species. The diagonal of α can also be thought of as defining the carrying capacity of each species. Vano et al. previously identified a chaotic attractor in this system using a brute-force parameter search [11]. Fig 1A shows the parameters identified, Fig 1B shows the resulting chaotic timeseries of the four species. We wanted to see if our ABC SMC methods could provide a posterior distribution for chaotic behaviour, capturing the findings of Vano et al. We identified a threshold of λ > 0.015 was sufficient for classifying chaotic behaviour and ran ABC SMC for this chaotic objective. We show the posterior of several parameters in Fig 1C. The black point corresponds to the parameters found by Vano et al, while the red scatter points correspond to chaotic behaviour we identified using ABC SMC. We can see that the interspecies interaction parameters, a 42 and a 43 , are constrained, indicating their importance for producing chaotic behaviour, given the prior parameter distributions. Conversely, the initial population of N 1 is not constrained, indicating the chaotic behaviours are robust to changing initial conditions. Similarly, r 3 , defining the growth rate of N 3 is not constrained. The full posterior parameter distribution is shown in S1 Fig. Mechanisms of interaction in microbial communities such as crossfeeding and toxin interactions would be subjected to time delays, accumulation of intermediate species and dynamic genetic regulation, contributing to non-linearity of these systems. gLV equations simplify these mechanisms and as such, are unable to capture chaotic behaviour with three species. In the next sections we move to studying more biochemically realistic systems.

Searching for chaos across synthetic microbial community models
In previous work we developed a model framework to describe QS regulated bacteriocin interactions in a three strain model space, and predicted topologies that form stable communities [17]. Here we use this same model space to investigate the existence of chaos in three strain synthetic microbial communities.  [11]. A Shows the parameters used in the chaotic attractor and an illustration of the interaction topology. B Time series of the chaotic attractor. C Posterior parameter distribution for chaotic objective, identified using ABC SMC (red) and the individual chaotic particle identified by Vano Fig 2A shows the pipeline we developed to search for chaos in synthetic three strain systems. Three strains, N 1 , N 2 , N 3 , optionally express bacteriocins, B 1 , B 2 , and B 3 under the control of optionally expressed QS molecules, A 1 , and A 2 . The QS molecules regulate expression of bacteriocins positively or negatively. Each strain can be optionally sensitive to a bacteriocin. The initial model space describes an enumeration of possible combinations of bacteriocin and QS systems, and forms the first uniform prior model space of 4182 models (Fig 2A(i)). Prior parameter distributions describe the range of characteristics for the different parts (Table 1). We expected the existence of chaos to be sparse in this three strain model space, and therefore computationally expensive to explore. Oscillations are a known route to chaos [1], therefore, in order to narrow down the search, we defined a novel set of three distances that are used to classify oscillatory behaviours. These were the period of the signal (defined through the Fourier transform), the number of expected peaks, and the amplitude of the signal (see Methods). We also define an extinction threshold of 10 −5 ; if a strain population falls below this it is classified as extinct. Using these distances, we performed ABC SMC for an oscillations objective  3 . 4182 models are generated forming our initial model space. A (ii) We then perform ABC SMC for an oscillatory objective, which yielded 117 models that were capable of producing oscillations. A(iii) These form the prior model space for the chaos objective, using a threshold of λ 1 > 0.003, we identify models capable of producing chaotic behaviour B The barchart shows the probability of models for the chaotic objective. The error bars represent the standard deviation. C An example time series representative of the chaos objective posterior distribution. Population densities as optical density (OD) show sustained, nonrepetitive oscillatory behaviour for the three species community.
https://doi.org/10.1371/journal.pcbi.1010548.g002 (Fig 2A(ii)). We identified 117 models capable of producing oscillations. These models become the new uniform prior model distribution for the next stage, where we perform ABC SMC for the previously described chaotic objective (Fig 2A(iii)). In this model framework we identified λ 1 > 0.003 as sufficient for classifying chaos.
The posterior probabilities of the models for the chaotic objective given the prior distributions used are shown in Fig 2B. Fig 2C shows a representative chaotic trajectory, demonstrating aperiodic non-repeating behaviour, satisfying the qualitative features of chaos.

Properties of chaotic models
We next explored some of the properties of chaotic topologies identified using ABC SMC. Fig  3A shows the top performing models when subsetting for complexity, based on the number of parts expressed. m k refers to the k-th model from the initial model space. m 850 contains four expressed parts and possesses the highest posterior probability for chaotic behaviour. Systems containing fewer parts all had a posterior probability of zero. As complexity increases to five and six parts (m 3177 and m 2547 ), the posterior probability decreases. Our previous work demonstrated that system stability increased with the number of parts [17]. The peak in the posterior probability for chaos at four parts reflects a balance that includes enough mechanisms to enable co-existence, without the tighter network of negative interactions that are associated with linear stability [17,24]. We highlight that these observations are limited to the small communities defined in our prior. These properties may not be reflective of larger communities, however, we hypothesise that a trade-off between stabilizing interactions that enable co-existence, and destabilising interactions to prevent linear stability, will remain important for producing chaotic population dynamics. Fig 3B provides summaries of how different parts contribute to chaotic behaviour in the three strain models. We can see that one QS system and positive regulation of bacteriocin is strongly favoured for producing chaos. This ensures all system bacteriocins are regulated in tandem. Expression rates are all dependent upon the same QS, resulting in stronger negative or positive correlations defined by the mode of regulation. Two bacteriocin systems also dominate the model posterior. Bacteriocin interactions can be categorised as either self-limiting (SL), whereby the strain is inhibited by the bacteriocin it produces, or other-limiting (OL) where a strain is inhibited by a bacteriocin produced by a different strain. Both SL only and a combination of SL and OL interactions are associated with producing chaotic behaviour. These observations are interesting in comparison to other work on ecological systems. Cooperative interactions were previously found to give rise to unstable systems, whereas competition was more indicative of stability [24]. The same effect might occur here in systems with one QS, rather than two, as the system would be expected to have increased correlation. While chaotic behaviour may seem to be very different from linear stability, both behaviours share the necessity for coexistence. Our previous work showed that SL interactions were important for producing linear stability, while OL interactions more frequently associated with extinction events and non-linear stability [17]. This may explain why we see tendencies for topologies to share a mixture of stability associated SL interactions, and instability associated OL interactions. We also find models with three bacteriocins, and hence higher suppression of growth, have a low posterior probability for chaos, given the prior distributions used.

Parameter importance for chaos in m 850
The model with the highest posterior probability for chaotic behaviour was m 850 ; the topology is shown in Fig 4A. It consists of a single QS system, produced by N 1 , that positively regulates two bacteriocins. B 1 is produced by N 1 and N 2 but it inhibits the growth of N 1 only. B 2 is produced by N 3 and inhibits the growth of N 3 only. The system in total consists of four expressed parts. m 850 also ranked highly for the oscillatory objective, ranking 3rd out of the initial 4182  Table 2 https://doi.org/10.1371/journal.pcbi.1010548.g004 models. This presents an interesting problem whereby a model that has promising use as an oscillator also has a high potential to produce chaos, relative to other candidate models. Identifying the parameters and initial conditions important for differentiating between chaotic and oscillatory behaviour gives us insight into how to control this behaviour when constructing genetic circuits or selecting chemostat settings.
As a first step, we analyzed the model to quantify the possible steady states and basins of attraction. Our analysis gave analytical conditions for the existence and stability for complete extinction and for single strain survival (See Methods). For three-strain co-existence, we find the following necessary conditions: This shows that for three-strain co-existence, the maximal growth rate of N 2 has to lie between certain upper and lower bounds. In particular, it has to be smaller than the maximal growth rate of N 1 or N 3 . We can see from the topology of m 850 (Fig 4A) that the growth of N 2 is not limited by any bacteriocin, therefore the only limitation on growth comes through  resource competition. If N 2 had a higher growth rate than N 1 or N 3 , it would out compete these strains and cause an extinction event.
We then wanted to explore the most important parameters that separate oscillatory and chaotic behaviours in m 850 only. We refer to a set of parameters and initial conditions as an input vector. Using ABC SMC, we performed parameter inference on m 850 for the chaotic and oscillatory objectives, generating 3750 input vectors for each objective. We can use this dataset of labelled input vectors to understand the importance of individual parameters, initial conditions and nearby steady states. Fig 4B shows multivariate parameter distributions for the oscillator and chaotic objectives for the experimentally accessible parameters. The dilution rate (D) is a directly controllable parameter of the chemostat. The production rate of A 1 (kA 1 ) can be tuned by using an inducible promoter to control expression of the AHL synthase species. Strain maximal growth rates (μ max1 , μ max2 , μ max3 ) can be controlled by using different base strains or through the combined use of auxotrophic strains and defined media. Finally, the initial population densities (N 1 , N 2 , N 3 ) can easily be set when inoculating the initial culture. Divergence between two parameter distributions indicates its importance in differentiating between the two objectives. We can see that the oscillatory objective distributions for D, N 1 and μ max2 are all constrained towards lower values relative to the prior. However, for all these distributions we can see that the chaotic and oscillatory regions overlap. This again implies that chaotic and oscillatory behaviour exist close to one another in parameter space, and highlights the multidimensional nature that determines the behaviour.
To further investigate the importance of parameters and initial conditions we trained a random forest classifier model using the input vectors as features [25]. We curated a label-balanced dataset of oscillating input vectors and chaotic input vectors. Using a train test ratio of 0.5, the trained classifier model was able to classify the test set with a *90% accuracy (Methods). Fig 4C shows the average information gain across all decision tree classifiers in the forest for all free parameters. This can be used as an indicator of feature importance in correctly classifying an input vector. K A 1 B 1 and K A 1 B 2 describe the concentration of A 1 required to produce half-maximal repression of bacteriocins B 1 and B 2 respectively. While the feature importance indicates these parameters are the most important, they are more difficult to tune compared with other parameters in this system. The error bars indicate the variability in the importance of a feature across all trees in the forest. Large error bars suggest single features are not essential for classification, and that redundancy exists between the features used [26].
From the set of chaotic input vectors, we used numerical methods to identify nearby steady states that could be reached by changing the initial state of the system only. Fig 4D shows the sensitivity analysis of a chaotic input vector. The black stars indicate stable steady states identified by numerical analysis. We perturbed the initial species values of either N 1 , B 1 or B 2 individually. The plots show how changing these initial states yields different Lyapounv exponents, highlighting the chaotic region in red. The range of Lyapunov exponents shown in Fig 4D suggest that by changing the initial conditions only we are able to produce a range of different behaviours. Perturbing N 2 , N 3 or A 1 did not produce chaotic behaviour. It is interesting that the initial state of N 1 as the A 1 producing strain appears to be more important whereas the initial concentration of A 1 itself is not.

Exploring the parameters in the transition to chaos
Being able to move a system from a chaotic state to a fixed point could be important in a bioprocess control scenario so we explored this in more detail. Previous studies have frequently identified the bioreactor dilution rate as an important parameter for transitioning between different population dynamics [27][28][29]. The posterior parameter distribution shown in Fig 4B strongly indicated the dilution rate, D, to be important for defining chaotic behaviour. We previously identified the QS production rate, k A1 and D as important parameters for transitioning between co-existence and extinction states [16]. We hypothesised that the antagonistic effect of k A1 to D would make it a useful parameter for controlling population behaviour.
First, we took an input vector known to produce chaotic behaviour and randomly sampled new values for k A1 and D from the prior and calculated the Lyapunov exponent of the new input vector. Fig 5A shows the results where filled colour indicates the maximal Lyapunov exponent calculated at each grid reference. The grid outline indicates the classification, the red grid region of Fig 5A shows the chaotic region. Fig 5A illustrates that, changing D and k A1 affects the Lyapunov exponent. The bifurcation diagrams in Fig 5B and 5C for k A1 and D respectively, illustrate the antagonistic transitions in behaviour that occur when changing the two parameters. Fig 5A and 5B show transitions through one strain extinctions (N x < 10 −5 ), stable steady state, oscillations and chaotic behaviour. Fig 5A and 5B both show that increasing k A1 results in transitions from stable co-existence, through oscillations and then to chaos, followed abruptly by an extinction event. Fig 5A and 5C and c both show that a lower dilution rate is associated with chaos; increasing the dilution rate reduces instability to produce oscillations, which abruptly transitions to a stable extinction state.
In a bioreactor control scenario it is interesting to understand if a community could be switched between states in real time. Fig 5D and 5E show how this is possible by modifying k A1 and D respectively. The red arrows on Fig 5A indicate the position of the single start point and two end points in these real-time transitions. It is important to note that when ramping up the dilution rate in real-time, we reach stable steady state in a region that would not be obtainable with a fixed dilution rate.

Discussion
We have developed a novel methodology to explore parameter regions that give rise to chaotic dynamics. We have applied it to the exploration of chaotic dynamics in synthetic microbial communities and found a high prevalence of such dynamics in these systems. This work is the first to query the existence of chaotic behaviour in synthetic microbial communities. We show that we can expect to find chaotic states in relatively small synthetic microbial systems, which has important ramifications for the field.
By first running ABC SMC for the oscillatory objective we were able to drastically reduce the model space for the search for chaos. However, the timecourse simulation and parameter sampling makes this pipeline computationally costly. In the future we can consider using eigenvalue stability methods to reject particles without simulation, improving the efficiency of our approach and therefore improving the number of samples available for posterior estimation [30,31].
We expect it will become increasingly important to consider the location of chaotic attractors in parameter space as the microbial communities we build or interact with become more complex. These methods can easily be applied to parametrise different models. It would be interesting to compare the existence of chaotic attractors in systems that use toxin-antitoxin systems [32], combinations of cooperative and competitive interactions [33], or mutualistic only interactions [34]. Genome scale metabolic models contain a large number of linear reactions [35]; they can be combined to describe microbial communities and used to model industrial bioprocesses [36,37]. Given the high dimensional nature of metabolic networks, it would be interesting to investigate whether these models yield chaotic behaviour in small community networks.

Conclusion
To conclude, we have developed methods for identifying chaotic parameter regions using ABC SMC. We have demonstrated the application of this method to resolve a previously identified chaotic attractor in a gLV model, and identified models susceptible to chaos in three strain synthetic microbial communities. Although chaotic attractors are generally thought to be sparse in low dimensional systems, we have shown their existence in realistic synthetic microbial systems. They may also exist in close proximity to stable steady state regions. This work demonstrates that deterministic chaos will be an important factor in microbial community design and should be studied in much more detail.

Three strain synthetic communities model space definition
Models are generated from a set of parts, that are expressed by different strains in the system. We represent an expression configuration through a set of options. We define the options for expression of A in each strain, where the options are: not expressed, expression of A 1 , or expression of A 2 (0, 1 or 2 respectively). We define the options for expression of bacteriocin as: no expression, expression of B 1 , expression of B 2 or expression of B 3 (0, 1, 2 and 3 respectively). Lastly we define the mode of regulation, R, for each bacteriocin, which can be either induced or repressed (0 and 1). This is redundant if a bacteriocin is not expressed.
This enables us to build possible part combinations that can be expressed by a population. Let P c be a family of sets, where each set is a unique combination of parts.
Each strain in a system can be sensitive to up to one bacteriocin. Let I represent the options for strain sensitivity. The options are: insensitive, sensitive to B 1 , sensitive to B 2 or sensitive to B 3 (0, 1, 2 and 3 respectively).
Each strain is defined by its sensitivities, and expression of parts. Let P E be all unique engineered strains: Which can be combined to form a model, yielding unique combinations: Finally, we use a series of rules to remove redundant models. A system is removed if: 1. Two or more strains are identical, concerning bacteriocin sensitivity and combination of expressed parts.
2. The QS regulating a bacteriocin is not present in the system.

3.
A strain is sensitive to a bacteriocin that does not exist in the system.

A bacteriocin exists that no strain is sensitive to.
This cleanup yields the options which are used to generate ODE equations for a system.

System equations
State variables in each system are rescaled to improve speed of obtaining numerical approximations. N X describes the concentration of a strain, B z describes the concentration of a bacteriocin and A y describes the concentration of a quorum molecule. C N , C B and C A are scaling factors: Each model is represented as sets where N defines the number of strains, B defines the set of bacteriocins and A defines the set of QS systems. The following differential equations are used to represent each model.
Growth is modelled by Monod's equation for growth limiting nutrient, S. m x max defines the maximal growth rate of the strain and K X defines the concentration of substrate required for half-maximal growth.
Killing by bacteriocin is defined by oðB 0 z Þ, where ω max defines the maximal killing rate which is set to 0 if the strain is insensitive. K ω defines the concentration at which half-maximal killing occurs.
Induction or repression of bacteriocin expression by QS, is defined by k B (z, y), where z defines the bacteriocin being expressed and y defines the quorum molecule regulating its expression. KB max z is the maximal expression rate of the bacteriocin and K B z is the concentration of quorum molecule at which bacteriocin is half-maximal. n z defines the cooperativity of the AHL binding.

Software packages and simulation settings
The ABC SMC model selection algorithm was written in python using Numpy [38], Pandas and Scipy [39]. ODE simulations were conducted in C++ with a Rosenbrock 4 stepper from the Boost library [40]. All simulations use an absolute error tolerance of 1e−9, and relative error tolerance of 1e−4. Simulations were conducted for 5000hrs, and were stopped early if the population of any strain fell below 1e−5 (extinction event). Simulations with an extinction event have distances set to maximum in order to prevent excessive time spent simulating collapsed populations.

Approximate Bayesian computation with sequential monte-carlo (ABC SMC)
Particles are first sampled from the prior distribution and simulated. A set of distances, d, are calculated from the simulation. If all distances are less than the intermediate threshold, ϵ t , the particle is accepted (d < ϵ t ). Accepted particles are weighted using importance sampling. The next population is sampled from the previous, and a new threshold is generated that is closer to the final threshold, ϵ F . This process is repeated until we reach a distance threshold of ϵ F . ABC SMC is highly parallel, allowing us to take advantage of high performance computing resources [21]. Updating ϵ t . At the end of each population of ABC SMC, the distance threshold ϵ t is updated to approach the final population, ϵ F . The quantile parameter, α, is defined. The distances of the population are sorted in ascending order and the distance at quantile α is used as the threshold for the next population. If ϵ t < ϵ F , we set ϵ t = ϵ F , marking the next population as the final population.
Algorithm 1: Algorithm for model selection with ABC SMC 1: Initialisation Set population indicator, t = 0 Set ϵ t Set final epsilon, ϵ F Set population size, N Set population particle count, i Set distance threshold quantile, α 2: Sample particle, consisting of a model (m) and parameters (θ): If If t > 0, perturb particle using perturbation kernel K t , yielding perturbed particle θ �� (m � ) 4: Simulate particle x � * f(x|θ � , m � ) 5: Calculate distance from objective d = ρ(x � , x 0 ) 6: Accept or reject particle If d > ϵ t , reject particle and go to 2 If d < ϵ t , accept particle, add θ �� , m � to population fyðm � Þ i t g 7: Set accepted particle weight Particle weight, w, is set to 1 for the initial population. For subsequent populations, the weight of a particle is equal to the probability of observing the particle given the prior, divided by the probability of observing the particle given the previous population.

Oscillatory population dynamic objective
We define the oscillatory population dynamic using three summary statistics for each strain. First, we use Fourier transform of the population signal to find the maximum frequency, f, and convert this to the period, T.
We set a minimum period of t/2 where t is the simulation time, giving us d o 1 . d o 1 . Any simulations in which T < t/2, d o 1 is set to 0, this distance ensures that all we have frequencies of oscillations that are on a scale relevant to the time period being measured. It was found that using the signal frequency alone resulted in acceptance of many simulations with very small oscillations, or simulations that rapidly dampen. We therefore generated two additional distances that account for oscillation amplitudes to select for sustained oscillations only. We can define the number of expected peaks in the simulation, p.
Peaks in the trajectory are identified by changes from a positive gradient to a negative gradient, and troughs via changes from negative gradient to positive gradient. The peak-to-peak amplitudes are calculated by differences between consecutive peaks and troughs. A K is the number of amplitudes above the threshold, K = 0.05. d o 2 is the difference between the number of expected oscillations in the simulation, and the count of above threshold oscillations. Because incomplete oscillations at the time the simulation ends can impact the distance measurement, we set a lenient final distance threshold for d o 2 :d o 3 compares the final amplitude A F in the simulation to the threshold. We set d o 3

Maximal Lyapunov exponent calculation
Lyapunov exponents can be used to measure chaotic behaviour; they describe the average exponential rate of divergence between two near trajectories of a dynamical system. The maximal Lyapunov exponent, λ 1 , can be used as determinant of chaotic behaviour. Using a method described by Sprott et al. [41], we evolve two nearby orbits and measure their average rate of separation. This directly investigates whether small changes to an initial state will produce a disproportionate separation. By periodically readjusting the distance of divergence after each time step we measure separation across a period of time, preventing a single event dominating subsequent states (Fig 6). The method is described in full by Algorithm 2. For all simulations we generate nearby orbits by perturbing one of the strain initial strain densities by 4 0 = 10 −10 . All simulations use a transient time equivalent to the first 10% of the time series.

Chaos population dynamic objective
d C 1 is the only distance for the chaotic objective. If d C 1 < 0, the particle is rejected. The final distance threshold, ϵ C , is equivalent to all λ 1 > 0.003.
For each sampled particle a prescreening process was performed to minimise time spent conducting the more computationally time consuming dual-orbit method. Simulations in which a strain fell below 10 −5 were rejected. The number of oscillations with an amplitude greater than 0.05 was counted for each strain signal. If any strain showed less than 2 oscillations the particle was rejected. ABC SMC was conducted with population sizes of 10, repeated 255 times yielding a combined final population of 2550 particles.

Random forest classifier model
Using the sci-kit learn (sklearn) python package [42], a random forest classifier was trained using 2000 estimators. The data used consisted of 3750 oscillatory input vectors, and 3750 chaotic input vectors. Training and test datasets were generated with a ratio of 0.5 by random sampling. Fig 7 shows the performance of the classifier model on the test data.
Algorithm 2: Description of dual-orbit method, demonstrated with two-dimensional system 9 S = S + log 2 (|4 1 /4 0 |) 10 Readjust y b 1 to align directionally with y a 1 , Repeat lines 6 to 11 for n iterations 13 Calculate maximal Lyapunov exponent as an average of the separation values, λ 1 = S/n Analysis of m 850 m 850 is described by the following equations By setting the left hand side of (1) to 0 we find a number of steady states P = (N 1 , N 2 , N 3 , S,  B 1 , B 2 , A 1 ).
The trivial steady state. P 0 = (0, 0, 0, S 0 , 0, 0, 0). The Jacobian of the linearisation has eigenvalues Consequently the trivial steady state always exists and is linearly stable for This shows that if the dilution rate is high enough, no strain can survive. One strain only steady states. There are three steady states where only one strain survives, P 1 , P 2 , P 3 . While P 2 , and P 3 can be calculated explicitly, P 1 is given implicitly (see below).
We start with P 2 : P 2 ¼ ð0; N 2 ; 0; S; 0; 0; 0Þ; where We see that P 2 exists provided D < m 2 max S 0 S 0 þ K : The linearisation at P 2 has eigenvalues This shows that P 2 exists and is linearly stable if D < m 2 max S 0 S 0 þ K ; and m 2 max > maxfm 1 max ; m 3 max g: The situation for P 3 is very similar: P 3 ¼ ð0; 0; N 3 ; S; 0; 0; 0Þ; where We see that P 3 exists provided The linearisation at P 3 has eigenvalues À D; This shows that P 3 exists and is linearly stable if D < m 3 max S 0 S 0 þ K ; and m 3 max > maxfm 1 max ; m 2 max g: The steady state P 1 = (N 1 , 0, 0, S, B 1 , 0, A 1 ) is more complicated and can not be expressed explicitly. Instead it is given as follows: Assume there exists a solution S to the following equation Proof: We need the solution to (2) to fulfil S < S 0 in order for A 1 to be positive. We interpret the right and left-hand-sides of (2) as a function of S, denoting them by R(S) and L(S) respectively. It is easy to see that A 1 (S) is a decreasing function of S, B 1 (S) increases as a function of A 1 and R(S) is an increasing function of B 1 . Consequently the R(S) is a decreasing function of S. We also see that R(0) = ω max > 0 and R(S 0 ) = 0. Further L(0) = −D and LðS 0 Þ ¼ m 1 max S 0 KþS 0 À D and L(S) increases as a function of S. This proves the statement. To summarise, the single-strain survival steady state requires the corresponding maximal growth rate to be large compared to other parameters.
The three-strain co-existence steady state. P 123 = (N 1 , N 2 , N 3 , S, B 1 , B 2 , A 1 ). From the equation for N 2 we obtain that S ¼ DK m 2 max À D m 2 max < minfm 1 max ; m 3 max g; Stability. Solved numerically using MATLAB. For each of the 3750 chaotic input vectors we used numerical root finding to calculate P 123 , and determined its stability by numerically determining the eigenvalues of the Jacobian. We found P 123 existed for all 3750 input vectors and was stable for 7.8% of them.