processes
Article
Application of Parameter Optimization to Search for
Oscillatory Mass-Action Networks Using Python
Veronica L. Porubsky * and Herbert M. Sauro
Department of Bioengineering, University of Washington, Seattle, WA 98105, USA; hsauro@uw.edu
* Correspondence: verosky@uw.edu; Tel.: +1-206-685-2119
Received: 21 February 2019; Accepted: 13 March 2019; Published: 18 March 2019
Abstract: Biological systems can be described mathematically to model the dynamics of metabolic,
protein, or gene-regulatory networks, but locating parameter regimes that induce a particular dynamic
behavior can be challenging due to the vast parameter landscape, particularly in large models. In the
current work, a Pythonic implementation of existing bifurcation objective functions, which reward
systems that achieve a desired bifurcation behavior, is implemented to search for parameter regimes
that permit oscillations or bistability. A differential evolution algorithm progressively approximates
the specified bifurcation type while performing a global search of parameter space for a candidate with
the best fitness. The user-friendly format facilitates integration with systems biology tools, as Python
is a ubiquitous programming language. The bifurcation–evolution software is validated on published
models from the BioModels Database and used to search populations of randomly-generated
mass-action networks for oscillatory dynamics. Results of this search demonstrate the importance of
reaction enrichment to provide flexibility and enable complex dynamic behaviors, and illustrate the
role of negative feedback and time delays in generating oscillatory dynamics.
Keywords: parameter optimization; differential evolution; evolutionary algorithm; bistable switch;
oscillator; turning point bifurcation; Hopf bifurcation; biological networks; mass-action networks;
BioModels Database
1. Introduction
Biological systems exhibit dynamic behaviors due to the regulation of metabolites, proteins,
or genetic components, and these dynamics are frequently represented by a series of nonlinear
equations for the purpose of computational modeling. Dynamical behaviors in biological systems
are dependent on motifs within the network, defined by the species interactions and rate laws which
construct the overall network topology. However, the behavior of the system is also heavily influenced
by the parameter values attributed to rate constants, regulatory elements, and initial concentrations
of floating and boundary species in the network, such that the behavior may shift depending on the
current parameter regime. When modeling these biological systems, it may be desirable to obtain a
particular dynamic behavior to approximate a physiologically-relevant result. Cell cycle oscillations
have been studied for decades but underlying mechanisms remain a topic of interest to systems
biologists [1]. Neuroscientists are constructing computational models that exhibit complex oscillatory
dynamics to explore the effects of parameter variation, which enriches their understanding of disorders
like Parkinson’s and could have implications for treatment [2]. Developing such models requires
knowledge of the parameter regimes that permit complex dynamic behaviors, and this knowledge is
not always available from experimental data. Searching the landscape which defines parameter space
can be a computationally-intensive task, as this landscape is N-dimensional, where N represents the
number of parameters in the model, causing the search space to expand dramatically as the number of
parameters defining the system increases. Algorithms to scan high-dimensional parameter spaces have
Processes 2019, 7, 163; doi:10.3390/pr7030163
www.mdpi.com/journal/processes
4
Article
Application of Parameter Optimization to Search for
Oscillatory Mass-Action Networks Using Python
Veronica L. Porubsky * and Herbert M. Sauro
Department of Bioengineering, University of Washington, Seattle, WA 98105, USA; hsauro@uw.edu
* Correspondence: verosky@uw.edu; Tel.: +1-206-685-2119
Received: 21 February 2019; Accepted: 13 March 2019; Published: 18 March 2019
Abstract: Biological systems can be described mathematically to model the dynamics of metabolic,
protein, or gene-regulatory networks, but locating parameter regimes that induce a particular dynamic
behavior can be challenging due to the vast parameter landscape, particularly in large models. In the
current work, a Pythonic implementation of existing bifurcation objective functions, which reward
systems that achieve a desired bifurcation behavior, is implemented to search for parameter regimes
that permit oscillations or bistability. A differential evolution algorithm progressively approximates
the specified bifurcation type while performing a global search of parameter space for a candidate with
the best fitness. The user-friendly format facilitates integration with systems biology tools, as Python
is a ubiquitous programming language. The bifurcation–evolution software is validated on published
models from the BioModels Database and used to search populations of randomly-generated
mass-action networks for oscillatory dynamics. Results of this search demonstrate the importance of
reaction enrichment to provide flexibility and enable complex dynamic behaviors, and illustrate the
role of negative feedback and time delays in generating oscillatory dynamics.
Keywords: parameter optimization; differential evolution; evolutionary algorithm; bistable switch;
oscillator; turning point bifurcation; Hopf bifurcation; biological networks; mass-action networks;
BioModels Database
1. Introduction
Biological systems exhibit dynamic behaviors due to the regulation of metabolites, proteins,
or genetic components, and these dynamics are frequently represented by a series of nonlinear
equations for the purpose of computational modeling. Dynamical behaviors in biological systems
are dependent on motifs within the network, defined by the species interactions and rate laws which
construct the overall network topology. However, the behavior of the system is also heavily influenced
by the parameter values attributed to rate constants, regulatory elements, and initial concentrations
of floating and boundary species in the network, such that the behavior may shift depending on the
current parameter regime. When modeling these biological systems, it may be desirable to obtain a
particular dynamic behavior to approximate a physiologically-relevant result. Cell cycle oscillations
have been studied for decades but underlying mechanisms remain a topic of interest to systems
biologists [1]. Neuroscientists are constructing computational models that exhibit complex oscillatory
dynamics to explore the effects of parameter variation, which enriches their understanding of disorders
like Parkinson’s and could have implications for treatment [2]. Developing such models requires
knowledge of the parameter regimes that permit complex dynamic behaviors, and this knowledge is
not always available from experimental data. Searching the landscape which defines parameter space
can be a computationally-intensive task, as this landscape is N-dimensional, where N represents the
number of parameters in the model, causing the search space to expand dramatically as the number of
parameters defining the system increases. Algorithms to scan high-dimensional parameter spaces have
Processes 2019, 7, 163; doi:10.3390/pr7030163
www.mdpi.com/journal/processes
4