Processes 2019, 7, 163
2.6. Machine Specifications
All computations with runtime calculations were performed on an Intel(R) Core(TM) i5-6300HQ
CPU (Intel Corporation, Santa Clara, CA, USA) at 2.30 GHz with 8.00 GB RAM.
2.7. Data Repository
All data required to construct the figures in the main text and the bifurcation evolution algorithm
sourcecode, evolveBifurcation.py, are publically available on Github. The repository location is
provided in the Supplementary Materials.
3. Results
3.1. Testing Bifurcation–Evolution Software on Models from the BioModels Database
To demonstrate the efficacy of the bifurcation–evolution software, models from the BioModels
Database underwent parameter optimization for turning point or Hopf bifurcations, depending on the
dynamic properties described for each model in the referenced publications. The rate laws describing
the models tested are not limited to mass-action kinetics, and demonstrate that the algorithm is
effective for optimization of more complex systems. The results of these test cases are shown in Table 1.
Graphical output of the optimized networks and model files are available in the supplementary data.
Most of these models were tested in the previous work, so we have demonstrated that the Pythonic
implementation maintains functionality for all previous test cases [6]. All models were evolved using
the default settings in the bifurcation–evolution software, with the exception of the threshold value,
which was set independently for turning point and Hopf bifurcations. All turning point models were
evolved with a fitness threshold of 10 −3 to induce bistability. All models capable of a Hopf bifurcation
were evolved with a fitness threshold of 5 to optimize for oscillatory dynamics. These threshold
values were chosen empirically. The runtime is the average number of seconds to complete a single
optimization, taken over 100 attempts.
The largest model tested, a negative feedback and bi-rhythmic oscillator containing 10 state
variables and 46 global parameters, could achieve oscillatory dynamics after optimization [21].
However, a single run with the default parameters in the bifurcation–evolution software lasted 22
minutes, with the majority of this time allocated to initializing a population of 50 members that could
achieve steady state before the differential evolution routine could begin. In some of the simulations
of the optimized model, chaotic oscillatory behavior described by the authors arose, which was not
observed in any of the other models tested.
Figure 2 shows the result of optimizing for a turning point bifurcation in the Hervagault bistable
switch model [22]. Before parameter optimization, the concentrations of S1 and S2 reach a single
steady state despite a parameter sweep in the global parameter, J1_k. After parameter optimization,
both S1 and S2 achieve two distinct steady states, demonstrating the bistable dynamics of the model
in a parameter regime with eigenvalues that satisfy the condition for a turning point bifurcation.
Figure 3 shows the result of optimizing for a Hopf bifurcation in the modified Edelstein relaxation
oscillator model [23]. Before parameter optimization, species A reaches a steady state upon simulation.
However, once optimized, species A achieves sustained oscillatory dynamics.
3.2. Oscillation Discovery in Randomly-Generated Networks
The bifurcation–evolution software was used to search populations of randomly-generated
networks for models exhibiting oscillatory dynamics. Table 2 shows the percentage of
sustained oscillators in populations of randomly-generated networks with variable network sizes.
Networks either have an equal ratio (1:1) of species to reactions, or are enriched with 50% more
reactions than species (1:1.5). Figure 4 shows the frequencies of oscillatory dynamics in networks of
variable size. The top row of Figure 4 contains frequency data of all oscillating systems, including
those with sustained and damped dynamics. The bottom row contains only the frequency of sustained
10
2.6. Machine Specifications
All computations with runtime calculations were performed on an Intel(R) Core(TM) i5-6300HQ
CPU (Intel Corporation, Santa Clara, CA, USA) at 2.30 GHz with 8.00 GB RAM.
2.7. Data Repository
All data required to construct the figures in the main text and the bifurcation evolution algorithm
sourcecode, evolveBifurcation.py, are publically available on Github. The repository location is
provided in the Supplementary Materials.
3. Results
3.1. Testing Bifurcation–Evolution Software on Models from the BioModels Database
To demonstrate the efficacy of the bifurcation–evolution software, models from the BioModels
Database underwent parameter optimization for turning point or Hopf bifurcations, depending on the
dynamic properties described for each model in the referenced publications. The rate laws describing
the models tested are not limited to mass-action kinetics, and demonstrate that the algorithm is
effective for optimization of more complex systems. The results of these test cases are shown in Table 1.
Graphical output of the optimized networks and model files are available in the supplementary data.
Most of these models were tested in the previous work, so we have demonstrated that the Pythonic
implementation maintains functionality for all previous test cases [6]. All models were evolved using
the default settings in the bifurcation–evolution software, with the exception of the threshold value,
which was set independently for turning point and Hopf bifurcations. All turning point models were
evolved with a fitness threshold of 10 −3 to induce bistability. All models capable of a Hopf bifurcation
were evolved with a fitness threshold of 5 to optimize for oscillatory dynamics. These threshold
values were chosen empirically. The runtime is the average number of seconds to complete a single
optimization, taken over 100 attempts.
The largest model tested, a negative feedback and bi-rhythmic oscillator containing 10 state
variables and 46 global parameters, could achieve oscillatory dynamics after optimization [21].
However, a single run with the default parameters in the bifurcation–evolution software lasted 22
minutes, with the majority of this time allocated to initializing a population of 50 members that could
achieve steady state before the differential evolution routine could begin. In some of the simulations
of the optimized model, chaotic oscillatory behavior described by the authors arose, which was not
observed in any of the other models tested.
Figure 2 shows the result of optimizing for a turning point bifurcation in the Hervagault bistable
switch model [22]. Before parameter optimization, the concentrations of S1 and S2 reach a single
steady state despite a parameter sweep in the global parameter, J1_k. After parameter optimization,
both S1 and S2 achieve two distinct steady states, demonstrating the bistable dynamics of the model
in a parameter regime with eigenvalues that satisfy the condition for a turning point bifurcation.
Figure 3 shows the result of optimizing for a Hopf bifurcation in the modified Edelstein relaxation
oscillator model [23]. Before parameter optimization, species A reaches a steady state upon simulation.
However, once optimized, species A achieves sustained oscillatory dynamics.
3.2. Oscillation Discovery in Randomly-Generated Networks
The bifurcation–evolution software was used to search populations of randomly-generated
networks for models exhibiting oscillatory dynamics. Table 2 shows the percentage of
sustained oscillators in populations of randomly-generated networks with variable network sizes.
Networks either have an equal ratio (1:1) of species to reactions, or are enriched with 50% more
reactions than species (1:1.5). Figure 4 shows the frequencies of oscillatory dynamics in networks of
variable size. The top row of Figure 4 contains frequency data of all oscillating systems, including
those with sustained and damped dynamics. The bottom row contains only the frequency of sustained
10