MI-Sim: A MATLAB package for the numerical analysis of microbial ecological interactions

Food-webs and other classes of ecological network motifs, are a means of describing feeding relationships between consumers and producers in an ecosystem. They have application across scales where they differ only in the underlying characteristics of the organisms and substrates describing the system. Mathematical modelling, using mechanistic approaches to describe the dynamic behaviour and properties of the system through sets of ordinary differential equations, has been used extensively in ecology. Models allow simulation of the dynamics of the various motifs and their numerical analysis provides a greater understanding of the interplay between the system components and their intrinsic properties. We have developed the MI-Sim software for use with MATLAB to allow a rigorous and rapid numerical analysis of several common ecological motifs. MI-Sim contains a series of the most commonly used motifs such as cooperation, competition and predation. It does not require detailed knowledge of mathematical analytical techniques and is offered as a single graphical user interface containing all input and output options. The tools available in the current version of MI-Sim include model simulation, steady-state existence and stability analysis, and basin of attraction analysis. The software includes seven ecological interaction motifs and seven growth function models. Unlike other system analysis tools, MI-Sim is designed as a simple and user-friendly tool specific to ecological population type models, allowing for rapid assessment of their dynamical and behavioural properties.


Introduction
Network motifs provide an approach to understand and characterise the behaviour of living systems at genomic, metabolic and ecological scales [1][2][3]. Food-webs, defined as a subset or module of larger, more complex networks, are used to analyse ecological interactions at the community or population level, as first described by mathematicians such as Lotka and Volterra, and have been widely used to explore phenomena observed at both macro-and microscales [4][5][6].
Mathematical modelling of ecological interactions is affected by the model objective (e.g., observation, prediction, control), the availability of existing knowledge and data, and the a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 structural complexity necessary to adequately describe the motif. For clarity, we define motif here to be analogous to interaction modules described by population ecologists, and the specific forms of these motifs are described widely in the literature (e.g. [7]).
The software presented here focuses on a mechanistic understanding of microbial interactions and, specifically, their analysis and simulation for two or three microbial species and associated substrates and products. The motif models are developed as systems of Ordinary Differential Equations (ODEs) used to describe the dynamics of and interactions between the individual organisms and their various components.
Mathematical analysis of such model structures is commonplace in fields such as chemostat theory [8][9][10][11], predator-prey system analysis [12,13], theoretical microbial ecology [14,15], and more recently in application to synthetic microbiology [16,17]. Methods that include steady-state analysis and basin of attraction characterisation are necessary to understand the stability, resilience and persistence of the modelled microbial populations. However, performing these analyses robustly requires a relatively high degree of competency with mathematical theory of dynamical systems. There are several tools available for the numerical analysis of dynamical ODEs (See Table 1 for details). Whilst often versatile, these tools are difficult for non-specialists to work with and are generally focused on users with some grounding in the mathematics of dynamical systems analysis. Furthermore, for use in systems with more than four ODEs, bifurcation and stability analysis is often problematic as finding exact solutions for higher-dimensional systems is non-trivial and often intractable.
We present here a mathematical analysis software, MI-Sim, specifically for use with mechanistically described microbial interactions described in terms of reactions, substrates/reactants and products. It provides a user-friendly environment in which microbiologists, microbial ecologists and biological mathematicians can rapidly and robustly characterise the dynamics of ecological motifs of up to three microbial species, without the requirement to develop their own code, programme models or rely on a detailed knowledge of the mathematics of dynamical systems analysis. MI-Sim has been developed using the proprietary software MATLAB (The Mathworks, Natick, USA).

Description of motifs
Foremost, we aimed to develop a tool that enables users to model and analyse their own species interactions by making the software as generic as possible. Here, we have taken six common ecological motifs describing interactions between two distinct species, plus one extended motif

MATCONT
A flexible MATLAB toolbox for continuation and bifurcation analysis of ODEs [20] MATDS Graphical MATLAB package for the interactive numerical study of ODE-based dynamical systems. Last modified in 2004 and still in beta version [21] PyDSTool Python-based environement for integrated simulation and analysis of dynamical models of physical systems. Although versatile, requires Python scripting to develop models and routines [22] XPPAUT Coupled solver for a range of equation types (XPP) with the widely-used bifurcation tool AUTO. It features a GUI but does require some coding of the system equations and parameters by the user that consists of three interacting species. The seven motifs, described in Table 2, are simple 'food-web' type networks that provide a theoretical basis by which scientists can test hypotheses in suitably sized community networks [7,15,24,25]. Whilst the interactions between microbial species are fixed by their motif, the substrates, products and reactions are not. The user may define these by changing the values ascribed to the parameters dictating the particular reaction kinetics of the system under investigation.

Development of the models
MI-Sim uses a deterministic rather than phenomenological approach for modelling and simulation of microbial species interactions. The described motifs are expressed as a series of ODEs, which describe the microbial growth, catabolic conversion processes, and species interactions within the system. The equations are developed using a standard mass-balance approach coupled with stoichiometric information describing the chemical transformation between reactants and products in the system. Whilst analytical approaches providing exact solutions are typically restricted to one or two species, numerical analysis allows extension to higher-dimensional models, albeit generating local rather than global solutions. The models currently available in MI-Sim take the following generalised form (shown here for one biomass and substrate pairing): where t is time, X = X(t) is the biomass concentration, S = S(t) is the substrate concentration, S in is the influent substrate concentration, μ(S) is the substrate dependent growth rate of the biomass and α is the flow of substrate in the system. In chemostat theory, for example, α typically represents the dilution rate or nutrient exchange, which is the inflow of substrate divided by the bioreactor volume. The seven motifs included in this version of the software have been developed according to this form. Users may manipulate the structure and behaviour of the different models by selecting appropriate parameter values. For example, all models contain a biomass decay rate term k dec , which may be removed by setting its value equal to zero. Bespoke models of n coupled ODEs (n being theoretically unlimited, although increasingly complex and larger systems may be restricted by CPU capacity), can be created by altering the code describing the default models, but this option is only recommended for users with some knowledge of MATLAB. Future versions of the software will include a Model Builder option to make this process more transparent.

Units
Molarity (M or mol/L) is the most commonly used unit for describing solution concentration in microbial systems. Kinetic studies of microorganism growth conventionally calculate parameters based on this unit for stoichiometric simplicity. However, in certain processes where solutes have an impact on the carbon oxidation state of the environment (e.g., wastewater treatment systems), solute strength is generally expressed in terms of Chemical Oxygen Demand (COD) for relevant components (i.e., those having a COD; M is used for those with no COD [26]). Accordingly, MI-SIM allows for the selection of units to be used by the model.

Description of the software User interface
The software is presented as a Graphical User Interface (GUI), thus removing any requirement for the user to be fluent with MATLAB coding. The user-friendly environment can be easily navigated using toolbar options for selection of motif and mathematical analysis routine, as well as visualisation and reporting options. A number of panels on the GUI enable the user to configure and parameterise the model and the simulation properties. The Growth Function panel is used to select the growth model (f n ) in the model (default: Monod). Parameterisation of the model is done using the Parameter Values panel. For dynamical analysis of the model the Initial Conditions panel provides an interface to select the starting concentrations of the model variables (substrate, S n and biomass, X n , where n is the index of the substrate/species). The Solver Options panel provides the user with a selection of in-built MATLAB ODE solver methods (default: ode23s for stiff systems) and absolute and relative tolerance values. These tolerances act to control the error estimation of the algorithm for the given solver method (e.g., Rosenbrock formula of order two for ode23s). It is not necessary for the user to understand the mechanics of the ODE solvers but, for nontrivial solutions, it may be necessary to adjust the tolerance values sacrificing speed for accuracy. A checkbox option Jacobian may reduce the stiff solver run-time by supplying a sparse Jacobian matrix to the solver. This avoids costly calls to the rate of change function. The Simulation Options panel is used to select parameters or variables and their value ranges for routines that analyse the models in a multiple-point simulation space (i.e., Multiple-point, Basin of attraction, and Phase portrait).
The remaining panels are used to display and visualise the various outputs from the mathematical analysis of the motifs. The Equations panel displays the set of differential equations for the selected motif. Running the multiple-point analysis generates a heatmap as the algorithm calculates the Jacobian matrix for each parameter pair, and after completion a table showing the stability of each steady-state will be displayed. The Progress panel displays information on percent completion, which is useful for monitoring routines requiring multiple iterations, together with some text on the task status. The Fixed Points & Stability panel displays the calculated fixed-point solutions and their stability for each variable. The output of this panel depends on the stability analysis selected in the Stability Check panel and is only available for single-point analysis. Two plot windows are available to display graphically the numerical outputs.

Analysis tools
MI-Sim comprises four numerical analysis tools used to fully characterise the motifs. They are presented here as pseudocode.
Single-point analysis. This routine is used for analysing the dynamics and steady-state stability of the motif at a single point in the model parameter space, Θ. The algorithm is summarised in Algorithm 1. Briefly, the system of ODEs, FðxðtÞÞ, are solved algebraically to find the fixed-point solutions,x g;Ã , where g represents the number of fixed-points found.x is the vector of variables (X and S) associated with the motifs. Two methods are available for analysis of the stability of the fixed-points; linear stability analysis and the Routh-Hurwitz criterion. In the former, defines the tolerance threshold for stability.

Algorithm 1: Single-point numerical analysis
Input: Θ, t,xð0Þ, q, Output:xðtÞ;x g;Ã 2 R n Â R g 1 Substitute numerical parameters Θ into system of symbolic ODEs FðxðtÞÞ 2 Solve FðxÞ numerically and plot solutionsxðtÞ 3 Assumex 2 R n ! 0 and solve algebraically for FðxðtÞÞ ¼ 0 to give fixed-point solutions,x g;Ã 4 Compute the Jacobian matrix J 5 Check stability of each fixed-point, x γ,Ã , for γ = 1, . . ., g: Method: Linear Stability Analysis Multiple-point analysis. This algorithm is an extension of the single-point routine in which the stability analysis is performed for a range of user-specified values (a i , b i ), with stepsize s i , for any given parameter pair, θ i for i = 1, 2. The resulting steady-state regions and their stability are visualised as a two-dimensional bifurcation phase plane for (θ 1 , θ 2 ). This is a powerful method for understanding the relationship between model parameters and the existence and stability of the system steady-states. In biological terms, the chosen parameters may represent operational properties of the system (e.g., dilution rate D, substrate input concentration S in ) or of the organisms themselves (e.g., growth rates, substrate yields), and thus can be used to test the effect of these parameters on the behaviour of the system itself. The method is described in Algorithm 2.

Algorithm 2: Multiple-point numerical analysis
Input: Θ 0 , fy i 2 Rjy i ! 0g, a i , b i , s i Output:xðtÞ;x g;Ã 2 R n Â R g 1 Assume ðx; y i Þ 2 R n ! 0 and solve algebraically for Fðxðt; y i ÞÞ ¼ 0 to give fixedpoint solutions,x g;Ã ðy i Þ 2 Define the Jacobian matrix J symbolically 3 for θ 1 ≔ a 1 to b 1 step s 1 do 4 for θ 2 ≔ a 2 to b 2 step s 2 do 5 Calculatex g;Ã and Jðy i ;x g;Ã Þ 6 Calculate characteristic polynomial of J : p J ðlÞ ¼ detðlI À JÞ 7 Calculate the eigenvaluesl by finding the roots of p J ðlÞ and determine stability T by: 8 if 8:l 2 Rðl < 0Þ then 9x g;Ã is stable 10 else 11x g;Ã is unstable 12 end 13 Determine steady-state, SS, from fX n 2 RjX n > 0g 14 Assign steady-state stability region as J j ¼ SS \ T 15 end 16 end 17 Plot J j ðy i Þ Basin of attraction analysis. This algorithm (Algorithm 3) is used to determine the basin of attraction for any pair of variables withinx, by assessing the influence of their initial conditions on the steady-state. Ecological systems can exhibit multiple stable states within the same parameter space. In other words, depending on the initial conditions, a system may tend towards different equilibrium points representing distinct biological behaviour. For example, the concept of resilience extends the purpose of investigating ecological systems beyond the goal of maintaining equilibrium (stability) to one in which system behaviour is reliant on the nature and character of the basins of attraction and its interaction with external factors [27]. In this way, it is informative to determine the effect that the initial conditions have in driving the system to alternate stable states. This allows for practical interpretation of operating and control strategies to ensure a system remains in the desired basin of attraction [28].
Phase portrait. The phase portrait algorithm (Algorithm 4) is an alternative approach for visualising the location of state attractors. Unlike the basin of attraction analysis, this routine shows the trajectories of the solutions for the system of ODEs. Critically, it can provide important information on the dynamics of the system in two or three dimensions.

Thermodynamics module
Mathematical modelling of microbial growth is typically tethered to empirical observations; kinetic models that are derived from experiments with parameters that may have no biological meaning (e.g., the half-saturation constant K S is merely an estimated shaping factor that is highly sensitive to both biological and physical characteristics of the process being modelled [29]). Metabolic thermodynamics have been proposed to resolve the weaknesses in purely kinetics based modelling, especially in systems with organisms growing under moderately exergonic conditions where a trade-off between kinetic and thermodynamically favourable conditions become apparent (e.g., anaerobic digestion [30]). Whilst kinetic models are restricted by their underlying mechanics, which govern the rate of change in reaction components, thermodynamics is immutable and describes the stability and direction of the same reaction. In this sense, for energy constrained systems, it defines the feasible free-energy conditions under which a species may grow. A number of models propose a framework by which single-and multiple-reactant thermodynamics may be integrated with classical growth kinetics [31][32][33].

Algorithm 4: Calculation of Phase Portrait
Input: Θ, fxð0Þ 2 R n jxð0Þ ! 0g; a n , b n , s n Output:xðtÞ 1 Solve FðxðtÞÞ for a n xð0Þ b n step s n 2 Plot FðxðtÞÞj b n a n We have implemented the simplest approach to incorporating thermodynamics, which only considers replacement of a kinetic inhibition function with one describing the effect of available energy on microbial growth. The function does not describe the thermodynamic energy required for anabolism, i.e. cellular maintenance or biomass production, and is only necessary when considering cases in which the free energy is limiting. We feel that this is an adequate method for implementing thermodynamics in models that describe an average behaviour of a microbial community, rather than individual based models where spatial diversity may also be considered [32,34]. Furthermore, more complex thermodynamic models, although useful in general microbial ecology modelling, may make the analyses provided here intractable. Extension of the thermodynamics is, therefore, not included here. The general form of the growth function using this inhibition term is given by: where k m,i is the maximum specific growth rate of species i and: where the Gibbs free energy change in the reaction is given by the standard equation: where ΔG 00 is the Gibbs free energy of all unmixed chemical species in the reaction (adjusted for temperature using Van't Hoff's equation), S π are the reaction products, S ν are the reactants/ substrates, R is the universal gas constant, and T is the temperature. This module is available upon selection of the thermodynamics option from the MI-SIM interface. The growth model functions, f n , default to Monod form in all cases.
The thermodynamic calculator comprises three steps: • System specification: The user should enter the operating temperature (K), reaction ΔG 0 (kJ/mol), and the index of reactants (r) and products (p) in the reaction. ΔG 0 can be set to zero if the user wishes to automatically calculate the ΔG value for each component of the reaction from the drop-down menus. For compounds not available in these menus, the cumulative ΔG values must be entered in the ΔG 0 box provided. Although r and p are set to one in all cases, it is necessary here to specify the actual reaction stoichiometry to properly define its thermodynamics. Hence, there may be more reactants or products involved that are not explicitly modelled by MI-SIM. A simple text display of the reaction is provided to indicate the substrates, S 0 n and biomass X n involved in the conversion. • Compound specification: In the central block of the calculator, the user may input stoichiometric information for each reactant and product featured in the reaction. This includes the number of moles of substrate and product, the molar mass of each compound (g/mol) and, if required, the option to calculate the Gibbs free energy for each compound. For the latter, the user must specify the compound using drop-down menus to choose the chemical system where it is found, and its state (e.g., solid, liquid, aqueous, gas). For the product compounds, the user must specify their stoichiometric coefficient, γ p . It should be noted that the molar mass is only required if the selected system units are mass of COD, whereby an additional conversion term (molarity to units COD, using the molar mass value) is required to calculate the ΔG rxn .
• Calculation: The calculation of the thermodynamic inhibition function is executed by pressing the Calculate button. The values for the compound specific DG 0 n at temperature T are taken from the ΔG 0 tables [35]. Polynomial curve fitting of order 2 is used to determine the ΔG 0 values for the user-specified temperature, based on the temperature dependent curves developed from the values in the cited tables.
The Add Sn to ODE function automatically defines equations for the additional substrates and products (S ν,new and S π,new ) not currently included in the selected motif. In the case additional reactants (S ν,new ) are included, a reconfiguration of the growth function (f new ) for biomass X n is required to account for its growth on these additional substrate (μ new ). The generic form of the additional equations are as follows: dX n dt ¼ À DX n þ Y new f new X n I th À k dec;n X n ð7Þ dS p;new dt ¼ À DS p;new þ g new ð1 À Y n Þf n X n I th ð8Þ The Finish button closes the thermodynamic calculator module and returns the user to the main GUI interface.

Example: A three-tiered food-web
To demonstrate the functionality of MI-Sim, we analysed a three-tiered microbial food-web with competitive and syntrophic interactions, as shown in Fig 1, and described by Eqs (10)- (19), with units expressed in COD. To parameterise the model, we took the values given in a previous analysis performed independently of MI-Sim [15], with a revised and more realistic value for the substrate affinity constant, K S,1 , shown in Table 3. The model is described by the six ODEs: where: We selected values of two operating parameters as dilution rate, D = 0.01 d −1 , such that D < max(k m,n Y m,n ), and substrate input concentration, S 1,in = 5 kgCOD/m 3 . Arbitrarily setting the initial conditions for the six variables to 0.1 and selecting the stiff ODE solver odes23s with default tolerances, we ran the single-point analysis for 100 days and observed that steady-state was reached after approximately 70 days, as shown in Fig 2a. A trajectory plot showing the dynamical relationship between species X 1 and X 3 is shown in Fig 2b. These plots can be compared to the results of the fixed-point analysis, which are displayed in the GUI window, and presented here in Table 4. The results show that four fixed points (FP),x g;Ã ! 0, are possible for the given parameters: i) complete washout of the three species (FP1), ii) washout of species X 3 (FP2/FP3), and iii) maintenance of all three species (FP4). Given the initial conditions (ICs), the steady-state reached is case iii), which is verified by examining both of the figures.
Additionally, stability analysis indicates that the system is bistable (FP1 and FP4) for the selected operating parameters, and the basin of attraction routine is subsequently used to determine which attractor (fixed-point) the system will converge to, given a range of ICs. This is shown in Fig 3a, in which the influence of changing the ICs for X 1 and X 3 is observed. The result indicates that there is a critical ratio for the two initial biomass concentrations (r 1,3 = X 1 (0)/X 3 (0)) defining a boundary between FP1 and FP4. If r 1,3 ! 0.707 and fXð0Þ 2 R 1;3 jðX 1 ð0Þ; X 3 ð0ÞÞ > 0g, then SS3 (steady-state 3 FP4) is assured. It should be noted that steady-states may comprise of one or two fixed-points in this example.
We also performed the analysis for multiple operating points (D, S 1,in ). This enables the user to visualise the transitions (bifurcations) between steady-state regions and their stability. It can be used to show the theoretical operating space necessary to attain a desired objective, such as community stability. In this example, three regions are identified corresponding to the fixed-points FP1-FP4. These regions also indicate the stability of the system at each point in the two-dimensional operating space and, together, we denote these regions as J 1 , J 2 and J 3 , as shown in Fig 3b and Table 5. It is observed that SS1 FP1 is always stable, whilst SS2 is associated with two fixed points (FP2, FP3) where it exists, one of which is always unstable.

Limitations and future work
Access. Whilst the software is currently only available for use with MATLAB, which is proprietary and may thus limit its accessibility, much of the functionality, such as access to analytical solvers, and integrated high-quality visualisation, provides an advantage over alternative open-source languages such as Octave. Development of a similar tool using C++ or Python may be plausible but limited, whilst MATLAB provides the opportunity for rapid development in a stable and multi-functional software. Furthermore, as Table 1 indicates, many of the available dynamical system analysis software are provided as MATLAB toolboxes, have versions available in MATLAB, or allow for downstream analysis with MATLAB compatible files. One drawback of the current version of is the inability of the user to define their own systems of equations. Taking MATCONT as an example [20], further development of the software will introduce a Model Builder functionality for bespoke systems of ODEs.
Computational performance. MI-Sim employs a number of algorithms that may be computationally expensive. This can be an issue in MATLAB, especially when nested loops are  Numerical analysis software for microbial interactions used. We have attempted to overcome performance concerns by using vectorisation instead of loops wherever possible and by making use of MATLAB specific commands such as matlab-Function for handling symbolic manipulation. In the case where loops are unavoidable, increasing the resolution of the output by specifying larger numbers of calculation points will significantly impact performance. This is most relevant for higher-dimensional models, such as the three-tiered motif we include in the software. This is a challenging system to solve using MATLAB, so we employ MuPad, which has a powerful symbolic engine that can solve complex systems of ODEs to generate all possible solutions. In MATLAB, the solutions may be incomplete and overriding the constraints for such systems results in failure. The MuPad option is computationally costly, but once solutions are found, subsequent analysis can be parallelised to significantly improve performance. Dimensionality. The current version of MI-Sim allows for the analysis of seven motifs comprised of two or three species, described by up to six ODEs. Theoretically, MATLAB can solve higher-dimensional systems of ODEs numerically, but this is dependent on the memory available to the CPU. For stiff problems, ODE solvers allocate memory for solution vectors of length n and a Jacobian matrix of size n × n. It is recommended to use the Jacobian solver option for these cases. The restriction to systems of six ODEs or less is necessary to ensure that the software will run on most standard CPUs, however experienced modellers may wish to test the algorithms with larger models. Future development of the software will aim to allow users to enter their equations manually via the interface and automatically configure the parameter, variable and solver options accordingly. Further, greater access to solver options and dimension reduction techniques may overcome some of the memory limitations of MI-Sim.

Software implementation
MI-Sim has been developed in MATLAB (R2015a, version 8.5.0) and tested on both R2016a and R2016b. In order to facilitate the use of the software and to motivate its further development, the MI-Sim package and documentation have been released at http://mi-sim.github.io. The software is also available for direct download on the Mathworks File Exchange at http:// uk.mathworks.com/matlabcentral/fileexchange/55492-mi-sim.

Conclusion
We have developed the MATLAB based software tool MI-Sim for numerical analysis of ecological interactions. The tool is built on a system of ODEs encompassing three to six dimensions. MI-Sim allows for the rigorous analysis of a range of standard two species interaction motifs as well as an extended motif of three species. As well as profiling of the motif dynamics for a user-defined set of parameters, MI-Sim also provides the capability to analyse the system in a two-or three-dimensional parameter space to enable investigation of system properties such as existence and stability of steady-states and bifurcation trajectories. In cases where the motif contains multi-stable regions, the basin of attraction tool is useful in understanding the role of the initial conditions on the ultimate equilibrium condition. Future releases of MI-Sim will include the possibility of including user-defined motifs for analysis, as well as a sensitivity analysis tool to explore the effect of parameter uncertainty on the behaviour and properties of the system. MI-Sim is a easy to use tool that does not require an in-depth understanding of the mathematics behind dynamical systems theory and can be used for determining the behaviour of interacting species under a desired operational parameter set. In this sense, as well as theoretical simulations, it may be useful for synthetic biologists and microbial ecologists wishing to understand the conditions or microbial properties that give rise to significant changes (performance, ecology) in simple multiple species systems.