Protocol for analyzing emergence dynamics of diabetes with obesity using numerical continuation and bifurcation analysis

Summary Type 2 diabetes (T2D) is a multifactorial disease that slowly and inconspicuously progresses over years. Here, we present a protocol for analyzing slow progression dynamics of T2D with obesity. We describe steps for using software to exploit the differences between the timescales of the metabolic variables and using numerical continuation and bifurcation analysis. We detail procedures to analyze bi-stable system dynamics and identify the thresholds that separate healthy and diabetic states. For complete details on the use and execution of this protocol, please refer to Yildirim et al. (2023).1


SUMMARY
Type 2 diabetes (T2D) is a multifactorial disease that slowly and inconspicuously progresses over years.Here, we present a protocol for analyzing slow progression dynamics of T2D with obesity.We describe steps for using software to exploit the differences between the timescales of the metabolic variables and using numerical continuation and bifurcation analysis.We detail procedures to analyze bi-stable system dynamics and identify the thresholds that separate healthy and diabetic states.For complete details on the use and execution of this protocol, please refer to Yildirim et al. (2023). 1

BEFORE YOU BEGIN
This section introduces necessary software packages and instructions for their installation.Here, we use the protocol to analyze the slow progression dynamics of T2D with obesity using a recently published model. 1 However, the protocol can easily be adapted for analyzing behaviors of other dynamical systems, whose dynamics are described with ordinary differential equations (ODEs).Although, some software packages are developed specifically for this protocol, they can easily be adapted for analyzing other dynamical systems.

Download and install XPPAUT
Timing: 20 min XPPAUT is a software package for simulating and analyzing mathematical models of dynamical systems. 2 XPP is another commonly used name for XPPAUT, which includes an implementation of AUTO numerical continuation software. 3Together with AUTO, XPPAUT provides a very robust and powerful platform for simulating and analyzing nonlinear dynamical systems.
1. Go to https://sites.pitt.edu/phase/bard/bardware/xpp/xpp.html to download the latest version of XPPAUT, which is freely available under GNU public license.a. Choose the suitable version for your operating system and download the folder that contains XPPAUT software.XPPAUT is available for Windows, Linux and MacOS operating systems.b.Follow the operating system specific step-by-step installation instruction.

Note:
The current version of XPPAUT runs on the X-window graphical environment.Installation guidelines for X-window servers for different operating systems are available online at the link given above.

Download and install anaconda python platform
Timing: 20 min One of the easiest ways to install Python and manage Python packages is through Anaconda platform, which provides a variety of software collections for scientific computing and machine learning applications.Anaconda also provides an advanced package and virtual environment control system (conda).
a. Download the Anaconda version that is suitable for your operating system.b.Install Anaconda following operating system specific instruction.
Note: Anaconda installation comes with Python along with a variety of packages such as NumPy, SciPy, Pandas and Matplotlib, which are necessary for this protocol.Anaconda platform can be used through both command prompt and an interactive desktop application (anaconda-navigator).
c. Locate the executable app shortcut located in your applications (programs) list, and start anaconda-navigator.

Download the model code and computer programs
Timing: 5 min The protocol is implemented for a recently published computational model that describes the dynamics of the obesity driven progression of type 2 diabetes and remission through weight loss. 1 The model describes system dynamics by a system of nonlinear ordinary differential equations (ODEs).An XPPAUT implementation of the model and other necessary programs are freely available under MIT license on GitHub repository.

Install an IDE for python
Timing: 15 min Python programs can be composed as scripts using any text editor.However, we recommend using an Integrated Development Environment (IDE).There are several freely available IDEs for Python.The choice of IDE makes no difference the way Python programs run, but it may contribute to the coding experience, and visualization of outputs, such as tables and figures.All Python installations STAR Protocols 5, 102880, March 15, 2024 Protocol come with an IDE (Idle).Here we used Spyder IDE, an open-source and freely available program, for its user-friendly interface and advanced features.
5. Download the latest version of Spyder for your operating system from https://www.spyder-ide.org/.a. Install Spyder IDE following the operating system specific instructions.
Note: If Python is installed through Anaconda platform, Spyder IDE can also be installed through Anaconda platform, which is listed in the applications on the anaconda-navigator.

MATERIALS AND EQUIPMENT
In this section, we describe the hardware requirements to run the protocol.We also provide details about the computational model used.We briefly describe XPPAUT implementation of a dynamical model and explain the theoretical basis for continuation and bifurcation analysis.

Computer hardware
The minimum hardware requirements for Python and XPPAUT are 4 GB of RAM and 5 GB of free disk space.Both Python and XPPAUT can be installed on MacOS, Windows and Linux machines.

T2D progression model
The protocol is implemented using a recently published computational model that describes the dynamics of the obesity driven progression of type 2 diabetes (T2D) and remission through weight loss. 1 Here we provide some background information and briefly describe the model.
In the model, T2D emerges as a result of insulin resistance in combination with insufficient insulin secretion, which is secondary to impaired b-cell function and reduced b-cell mass.b-cells secrete insulin in response to elevated plasma glucose levels, which induces glucose uptake in liver, muscle and adipose tissue cells.In the model, system dynamics are described by a set of coupled ordinary differential equations that determine the rates of change of body weight (body mass index, BMI), systemic inflammation level (q), insulin sensitivity (S i ), plasma free fatty acids (FFA), b-cell function The interactions between the model components are summarized in Figure 1.The details of the definitions of model variables and their mathematical formulations are described in Yildirim et al. 1 In short, the system is at a healthy state at the baseline, where the total daily energy intake and expenditure are in balance, and a healthy BMI is preserved.A perpetual imbalance between energy intake and expenditure results in weight gain that initiates a chain of events that pushes the system out of the healthy state and eventually results in insulin resistance and diabetes.In the model, obesity is linked to insulin resistance and diabetes through increased

STEP-BY-STEP METHOD DETAILS
In this section we describe a step-by-step implementation of the protocol for analyzing the way metabolic system slowly progresses into the diabetic state as a result of increased body weight, and the way weight loss interventions result in diabetes remission.We start with explaining how to execute the model on XPPAUT and generate model simulations.We then use AUTO software to perform a numerical continuation analysis and generate bifurcation diagrams of the system.We then show how to export model and bifurcation diagram data into Python and identify the thresholds that separate healthy and diabetic states.We show how to use bifurcation diagram to calculate the calorie restriction necessary for successful remission, and how to predict the longterm behavior of the metabolic system under perturbations, such as weight gain and weight loss.We also introduce software packages for generating model simulations and animations.
Executing the model on XPPAUT

Timing: 15 min
In this step, we provide details of how to execute the model on XPPAUT.We introduce the XPPAUT interface and necessary steps to run and simulate the model.For illustrations and figures, we use the model.odescript provided in KRT.
1. Open a terminal session and move to the directory where the XPPAUT is located.When the above command is run, X-window server will start automatically and XPP window in Figure 2 will pop-up.

Note:
The default background color and text size might be different on your platform.You may change the appearance of the XPPAUT either directly by declaring options within the ode file, or you can create a separate options script and call this script within the ode file.
For details see online XPPAUT documentation.Note: XPPAUT saves the time, dynamical variables and auxiliary variables data at each integration point.
c. Simulation window: When model is integrated the solution trajectories appear in this window.Currently, the window shows time in moths (TMONTH) on the x-axis and plasma glucose (G) on the y-axis.This can be changed in Viewaxes in the XPP Menu.d.XPP Menu: Main XPP menu items.Menu items can be accessed by clicking on them or by keyboard shortcuts, where the capital letter in each menu item (or the letter inside the parentheses in sub menus) is the keyboard shortcut for the item.
Note: Consecutive keyboard shortcuts can be used to access the menu>submenu items.
e.A few important XPP Menu items for the execution of the protocol are; i. Initialconds: Provides access to the menu items to initiate the integration.To initiate the integration with current parameter values and initial conditions, follow Initialconds > (G)o, or press I > G on the keyboard.The following trajectory will emerge on the simulation window that shows the plasma glucose time course for 120 moths (Figure 3).Here we set the parameter inc_i1 = 0.75 and it1 = 300 in the Param quick access menu, to represent 75% increase at daily energy intake at T = 300 days.For a list of all model parameters and their explanations, please refer to Yildirim et al. 1 Note: The (R)ange option allows initiating multiple solutions from a range of initial conditions.( 2)par range option allows generating multiple simulations for a range of parameter ii.Continue: This option allows to continue simulation until an entered new time point.After choosing Continue (or pressing C), the new time point can be entered to the XPP command line.iii.Nullcline: Allows plotting and managing nullclines of the system on the phase plane.iv.nUmerics: Allows setting up numerical solver and associated parameters.v. File: This menu allows access to the main options such as changing model.odefile elements inside the XPPAUT interface (Edit), saving model information (Save info), and access to the AUTO interface (Auto).vi.Makewindow: Menu options for creating or destroying additional simulation windows (Figure 4).

Note:
The active window is indicated by a small square on the top left corner of the window.Any window can be made the active window by clicking anywhere on the window.The simulation on the active window can be refreshed by pressing R (Restore) and the active window can be cleared by pressing E (Erase).
vii.Viewaxes: This menu allows changing the active simulation window axes in 2D or 3D.
Note: Although XPPAUT provides excellent tools for analyzing system's dynamics, its visualization capabilities are limited.Therefore, we will use Python for generating model simulations and further analysis.This can be done manually, by creating a python implementation of the model.However, here we provide a package that contains a parser function written in Python to create a python script from the mode information exported from XPPAUT.The Protocol details of this function will be discussed at step 15, but for now we will describe how to save model information data on XPPAUT.
3. On the XPP Menu, go to File > Save info.This will save the model information as the model.ode.pars file.
Note: The model information data file is provided in KRT (model.ode.pars).
Using BMI as a parameter for the continuation analysis

Timing: 30 min
In this step, we describe how to use BMI (a dynamically changing model variable) as a parameter to perform a continuation analysis.This way the effect of weight gain on the emergence dynamics of type 2 diabetes can be observed through a bifurcation diagram.This requires creating a new model, where BMI is defined as a model parameter that is used to calculate body weight.For a detailed description of bifurcation analysis and numerical continuation, please refer to the (Guckenheimer,  1980). 4te: Here we use BMI instead of body weight, since BMI provides more information about the obesity status of the subject regardless of height.
Note: In the model depicted in Figure 1, variables have fast, intermediate and slow timescales.G and I are the fast variables that act on the order of hours-days; FFA, S i , q, s, and BMI are intermediate timescale variables that act on the order of weeks-months; and b is a slow variable that acts on the order of years.The timescales of the variables determine the way they dynamically respond to perturbations, such as a change in daily energy intake (DE i ).Since we are interested in long-term slow dynamics of the system in relation to weight gain, we focus on the slow variable b and the way its steady state value changes in response to the changes in body weight or BMI.

4.
Create a new copy of model.odefile and name it model_bif.ode.a. Open model_bif.odewith a text editor, and define BMI as a parameter of the system.b.To represent W with body mass index (BMI), define following relation in the model; Where h is the height in meters, and it is a parameter of the model.This reduces the number of variables in the model by one.
5. Save model_bif.odeinto the directory, where other protocol files are located.
Note: A copy of the model_bif.odefile is included in KRT.
CRITICAL: The model presented in model_bif.odeshould only be used for generating the bifurcation diagrams of the system using BMI as the bifurcation parameter, and not for model simulations.Because the model presented in model_bif.odedoes not include body weight (W) or BMI as a dynamic variable.

Timing: 30 min
In this step, we describe how to start AUTO, and generate a bifurcation diagram of the system using model_bif.odefile and BMI as the bifurcation parameter.We will explain important numerical parameters, and propose solutions to potential problems.
6. Execute the model_bif.odefile on XPPAUT as described above.7. On the XPP menu, go to File and select Auto to start the AUTO interface (Figure 5).
Note: When AUTO is initiated, the XPP window remains open and active, where any changes made in XPP affects the AUTO.
The AUTO menu items are located on the left side of the main simulation window.Important AUTO menu items for the protocol are as follows.a. Parameter: On this menu you can specify the bifurcation parameter, which is the first parameter on the list.You should type the name of the parameter as it is declared in the ode file.AUTO defaults to the first parameter defined on the ode file (Eg0 in this case).b.Axes: On this menu, you can set the variable you want to draw on the y-axis.In other words, the major variable of interest.On this menu, you can also set the axes ranges.By default, AUTO sets the y-axis to the first dynamic variable on the ode file (G in this case).c.Numerics: On this meu you can set the numeric parameters for running continuation analysis.d.Run: This initiates the continuation calculations from the current steady-state solution.e. Grab: Selecting this option allows grabbing a point on the diagram to continue to track the solutions.
CRITICAL: If the grabbed point is a bifurcation point, this allows you to select on which branch the numerical solver to continue.For instance, when a Hopf-Bifurcation point is grabbed, you can tell AUTO to keep track of the branch of the steady state solutions or the branch of the periodic solutions.
f. Clear and reDraw options allow clearing the screen and redrawing the latest diagram, respectively.CRITICAL: When choosing a suitable Dsmax, the scale of the bifurcation parameter should be taken into account.If Dsmax is set to a large value, AUTO might miss some important points.So, if it seems to be missing a stability transition, or if the diagram is jagged, reducing this number might help.
f. Par Min[25]: This is the minimum value for the bifurcation parameter.The calculation stops if the parameter value falls below his value.g.Par Max[50]: This is the maximum value for the bifurcation parameter.The calculation stops if the parameter value increases above this value.h.Norm Max[1e+6] The upper bound for the L2 norm of the solution.If the L2 norm of the solution is greater than this number, the calculation will stop.For large models with variables on large scales, Norm Max should be set to a large value.i. Click Ok or press tab to save the changes.
Note: These internal AUTO options can be directly declared within the ode file using the code block shown below.
11. Now the numeric parameters are all set and continuation analysis can be initiated.
CRITICAL: Continuation calculations are initiated from the current initial condition, which is set on the XPP.Therefore, the current initial condition on XPP should be a steady state solution of the system.If the current initial condition of the model is not set to a steady state on XPP, AUTO may fail to converge.Please see the problem 1 and proposed solution in troubleshooting section.
a. Press the Run button on AUTO Menu, and a submenu will appear (Figure 6, top left).b.Select the 'Steady state' option on that submenu for auto to start.c.A steady state solution curve should appear on the AUTO's simulation screen (Figure 6).Note: This data will be used to continue the analysis on Python, which provides more flexibility and better visualization options.We will explain the details of the diagram in step 16.A dat file of the bifurcation diagram data is provided in KRT as bifurcation_data.dat.

Timing: 15 min
This section describes how to import model data into Python.Importing model and bifurcation diagram data into Python allows using Python's more advanced and versatile features for simulation and analysis.

Open Python IDE.
Note: On all operating systems, Spyder can be started from Anaconda-Navigator or by running Spyder command on a Terminal session.If you do not prefer using an IDE, you can use any python console to run the protocol on all operating systems.You may also use an interactive Python notebook such as Jupyter Notebook.
a. Move to the folder that contains model.ode.parsfile you saved at step 3a and bifurcation_ data.datfile saved at step 14. b.Make sure xpp2python.py and dynamics.pyfiles (see key resources table) are also in this directory.c.Use following command to import functions stored in xpp2python.pymodule into the workspace.
This module includes 3 functions.i. read_info: This function takes the name of the file that contains the model information exported from XPPAUT (model.ode.pars) as input.The function parses the information saved in the file, and returns this information as a python dictionary object.This information is used to create a python script for the model in the next step.ii.create_script: This function uses model information exported from XPPAUT to create a python script for the model.iii.simul: This function integrates the model and generates plots.d.Use following command to create a python executable script for the model.
Note: This will create a python script with the name model.py and save it to the current working directory.The script contains a function (odde) that defines the model equations, parameter values, and initial conditions as exported from the XPPAUT.This function call also saves the model information, such as variable names, parameter names, initial conditions, etc. into the info dictionary object.
Note: Model information can also directly be stored into the info dictionary object by using read_info function; e. Run following command to import the model to the workspace.

Timing: 15 min
This section describes how to plot the bifurcation diagram on Python and summarize the dynamic behavior of the system.We also show how the branch of unstable steady states serves as a threshold that separates the healthy state from the diabetic state.solutions either converge to a lower or higher b steady state depending on which side of the unstable steady state the solution was initiated.Therefore, for BMI values between 25-39 kg/ m 2 branch of unstable steady states separates the phase plane and serves as a threshold.

Timing: 15 min
The model behaves as described above for fixed BMI.However, in the model, BMI is a dynamic variable that changes over time.In the model, an increase in DE i results in weight gain, whereas a decline in DE i represents a calorie restriction or a diet intervention that results in weight loss.We now show how to use bifurcation diagram to capture the association between level of BMI and b for the emergence of T2D.We also show how to calculate the level of calorie restriction necessary for successful remission at each region on the BMI-b plane.17.Set dec_bound option to True in bif_fig function to create a bifurcation diagram with a color map that shows the necessary level of calorie restriction for remission for each BMI and b value (Figure 10).
Note: In Figure 10, the color map shows the levels of BMI, where T2D emerges due to a decline in b.The color map also shows the level of calorie restriction necessary for successful remission in each region.In the yellow region labeled as no restriction (NR) the solution trajectories are attracted by the high b stable steady state and hence T2D does not emerge.However, in the region marked by the shades of green to blue, T2D emerges due to the decline in b.Shades of green and blue also show the level of the calorie restriction necessary for remission in each region, which is calculated as the maximum DE i over the baseline necessary for successful remission in the region.
Note: Setting dec_bound = True in the bif_fig function calculates the calorie restriction boundaries for successful remission.For this, the algorithm generates a grid of initial CRITICAL: Repeating the analysis for different sets of parameters is possible, but a new bifurcation diagram also needs to be generated for the new parameter set on XPPAUT following steps 6 through 17.
18. Use animation function to create dynamic animations that show the way the phase point moves through the phase space.A GIF file of a representative animation is provided in KRT (animation.gif),where Figure 11 shows a snapshot of the animation.

> ani=animation()
Note: The trajectory is initiated from (BMI = 25, b = 1000), where DE i is increased by 75% over the baseline (BL) at T = 10 months.Initially, the phase point moves right and upward (Figure 11, green arrows).In the yellow region, BMI is increased due to the positive imbalance between energy intake and expenditure, meanwhile b is increased since the phase point is in the basin of attraction of the high b stable steady states.When the phase point passes through the threshold, the trajectory switches direction in the b-axis by entering to the basin of attraction of the lower b stable steady state.In this region the phase point keeps moving to the right, but it starts moving downward (Figure 11, magenta arrows).Inside the region marked as DE i = 150% BL, DE i must be reduced from 175% BL to a value less than or equal to the 150% of the BL for successful remission.Over time the phase point keeps moving downward towards darker regions, which implies a stricter calorie restriction for remission.At T = 70 months, the trajectory takes a sharp left turn indicating a diet intervention imposed by setting DE i to 110% of the BL (Figure 11, orange arrows).Since the phase point is inside the region marked by DE i = 120% BL, this calorie restriction is able to move the phase point to the left side of the threshold, where it starts moving upward in the b-axis (Figure 11, blue arrows).Notice that the trajectory moves left and downward after the intervention, because decline in b continues until phase point reaches to the threshold.This implies the importance of the pace of the weight loss, where weight loss on a slower pace would not be able push the phase point beyond the threshold and remission would fail.

Note:
The animation function can also be used to create custom animations for different scenarios by changing the default values of its parameters; Where, fignum is the figure number, it1 and it2 are the time points where changes in daily energy intake inc_i1 and inc_i2 are initiated, respectively.The parameters, t0 and tend determines simulation start and end times in days, respectively.Setting dec_bound = True plots the color map that shows the necessary level of calorie restriction for remission for each BMI and b value as shown in Figure 11.

Note:
The animation function returns a matplotlib.animation.FuncAnimation object (ani), which can be saved in gif file format using following command; CRITICAL: Protocol can also be run on interactive Python notebooks such as Jupyter Notebook.On these notebooks, animations may not run inline.In this case, refer to the problem 5 and suggested solution in troubleshooting section.

EXPECTED OUTCOMES
Figures 9, 10, and 11 show that for each BMI value between 25-39 kg/m 2 , branch of unstable steady states separates the phase plane and serves as a threshold.On the right side of (or below) the threshold, solutions converge to the lower b stable steady state and T2D emerges.On the other hand, on the left side of (or above) the threshold, solutions converge to a higher b stable steady state.In this section, we use the protocol to analyze the model simulations under different weight gain and weight loss scenarios.We show how to use the protocol to analyze the long-term dynamic behavior of the system.

Protocol
Running outputs function with fn = 12 option crates a figure with 3 solution trajectories superimposed onto the bifurcation diagram (Figure 12).
Figure 12 shows the bifurcation diagram along with superimposed solution trajectories.Green trajectory remains on the left-hand side of the threshold and converges to a high b stable steady state (Figure 12D).Therefore, diabetes does not emerge at all, where glucose levels always remain below diabetes cut off (Figure 12C, green).However, the orange trajectory passes through the threshold and converges to a low b stable steady state (Figure 12D, orange), which results in high glucose and diabetes (Figure 12C, orange).On the other hand, magenta trajectory passes to the right hand-side of the threshold, where it starts moving downward towards a lower b stable steady state.In magenta simulations, diabetes also Figure 12D shows that trajectories need to pass through the threshold to switch direction in the b-axis.Once the trajectory passes to the right side of the threshold, diabetes emerges on the time-scale of b.Therefore, remission can only take place if a weight loss intervention pushes the phase point back to the left side of the threshold.This depends on the time spent on the right side of the threshold and the intensity of the calorie restriction.We will now show how the level of calorie restriction is related to the success remission.
Running outputs function with fn = 13 creates Figure 13.
Figure 13 shows simulations for 3 different scenarios.For green simulation, DE i is increased by 175% over the baseline (BL) at T = 10 months and set back to 120% of the BL at T = 70 months.For orange simulation, DE i is increased by 175% over the baseline (BL) at T = 10 months and set back to 130% of the BL at T = 70 months.For magenta simulation, DE i is increased by 175% over the baseline (BL) at T = 10 months and set back to the BL at T = 88 months.In green simulation, the remission is achieved because the DE i is reduced to 120% of the baseline, while the phase point is inside the region marked by the 120% BL.On the other hand, at this same point, reducing DE i to 130% of the baseline fails in remission, since this intervention cannot push the phase point to the left side of the threshold (Figure 13D, orange).Finally setting DE i to 100% BL after 88 months fails in remission as well (Figure 13 it1: The time point to initiate the increase in daily energy intake over the baseline (days, the default is 300 days).
inc_i1s: The list of increase in daily energy intake over the baseline (the default is [0.7, 0.75]).
it2s: The list of intervention time points (the default is [2100,2400] days).
inc_i2s: The list daily energy intake at the intervention time-point (it2) calculated over the baseline (the default is [0.5, 0.2]).An inc_i2 value above inc_i1 indicates a further increase in daily energy intake at it2, and weight gain.
dec_b: the option whether to draw the color map (the default is False).
> outputs(fn=13) labels: The labels for each simulation (the default is a list of simulation number, where a list of strings should be provided for custom labels).
Note: Running the main method in main.pyscript in KRT generates all the Python figures and animations.

LIMITATIONS
This protocol describes how to use bifurcation analysis to understand the slow dynamics of the emergence and progression of diabetes with obesity using a recent model. 1 However, the protocol and provided software can be used to analyze any models that use ordinary differential equations to describe the system's dynamics.Although the software packages were designed and described with general use in mind, the software might need to be adjusted for different use cases.The other limitation of the protocol is that it gives best results if system variables operate on different timescales.This is particularly true for high dimensional models, since a two-dimensional projection of the bifurcation diagram is studied for the analysis.A larger separation between the time scales of the system variables makes analysis and results more robust.

TROUBLESHOOTING Problem 1
The AUTO needs to be from a steady state to keep track of the steady state solutions as the bifurcation parameter changed in either direction.AUTO uses the initial conditions set in XPP as the starting point.Therefore, initial conditions in XPP must be set to a steady state.If AUTO is initiated from a point that is far from a steady state, the convergence fails and a single point appears in AUTO simulation screen (Figure 14) (related to Step 11).

Potential solution
Go to the XPP window.Switch the y-axis on the main simulation screen to B using Viewaxes menu, and set the axes limits.Press P on the keyboard to allow parameter input.Type in bmi and press enter.Now type in the desired value for bmi.

Problem 2
If AUTO exceeds the maximum number of steps (Nmax) and stops prematurely before completing a solution branch (Figure 15A) you may continue from where it left by using the Grab option (related to Step 11).To better describe this problem and possible solution, we set Nmax to 800, a relatively small value for following illustration.

Potential solution
Click the Grab button on the menu, and a blue cross will appear on the initial point (Figure 15A).This allows you to navigate through the branch.You can use the arrow keys to move on step at a time or you can press tab to jump through Npr number of steps set in Numerics.
Move the blue cross to the end of the branch.Notice that information down the screen will change as you move the blue cross through the curve.When you reach to the end of the curve, Ty will say EP (Figure 15B), which stands for end point indicating the last point of a run.Press enter to grab the EP and then click Run for AUTO to continue.AUTO takes another Nmax steps and stops again (Figure 15C).Notice system goes through a saddle node bifurcation (SNB), where curve turns down and leftward.
(s), b-cell mass (b), plasma glucose (G) and plasma insulin (I).The timescales of the model variables range from days to years and their interactions give rise to the emergence of insulin resistance and diabetes.
FFA and inflammation levels, which reduce insulin sensitivity and b-cell function.Reduced b-cell function causes a decline in insulin levels.Together with reduced insulin sensitivity, reduced insulin levels lead to severe hyperglycemia.Prolonged hyperglycemia pushes b-cell metabolic rate to its limits and results in a further decline in b-cell function, and eventually results in b-cell death.In the model T2D emerges as a result of reduced insulin sensitivity and insufficient insulin secretion that results from reduced b-cell function and b-cell mass.

Figure 1 .
Figure 1.Model diagram DE i : Daily energy intake, W: Body weight, q: Systemic inflammation, FFA: Plasma free fatty acids, S i : Insulin sensitivity, G: Plasma glucose, I: Plasma insulin, s: b-cell function, b: b-cell mass.Green arrows indicate a positive/stimulatory relation, whereas red arrows indicate an inhibitory/negative relation.The timescales of the variables are color coded; green: slow variables, blue: intermediate timescale variables, orange: fast variables, white: parameters/constants.

>
Figure 2. XPPAUT interface The screenshot of XPPAUT interface after model execution.

2 .
Major items on the XPPAUT's main window are; a. Title: The version of XPPAUT and the name of the executed ode file.b.Quick access menu: A list of items that can be accessed quickly, which include initial conditions (ICs), boundary conditions (BCs), delays, parameter values (Param), equations (Eqns), and Data.

Figure 4 .
Figure 4.Additional XPP simulation windows can be created using Makewindow menu The currently active window is indicated by a small square on the top left corner of the window.

Figure 5 .
Figure 5.The AUTO interface Figure 6.First segment of the bifurcation diagram on AUTO Upper stable (red) and unstable (white/black) branches of the bifurcation diagram.

Figure 7 .
Figure 7.The bifurcation diagram Stable steady state solutions (red) and unstable steady state solutions (white/black) of the system.

>
Figure 8. Model simulations generated using the simul function For blue simulations default parameters were used.For orange simulations inc_i1 is increased to 0.75 and it1 is set to 300 to simulate the effect of increasing daily energy intake by 75% at T = 300 days (10 months).The second set of simulations are superimposed on the existing figure by setting keep_old = True and fign to the current figure number.
Figure 9.The bifurcation diagram of the system projected onto BMI-b axes Red curves show the stable steady states of b, whereas black curve shows the unstable steady states of b for corresponding BMI value.Vertical dotted lines show BMI = 35, 39 and 45 cross sections, where red dots show stable steady states, and black dot shows unstable steady state.The black arrows show the way solutions behave near each steady state.The vertical lines, arrows and text do not appear on the figure automatically, and they are added later for illustration purposes.

Figure 10 .
Figure 10.Bifurcation diagram with calorie restriction boundaries for successful remission The color bar shows the maximum DE i levels relative to baseline (percentage of BL) necessary for successful remission within the corresponding region.NR: no restriction; F: failure.

Figure 11 .
Figure 11.A snapshot of the animation created using the animation() function BMI (A), b (B) and Glucose (C) time courses are shown.(D) Solution trajectory (blue) is superimposed onto the bifurcation diagram.The DE i is increased by 175% over the baseline (BL) at T = 10 months and set back to 110% of the BL at T = 70 months.In animations, the filled circle represents the phase point and it navigates through the trajectories.The color of circle is mapped to the glucose level shown on the color bar as an indicator of diabetes status.The arrows show the pre-diabetes (green), diabetes emergence (magenta), weight loss (orange) and remission (blue) phases.In panel D, the size of the circle is proportional to BMI.

Figure 12 .
Figure 12.Solution trajectories against the bifurcation diagram BMI (A), b (B) and Glucose (C) time courses are shown for different parameter sets.(D) Solution trajectories are superimposed onto the bifurcation diagram.For green simulation, DE i is increased by 50% over the baseline (BL) at T = 10 months.For orange simulation, DE i is increased by 75% at T = 10 months.For magenta simulation, DE i is increased by 75% at T = 10 months and set back to the baseline at T = 70 months to simulate the effect of a calorie restriction.Purple and green arrows in panels B and D show the peak and nadir of green trajectory from a different perspective, respectively.In panel C, black and red horizontal dotted lines mark the glucose cut offs for emergence of prediabetes (100 mg/dL) and T2D (125 mg/dL).
Figure13shows simulations for 3 different scenarios.For green simulation, DE i is increased by 175% over the baseline (BL) at T = 10 months and set back to 120% of the BL at T = 70 months.For orange simulation, DE i is increased by 175% over the baseline (BL) at T = 10 months and set back to 130% of the BL at T = 70 months.For magenta simulation, DE i is increased by 175% over the baseline (BL) at T = 10 months and set back to the BL at T = 88 months.In green simulation, the remission is achieved because the DE i is reduced to 120% of the baseline, while the phase point is inside the region marked by the 120% BL.On the other hand, at this same point, reducing DE i to 130% of the baseline fails in remission, since this intervention cannot push the phase point to the left side of the threshold (Figure13D, orange).Finally setting DE i to 100% BL after 88 months fails in remission as well (Figure13, magenta).Because at this point, the phase point has already entered the dark blue region marked by failure (Figure 13D, magenta).

Figure 13 .
Figure 13.Bifurcation diagram with calorie restriction boundaries BMI (A), b (B) and Glucose (C) time courses are shown for different parameter sets.(D) Solution trajectories are superimposed onto the bifurcation diagram along with decision boundaries that show the necessary calorie restriction for remission.For green simulation, DE i is increased by 175% over the baseline (BL) at T = 10 months and set back to 120% of the BL at T = 70 months.For orange simulation, DE i is increased by 175% over the baseline (BL) at T = 10 months and set back to 130% of the BL at T = 70 months.For magenta simulation, DE i is increased by 175% over the baseline (BL) at T = 10 months and set back to the BL at T = 88 months.Purple and green arrows in panels B and D show the peak and nadir of green trajectory from different perspectives, respectively.In panel C, black and red horizontal dotted lines mark the prediabetes (100 mg/dL) and diabetes (125 mg/dL) glucose cut offs.The color bar shows the maximum DE i levels relative to baseline (percentage of BL) necessary for successful remission within the corresponding region.BL: baseline; NR: no restriction; F: failure.

Note:
You may also click on Param on the Quick Access Menu on XPP, which lists all parameters and their values.Set bmi to desired value and click Go to initiate the simulation.Run the simulation by choosing Initialcond > (G)o (or press I > G on keyboard).Then repeat the simulation using the last point from previous simulation as new initial condition by choosing Initialcond > (L)ast (or press I > L on keyboard).Repeat the last step a few times until simulation results in a horizontal line, indicating variable does not change and system is at equilibrium.Go back to AUTO, and go to File > Reset Diagram to restart the process.

Figure 14 .
Figure 14.AUTO fails to converge if it is not started from a steady state

Figure 15 .
Figure 15.AUTO stops if it exceeds maximum step size Nmax (A) Pressing Grab a blue cross appears on the initial point.(B) Moving blue cross to the end of the curve (EP) to grab that point.(C) Restarting AUTO from the end of the first run AUTO continues for another Nmax.Curve takes a turn at a saddle node bifurcation (SNB).(D) Restarting AUTO from the end of the second run with a negative Ds to have bifurcation parameter to change in negative direction.AUTO completes the branch after 3 grabs and runs.