Acceleration Strategies to Enhance Metabolic Ensemble Modeling Performance

Abstract

Developing reliable, predictive kinetic models of metabolism is a difficult, yet necessary, priority toward understanding and deliberately altering cellular behavior. Constraint-based modeling has enabled the fields of metabolic engineering and systems biology to make great strides in interrogating cellular metabolism but does not provide sufficient insight into regulation or kinetic limitations of metabolic pathways. Moreover, the growth-optimized assumptions that constraint-based models often rely on do not hold when studying stationary or persistor cell populations. However, developing kinetic models provides many unique challenges, as many of the kinetic parameters and rate laws governing individual enzymes are unknown. Ensemble modeling (EM) was developed to circumnavigate this challenge and effectively sample the large kinetic parameter solution space using consistent experimental datasets. Unfortunately, EM, in its base form, requires long solve times to complete and often leads to unstable kinetic model predictions. Furthermore, these limitations scale prohibitively with increasing model size. As larger metabolic models are developed with increasing genetic information and experimental validation, the demand to incorporate kinetic information increases. Therefore, in this work, we have begun to tackle the challenges of EM by introducing additional steps to the existing method framework specifically through reducing computation time and optimizing parameter sampling. We first reduce the structural complexity of the network by removing dependent species, and second, we sample locally stable parameter sets to reflect realistic biological states of cells. Lastly, we presort the screening data to eliminate the most incorrect predictions in the earliest screening stages, saving further calculations in later stages. Our complementary improvements to this EM framework are easily incorporated into concurrent EM efforts and broaden the application opportunities and accessibility of kinetic modeling across the field.

Figure 1

Modifications to base ensemble modeling framework. Given here is a general overview of base ensemble modeling framework highlighting where to apply acceleration strategies. In this work, we have implemented conservation analysis to reduce our metabolic network and only track independent metabolite and enzyme fraction species (Acceleration Strategy 1). We have also characterized the local stability check after implementing it after the initial kinetic parameter sampling process (Acceleration Strategy 2). Lastly, we have described the benefits of preranking perturbation data sets before screening the ensemble (Acceleration Strategy 3). To see this figure in color, go online.

Figure 2

Test metabolic networks used throughout study. (A) Given here is a small model with six enzymatic reactions adapted from a preexisting toy model (17). (B) Given here is a medium model with 34 enzymatic reactions and three regulation reactions created to resemble a simplified central carbon metabolism. (C) Shown here is a large model with 138 enzymatic reactions and 60 regulation reactions adapted from the core E. coli model provided by Khodayari et al. (23) The line color denotes the reaction type: blue (internal reaction), green (exchange reaction), and red (reaction knocked out during screening and time trial testing). Metabolites are represented by blue circles, and common cofactors present in the medium and large models are uniquely colored to distinguish their repetitive occurrence throughout the models. To see this figure in color, go online.

Figure 3

Conservation analysis and local stability check improve solve time. Given here is the total CPU time for integration of 700 ODE calculations (100 kinetic parameter sets for large model perturbed for seven unique enzyme knockouts) using different ensemble modeling frameworks: base framework, reduced network after conversation analysis, and reduced network after conservation analysis with local stability check. Solve times compared using Student’s t-test indicated p < 0.00001 between all condition pairs. Error bars are SD (n = 3).

Figure 4

Solve-time improvements conferred by conservation analysis and local stability check scale with model size. Given here is the percent improvement in solve time from the base EM framework for the small, medium, and large models after implementing conservation analysis and the local stability check. Small model is 100 kinetic parameter sets screened against one knockout. Medium model is 100 kinetic parameter sets screened against seven knockouts. Large model is 100 kinetic parameter sets screened against seven knockouts. Student’s t-test was performed on average solve times across three trials. Average solve times between implementing conservation analysis treatment and implementing the conservation analysis and stability test for the medium model (p < 0.05) and the large model (p < 0.00001) were significant. The small model did not see significant improvement when adding the stability test. Error bars are SDs using propagation of error for percentage change calculations (n = 3). To see this figure in color, go online.

Figure 5

Screening KOs in rank order eliminates unpredictive models earlier. (A) The solution fitness averaged for large models across predictions for all seven knockouts of the initial ensemble of 10,000 locally stable models is 0.042. Moving across the figure from left to right, the initial ensemble is screened against knockouts one at a time, ranked from most to least different (blue) and least to most different (green) flux distribution determined by clustering. The average solution fitness for all seven knockouts is plotted after each additional knockout screen. (B) Total additive solve time after each screening step is recorded in total CPU hours. When the knockouts are fed from least to most different from the WT for the large model, the total time to screen all seven knockouts is 118 total CPU hours. The time to screen all seven knockouts when fed from most to least different is 225 total CPU hours. The simulations were run in parallel across 12 nodes. Error bars are SD (n = 3). (C) Regardless of knockout screening order, the initial ensemble of 10,000 kinetic parameter sets for the large model is reduced to 26 kinetic parameter sets. However, when knockouts are fed from highest screening power to lowest (blue), the ensemble size is reduced by >98% to 152 kinetic parameter sets after the first screening step. Alternatively, when the knockouts are fed from lowest to highest screening power (green), the ensemble size is only reduced 4% to 9580 parameter sets. Inset plot zooms in on the bottom 2% of remaining models. To see this figure in color, go online.