7. Non-linear dynamic simulations

Dynamic simulations can directly be perform by importing a model that contains parameter values. For large-scale kinetic models though it is recommended to parametrize the initial-conditions to a known / assumed reference steady-state as high dimensional non-linear often exhibit multiple steady-states.

from pytfa.io.json import load_json_model
from skimpy.io.yaml import  load_yaml_model
from skimpy.analysis.oracle.load_pytfa_solution import load_concentrations, load_fluxes
from skimpy.core.parameters import ParameterValues
from skimpy.utils.namespace import *

# Units of the parameters are muM and hr
CONCENTRATION_SCALING = 1e6
TIME_SCALING = 1 # 1hr to 1min
DENSITY = 1200 # g/L
GDW_GWW_RATIO = 0.3 # Assumes 70% Water

kmodel =  load_yaml_model('./../../models/varma_strain_1.yml')
tmodel = load_json_model('./../../models/tfa_varma.json')

# Reference steady-state data
ref_solution = pd.read_csv('./../../data/tfa_reference_strains.csv',
                           index_col=0).loc['strain_1',:]

ref_concentrations = load_concentrations(ref_solution, tmodel, kmodel,
                                         concentration_scaling=CONCENTRATION_SCALING)

To run dynamic simulations the model need to contain compiled ODE expressions, by calling kmodel.prepare() and kmodel.compile_ode(). To compute fluxes at specific time points it can be useful to build a FluxFunction using the make_flux_fun(kmodel, QSSA) function.

kmodel.compile_ode(sim_type=QSSA,ncpu=8)
# make function to calculate the fluxes
flux_fun = make_flux_fun(kmodel, QSSA)

for k in kmodel.initial_conditions:
    kmodel.initial_conditions[k] = ref_concentrations[k]

desings = {'vmax_forward_LDH_D': 2.0,
           'vmax_forward_GAPD': 2.0}

fluxes = []
for p,v in desings.items():
    kmodel.parameters = parameter_values
    # Implement parameter desing
    kmodel.parameters[p].value = kmodel.parameters[p].value*v

    sol = kmodel.solve_ode(np.logspace(-9,0, 1000),
                             solver_type='cvode')

    # Calculate fluxes:
    this_fluxes = []
    for i, concentrations in sol.concentrations.iterrows():
        t_fluxes = flux_fun(concentrations, parameters=parameter_values)
        this_fluxes.append(t_fluxes)

    # Make it a DataFrame
    this_fluxes = pd.DataFrame(this_fluxes)/ref_fluxes
    this_fluxes.index = sol.time
    fluxes.append(this_fluxes)

ldh_fluxes = np.array([fluxes[0]['LDH_D'].values, fluxes[1]['LDH_D'].values ]).T

timetrace_plot(sol.time,ldh_fluxes ,
               filename='ldh_flux.html',
               legend=['2 x [LDH]','2 x [GAPD]'],
               backend='svg'
               )