Processes 2019, 7, 163
4. Discussion
4.1. Algorithm Evaluation
In this paper, a Pythonic interpretation of existing objective functions to evaluate the fitness of
system parameters for bifurcation behavior is presented, and a differential evolution algorithm is
implemented for progressive parameter optimization. The bifurcation–evolution software relies on
iterative minimization of an objective function that rewards networks with eigenvalues approximating
the specified bifurcation behavior. Several published models from the BioModels Database were
tested, and demonstrated that the algorithm performs well using default conditions for models
with few parameters, evolving the desired bifurcation behavior with a short runtime and within
a few generations. The bifurcation–evolution software could produce oscillatory dynamics in an
mitogen-activated protein kinase feedback oscillator model with 30 parameters (including 22 global
parameters and eight state variables) in less than 45 seconds. However, the largest model tested,
which contained 56 total parameters for optimization with a Hopf bifurcation objective, required
a 22 minute runtime, suggesting that the bifurcation–evolution software may not be efficient for
larger models. This extended runtime was attributed to the task of generating an initial population of
suitable members, indicating that the randomized parameter value selection created many individuals
which could not reach steady state, excluding these parameter sets from the population. It is likely
that this model is also sensitive, only reaching steady state and achieving oscillatory dynamics for a
small subset of parameter space, such that small changes in parameter values dramatically alter the
model dynamics.
The bifurcation–evolution software presents several optional arguments that the user can alter to
increase the probability of finding a suitable solution or to expedite convergence, as the stiffness of
the model will affect the efficiency of the algorithm converging on a solution. Increasing the mutation
constant and population size, or lowering the recombination constant will improve the chance of
converging on a global minima. However, these actions will also increase the time required for
convergence. For larger models which complicate the steady state solver computation, it may be
desirable to decrease the population size and increase the maximum number of iterations so that
computation time is shifted from initializing the population to performing the evolution routine.
Additionally, differential evolution is sensitive to parameter ranges. The default settings for parameter
range selection attempts to create narrow, parameter-specific ranges while still providing sufficient
space for non-trivial randomization. Providing broad ranges will generally reduce the efficiency of the
differential evolution algorithm and greatly reduce the chance of converging on a suitable solution.
However, the user can further restrict or expand these ranges when appropriate to improve the chance
of convergence.
While the results of testing the bifurcation–evolution software on multiple published biological
models and on randomly-generated mass-action networks suggest that the algorithm previously
implemented by Chickarmane et al. is suitable for detecting turning point and Hopf bifurcations,
there are some limitations. The objective function for turning point bifurcations cannot detect bistable
systems automatically. Instead, the user must pass the optimized model to external software that
permits bifurcation analysis or perform manual parameter scans to search for bistability, as in Figure 1.
The Pythonic implementation presented does not distinguish between systems with turning point
bifurcations and pitchfork or transcritical bifurcations, so it would be desirable to add an additional step
to eliminate parameter sets that correspond to these bifurcation types as in the previous implementation
in C#. Similarly, the objective function for Hopf bifurcations still permits damped oscillators to be
represented as solutions with good fitness values as determined by the eigenvalues of the optimized
parameter set, depending on the stringency of the threshold value. Chaotic oscillatory dynamics also
resulted following parameter optimization in a subset of the tested cases in the model of drosophila
circadian rhythms by Leloup et al. [21]. While chaos and birhythmicity are known to arise in this model,
the presence of these dynamics in the parameter-optimized model suggests that the algorithm does not
16
4. Discussion
4.1. Algorithm Evaluation
In this paper, a Pythonic interpretation of existing objective functions to evaluate the fitness of
system parameters for bifurcation behavior is presented, and a differential evolution algorithm is
implemented for progressive parameter optimization. The bifurcation–evolution software relies on
iterative minimization of an objective function that rewards networks with eigenvalues approximating
the specified bifurcation behavior. Several published models from the BioModels Database were
tested, and demonstrated that the algorithm performs well using default conditions for models
with few parameters, evolving the desired bifurcation behavior with a short runtime and within
a few generations. The bifurcation–evolution software could produce oscillatory dynamics in an
mitogen-activated protein kinase feedback oscillator model with 30 parameters (including 22 global
parameters and eight state variables) in less than 45 seconds. However, the largest model tested,
which contained 56 total parameters for optimization with a Hopf bifurcation objective, required
a 22 minute runtime, suggesting that the bifurcation–evolution software may not be efficient for
larger models. This extended runtime was attributed to the task of generating an initial population of
suitable members, indicating that the randomized parameter value selection created many individuals
which could not reach steady state, excluding these parameter sets from the population. It is likely
that this model is also sensitive, only reaching steady state and achieving oscillatory dynamics for a
small subset of parameter space, such that small changes in parameter values dramatically alter the
model dynamics.
The bifurcation–evolution software presents several optional arguments that the user can alter to
increase the probability of finding a suitable solution or to expedite convergence, as the stiffness of
the model will affect the efficiency of the algorithm converging on a solution. Increasing the mutation
constant and population size, or lowering the recombination constant will improve the chance of
converging on a global minima. However, these actions will also increase the time required for
convergence. For larger models which complicate the steady state solver computation, it may be
desirable to decrease the population size and increase the maximum number of iterations so that
computation time is shifted from initializing the population to performing the evolution routine.
Additionally, differential evolution is sensitive to parameter ranges. The default settings for parameter
range selection attempts to create narrow, parameter-specific ranges while still providing sufficient
space for non-trivial randomization. Providing broad ranges will generally reduce the efficiency of the
differential evolution algorithm and greatly reduce the chance of converging on a suitable solution.
However, the user can further restrict or expand these ranges when appropriate to improve the chance
of convergence.
While the results of testing the bifurcation–evolution software on multiple published biological
models and on randomly-generated mass-action networks suggest that the algorithm previously
implemented by Chickarmane et al. is suitable for detecting turning point and Hopf bifurcations,
there are some limitations. The objective function for turning point bifurcations cannot detect bistable
systems automatically. Instead, the user must pass the optimized model to external software that
permits bifurcation analysis or perform manual parameter scans to search for bistability, as in Figure 1.
The Pythonic implementation presented does not distinguish between systems with turning point
bifurcations and pitchfork or transcritical bifurcations, so it would be desirable to add an additional step
to eliminate parameter sets that correspond to these bifurcation types as in the previous implementation
in C#. Similarly, the objective function for Hopf bifurcations still permits damped oscillators to be
represented as solutions with good fitness values as determined by the eigenvalues of the optimized
parameter set, depending on the stringency of the threshold value. Chaotic oscillatory dynamics also
resulted following parameter optimization in a subset of the tested cases in the model of drosophila
circadian rhythms by Leloup et al. [21]. While chaos and birhythmicity are known to arise in this model,
the presence of these dynamics in the parameter-optimized model suggests that the algorithm does not
16