Processes 2019, 7, 163
provide researchers with a greater understanding of the impact of small perturbations on overall
system dynamics. The presented bifurcation evolution tool would facilitate such a study. With this
tool, modeling of complex biological systems which exhibit oscillatory dynamics or turning point
bifurcations can be readily optimized to achieve such behavior given the appropriate network topology.
The tool also facilitates exploration of parameter space in models for which the full range of dynamic
behavior is unknown, which could enable detection of rare dynamic behavior inside a small subspace of
the parameter landscape and inform experimental studies with possible implications for understanding
disease and the molecular underpinnings of life.
5. Conclusions
Biological systems are often described with systems of nonlinear equations to model the dynamics
of metabolic, protein, or gene-regulatory networks. While the topology of such models constrains
the range of dynamic behavior that can be achieved, parameter regimes that fully define the reaction
rates and interactions, as well as the state variables of the system, also determine whether a given
dynamic behavior will arise. Locating parameter regimes that induce bifurcations can be challenging
due to the size of the parameter landscape. In this work, a Pythonic bifurcation–evolution software
is presented which employs existing bifurcation objective functions and a differential evolution
algorithm that minimizes these objectives to improve the fitness of the best candidate parameter regime,
progressively optimizing the system for either a turning point or Hopf bifurcation. The objective
functions reward systems with steady state eigenvalues approximating those characteristic of the
desired bifurcation. The bifurcation–evolution software is validated using published models from
the BioModels Database, confirming that the algorithm performs well for evolving bistability and
oscillations. The bifurcation-evolution software was subsequently used to determine the frequency
of oscillators in randomly-generated mass-action network populations, and the results of this search
indicate that populations of large random networks with 50% reaction enrichment achieve oscillatory
dynamics more frequently than networks with fewer species or with an equivalent number of species
and reactions. This demonstrates the importance of reaction enrichment for flexibility that enables
complex dynamic behaviors. An analysis of selected randomly-generated mass-action networks from
these populations shows that negative feedback and time delays are involved in the generation of
oscillatory dynamics. Ultimately, the studies presented in the current work demonstrate the utility of
the bifurcation-evolution software for efficiently exploring parameter space for a solution that satisfies
the bifurcation objective. Due to the ubiquity of Python in computational biology, this Pythonic
bifurcation-evolution software can be readily understood and easily-integrated with existing modeling
software, providing the modeling community with a tool that could help detect of rare dynamics and
inform experimental studies of disease and the molecular mechanisms underlying a range of biological
phenomena.
Supplementary Materials: All supplementary materials required to generate the data and figures in the main text
are available at https://github.com/vporubsky/evolve-bifurcation. evolveBifurcation.py contains the sourcecode
for the bifurcation evolution algorithm described in the text. Optimized bifurcated BioModels represented in an
Antimony string format and time-course simulation output is included for all test cases. All randomly-generated
networks are provided in Antimony formats and time-course simulation output for each network is included.
Author Contributions: V.L.P. was responsible for data curation, formal analysis, investigation, methodology,
software, validation, visualization, and writing the original manuscript draft. H.M.S. was responsible for
conceptualization, funding acquisition, project administration, resources, and assisted in formal analysis and
visualization. Both V.L.P. and H.M.S. reviewed and edited the manuscript.
Funding: H.M.S. is very grateful for the support provided by National Institutes of Health grants GM123032-01,
U01HL122199-02, and P41EB023912. V.L.P. is grateful for the support of National Institutes of Health
grant P41EB023912.
Acknowledgments: V.L.P. would like to thank Kiri Choi for valuable discussions about integrating steady
state solvers and the functionality provided by Tellurium and libRoadRunner for implementing the bifurcation
evolution algorithm. V.L.P. would also like to thank J. Kyle Medley for providing insight into the utility and
shortcomings of several global and local optimization algorithms.
18
provide researchers with a greater understanding of the impact of small perturbations on overall
system dynamics. The presented bifurcation evolution tool would facilitate such a study. With this
tool, modeling of complex biological systems which exhibit oscillatory dynamics or turning point
bifurcations can be readily optimized to achieve such behavior given the appropriate network topology.
The tool also facilitates exploration of parameter space in models for which the full range of dynamic
behavior is unknown, which could enable detection of rare dynamic behavior inside a small subspace of
the parameter landscape and inform experimental studies with possible implications for understanding
disease and the molecular underpinnings of life.
5. Conclusions
Biological systems are often described with systems of nonlinear equations to model the dynamics
of metabolic, protein, or gene-regulatory networks. While the topology of such models constrains
the range of dynamic behavior that can be achieved, parameter regimes that fully define the reaction
rates and interactions, as well as the state variables of the system, also determine whether a given
dynamic behavior will arise. Locating parameter regimes that induce bifurcations can be challenging
due to the size of the parameter landscape. In this work, a Pythonic bifurcation–evolution software
is presented which employs existing bifurcation objective functions and a differential evolution
algorithm that minimizes these objectives to improve the fitness of the best candidate parameter regime,
progressively optimizing the system for either a turning point or Hopf bifurcation. The objective
functions reward systems with steady state eigenvalues approximating those characteristic of the
desired bifurcation. The bifurcation–evolution software is validated using published models from
the BioModels Database, confirming that the algorithm performs well for evolving bistability and
oscillations. The bifurcation-evolution software was subsequently used to determine the frequency
of oscillators in randomly-generated mass-action network populations, and the results of this search
indicate that populations of large random networks with 50% reaction enrichment achieve oscillatory
dynamics more frequently than networks with fewer species or with an equivalent number of species
and reactions. This demonstrates the importance of reaction enrichment for flexibility that enables
complex dynamic behaviors. An analysis of selected randomly-generated mass-action networks from
these populations shows that negative feedback and time delays are involved in the generation of
oscillatory dynamics. Ultimately, the studies presented in the current work demonstrate the utility of
the bifurcation-evolution software for efficiently exploring parameter space for a solution that satisfies
the bifurcation objective. Due to the ubiquity of Python in computational biology, this Pythonic
bifurcation-evolution software can be readily understood and easily-integrated with existing modeling
software, providing the modeling community with a tool that could help detect of rare dynamics and
inform experimental studies of disease and the molecular mechanisms underlying a range of biological
phenomena.
Supplementary Materials: All supplementary materials required to generate the data and figures in the main text
are available at https://github.com/vporubsky/evolve-bifurcation. evolveBifurcation.py contains the sourcecode
for the bifurcation evolution algorithm described in the text. Optimized bifurcated BioModels represented in an
Antimony string format and time-course simulation output is included for all test cases. All randomly-generated
networks are provided in Antimony formats and time-course simulation output for each network is included.
Author Contributions: V.L.P. was responsible for data curation, formal analysis, investigation, methodology,
software, validation, visualization, and writing the original manuscript draft. H.M.S. was responsible for
conceptualization, funding acquisition, project administration, resources, and assisted in formal analysis and
visualization. Both V.L.P. and H.M.S. reviewed and edited the manuscript.
Funding: H.M.S. is very grateful for the support provided by National Institutes of Health grants GM123032-01,
U01HL122199-02, and P41EB023912. V.L.P. is grateful for the support of National Institutes of Health
grant P41EB023912.
Acknowledgments: V.L.P. would like to thank Kiri Choi for valuable discussions about integrating steady
state solvers and the functionality provided by Tellurium and libRoadRunner for implementing the bifurcation
evolution algorithm. V.L.P. would also like to thank J. Kyle Medley for providing insight into the utility and
shortcomings of several global and local optimization algorithms.
18