Processes 2019, 7, 163
been developed and extensively researched, using combinations of global and local searches to define
the landscape of computational models and estimate model parameters [3–5]. Still, there is a need for
efficient parameter optimization tools to search for hallmark dynamic behaviors in systems biology.
A tool implemented in C# was previously developed by Chickarmane et al. to optimize parameter
values of biological network models defined by systems of nonlinear equations for bifurcation
behavior [6]. Using information about the eigenvalues of bifurcated systems, the authors developed
objective functions to independently optimize parameters for either Hopf bifurcations, characteristic
of oscillatory systems, or for turning point bifurcations, which can lead to bistability [6]. Such
functionality would be desirable in a Pythonic computing environment for those interested in modeling
biological systems, as Python is a more ubiquitous computer language in the biological sciences,
implemented by expert and novice computer scientists alike, is easily-interpreted by the user, and
facilitates integration with existing software for modeling and simulation in systems biology. In the
current work, bifurcation–evolution software (evolveBifurcation v1.0.0, Seattle, WA, USA, 2019) is
developed in which these objective functions are adapted from C# into user-friendly Python code, and
global and local optimization algorithms are implemented for parameter evolution in computational
models available through the BioModels Database [7]. The bifurcation–evolution software is then
employed to search for oscillatory dynamics in populations of randomly-generated mass-action
kinetic models of variable size, and oscillatory models discovered during this search are analyzed to
understand how a reduced network topology generates oscillations.
2. Materials and Methods
This bifurcation–evolution software relies on standard biological network manipulation and
analysis tools available through Tellurium, a Python environment for dynamical modeling of biological
networks, and the associated library for simulation of biological models, libRoadRunner [8,9].
The algorithm implemented relies on progressively approximating an acceptable solution to
the bifurcation-specific objective function by evolving a population of parameter value vectors.
Each parameter vector represents a single point in the landscape of available parameter space that
the network can occupy, and vectors which minimize the objective function approximate the global
minimum of parameter space, where the desired bifurcation is achieved.
2.1. Objective Function
The objective functions introduced by Chickarmane et al. are re-implemented in the current work,
and enable optimization for either switch-like or oscillatory behavior, depending on the bifurcation
type selected by the user [6]. Both objective functions rely on intrinsic properties of eigenvalues
corresponding to the parameter set governing a system of nonlinear equations at steady state.
2.1.1. Optimization for Turning Point Bifurcations
Turning point bifurcations, capable of introducing bistability and switch-like behavior, can be
discovered by minimizing the following objective function as previously described [6]:
ǫ =
∏ λ i
(1 − 0.99 × e −| ∏ λ Min | )
.
(1)
A turning point bifurcation requires that one eigenvalue is zero. This objective function is effective
for evolving turning point bifurcations because the numerator, which is the product of all eigenvalues
of the system, will force the system to assume eigenvalues that approximate zero during minimization.
The denominator introduces a penalty for systems in which all eigenvalues are becoming very small,
suggesting they are all moving towards the imaginary axis [6]. λ Min includes all eigenvalues except the
smallest eigenvalue, so that no penalty results from the system achieving one zero-valued eigenvalue.
5
been developed and extensively researched, using combinations of global and local searches to define
the landscape of computational models and estimate model parameters [3–5]. Still, there is a need for
efficient parameter optimization tools to search for hallmark dynamic behaviors in systems biology.
A tool implemented in C# was previously developed by Chickarmane et al. to optimize parameter
values of biological network models defined by systems of nonlinear equations for bifurcation
behavior [6]. Using information about the eigenvalues of bifurcated systems, the authors developed
objective functions to independently optimize parameters for either Hopf bifurcations, characteristic
of oscillatory systems, or for turning point bifurcations, which can lead to bistability [6]. Such
functionality would be desirable in a Pythonic computing environment for those interested in modeling
biological systems, as Python is a more ubiquitous computer language in the biological sciences,
implemented by expert and novice computer scientists alike, is easily-interpreted by the user, and
facilitates integration with existing software for modeling and simulation in systems biology. In the
current work, bifurcation–evolution software (evolveBifurcation v1.0.0, Seattle, WA, USA, 2019) is
developed in which these objective functions are adapted from C# into user-friendly Python code, and
global and local optimization algorithms are implemented for parameter evolution in computational
models available through the BioModels Database [7]. The bifurcation–evolution software is then
employed to search for oscillatory dynamics in populations of randomly-generated mass-action
kinetic models of variable size, and oscillatory models discovered during this search are analyzed to
understand how a reduced network topology generates oscillations.
2. Materials and Methods
This bifurcation–evolution software relies on standard biological network manipulation and
analysis tools available through Tellurium, a Python environment for dynamical modeling of biological
networks, and the associated library for simulation of biological models, libRoadRunner [8,9].
The algorithm implemented relies on progressively approximating an acceptable solution to
the bifurcation-specific objective function by evolving a population of parameter value vectors.
Each parameter vector represents a single point in the landscape of available parameter space that
the network can occupy, and vectors which minimize the objective function approximate the global
minimum of parameter space, where the desired bifurcation is achieved.
2.1. Objective Function
The objective functions introduced by Chickarmane et al. are re-implemented in the current work,
and enable optimization for either switch-like or oscillatory behavior, depending on the bifurcation
type selected by the user [6]. Both objective functions rely on intrinsic properties of eigenvalues
corresponding to the parameter set governing a system of nonlinear equations at steady state.
2.1.1. Optimization for Turning Point Bifurcations
Turning point bifurcations, capable of introducing bistability and switch-like behavior, can be
discovered by minimizing the following objective function as previously described [6]:
ǫ =
∏ λ i
(1 − 0.99 × e −| ∏ λ Min | )
.
(1)
A turning point bifurcation requires that one eigenvalue is zero. This objective function is effective
for evolving turning point bifurcations because the numerator, which is the product of all eigenvalues
of the system, will force the system to assume eigenvalues that approximate zero during minimization.
The denominator introduces a penalty for systems in which all eigenvalues are becoming very small,
suggesting they are all moving towards the imaginary axis [6]. λ Min includes all eigenvalues except the
smallest eigenvalue, so that no penalty results from the system achieving one zero-valued eigenvalue.
5