diff --git a/doc/ORACLE_parameter_sampling.rst b/doc/ORACLE_parameter_sampling.rst index f329178..ecf9609 100644 --- a/doc/ORACLE_parameter_sampling.rst +++ b/doc/ORACLE_parameter_sampling.rst @@ -1,10 +1,9 @@ ORACLE parameter sampling =========================== -We here present the first open-source implementation of the ORACLE parametrization method [3-6]. ORACLE makes use of sampling alternative steady states and parameter backcalculation to effeciently sample large sets of locally stable parameters. +We here present the first open-source implementation of the ORACLE parametrization method [3-6]. ORACLE makes use of sampling alternative steady states and parameter backcalculation to efficiently sample large sets of locally stable parameters. -The first step in the ORACLE workflow is to integrate physological data. Here we integrare observed concentrations and we assume that the glucose transporter operates far from equilibrium. - +The first step in the ORACLE workflow is to integrate physiological data. Here we integrare observed concentrations and we assume that the glucose transporter operates far from equilibrium. .. code-block:: python @@ -18,7 +17,6 @@ The first step in the ORACLE workflow is to integrate physological data. Here we tmodel.solver = GLPK tmodel.solver.configuration.tolerances.feasibility = 1e-9 - #tmodel.solver.configuration.tolerances.optimality = 1e-9 # Define the reactor medium #Glucose to 10 g/L = 0.056 mol/l 1/2 Reactor concentration @@ -42,19 +40,6 @@ The first step in the ORACLE workflow is to integrate physological data. Here we # Enforce glucose transporter displacements tmodel.thermo_displacement.GLCptspp.variable.ub = -2.0 - # Reduce the other - # Constraint non-medium concentrations to be lower then muM - LC_SECRETION = np.log(1e-6) - secretions = [r for r in tmodel.boundary if r.upper_bound <= 0] - for sec in secretions: - for met in sec.metabolites: - if met.id in ['h2o_2', 'h_e']: - continue - try: - tmodel.log_concentration.get_by_id(met.id).variable.upper_bound = LC_SECRETION - except KeyError: - pass - # Test feasiblity print(tmodel.optimize()) @@ -70,8 +55,24 @@ In the next step we sample the flux and concentration space: samples = sample(tmodel, NUM_TFA_SAMPLES, method='achr') samples.to_csv('./samples_fdp1_1000.csv'.format()) - -Finally we call the parameter sampler for each sampled reference steady-state. The parameter sampler provides then information on the charateristic timescales this information can then later be used to discard parameter samples exhibting unphysiological charateristics. + +Fluxes sampled in the TFA problem are mass fluxes in units of mmol/gDW/hr we convert these +to reaction rates in units of [concentration]/[time] therefore we require the cell density and +the wet weight to dry weight ratio. It is also useful to choose a different concentration unit than [M] since +most intracellular metabolite concentration are between 1 mM and 1 uM, which can be achived by choosing a different +concentration scaling factor: ``[C] = [M] x CONCENTRATION_SCALING``. + +.. code-block:: python + + N_PARAMETER_SAMPLES = 10 + CONCENTRATION_SCALING = 1e6 # 1 umol + TIME_SCALING = 1 # 1hr + DENSITY = 1200 # g/L + GDW_GWW_RATIO = 0.3 # Assumes 70% Water + + +Finally we call the parameter sampler for each sampled reference steady-state. The parameter sampler provides then information on the charateristic timescales this information can then later be used to discard parameter samples exhibiting unphysiological characteristics. +The parameter sampler will requires that the Jacobian is compiled ``kmodel.compile_jacobian()`` to evaluate the stablity of the samples. .. code-block:: python @@ -84,17 +85,6 @@ Finally we call the parameter sampler for each sampled reference steady-state. T NCPU = 8 - N_PARAMETER_SAMPLES = 10 - CONCENTRATION_SCALING = 1e6 # 1 mol to 1 mmol - TIME_SCALING = 1 # 1hr - # Parameters of the E. Coli cell - DENSITY = 1200 # g/L - GDW_GWW_RATIO = 0.3 # Assumes 70% Water - - flux_scaling_factor = 1e-3 / (GDW_GWW_RATIO / DENSITY) \ - * CONCENTRATION_SCALING \ - / TIME_SCALING - path_to_kmodel = './../../models/kin_varma.yml' path_for_output = './paramter_pop_{}.h5' @@ -118,27 +108,24 @@ Finally we call the parameter sampler for each sampled reference steady-state. T fluxes = [] for i, sample in samples.iterrows(): - # Load fluxes and concentrations + # Load fluxes and concentrations and scale them to the desired units fluxes = load_fluxes(sample, tmodel, kmodel, density=DENSITY, ratio_gdw_gww=GDW_GWW_RATIO, concentration_scaling=CONCENTRATION_SCALING, time_scaling=TIME_SCALING) + concentrations = load_concentrations(sample, tmodel, kmodel, concentration_scaling=CONCENTRATION_SCALING) - # Calibrate the glucose transporter to have a Vmax of 10mmol/gDW/h - pep_c = concentrations['pep_c'] - pyr_c = concentrations['pyr_c'] - g6p_c = concentrations['g6p_c'] - # Fetch equilibrium constants + # Fetch equilibrium constants and scale them to the desired units load_equilibrium_constants(sample, tmodel, kmodel, concentration_scaling=CONCENTRATION_SCALING, in_place=True) - # Generate sampels and fetch slowest and fastest eigenvalues + # Generate samples and fetch slowest and fastest eigenvalues params, lamda_max, lamda_min = sampler.sample(kmodel, fluxes, concentrations, only_stable=True, min_max_eigenvalues=True) diff --git a/doc/dynamic_simulations.rst b/doc/dynamic_simulations.rst index 8b13789..6de65f2 100644 --- a/doc/dynamic_simulations.rst +++ b/doc/dynamic_simulations.rst @@ -1 +1,78 @@ +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. + + +.. code-block:: python + + 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. + +.. code-block:: python + + 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' + ) diff --git a/doc/metabolic_control_analysis.rst b/doc/metabolic_control_analysis.rst index 8b13789..58d9f7d 100644 --- a/doc/metabolic_control_analysis.rst +++ b/doc/metabolic_control_analysis.rst @@ -1 +1,77 @@ +Metabolic control analysis +=========================== + + +Further we need to compute the steady state concentrations and fluxes for the steady-state +we aim to perform the MCA for. Here we use the concentrations and fluxes +from we used for the respective ORACLE parametrization. + +.. code-block:: python + + 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) + ref_fluxes = load_fluxes(ref_solution, tmodel, kmodel, + density=DENSITY, + ratio_gdw_gww=GDW_GWW_RATIO, + concentration_scaling=CONCENTRATION_SCALING, + time_scaling=TIME_SCALING) + + # Extract the current parameter set from the model + parameter_values = {p.symbol:p.value for p in kmodel.parameters.values()} + parameter_values = ParameterValues(parameter_values, kmodel) + + +With a set of parameters, fluxes and concentration we can then calculate the control coefficients. +To perform metabolic control analysis the control-coeffcient expressions +need to be compiled by calling ``kmodel.prepare()`` and ``kmodel.compile_mca(parameter_list=parameter_list)``, +where ``parameter_list`` is a TabDict contaning the symbols of the parameters we +want to the control-coeffcients for. Here we calculate the control-coeffcients with respect +to all maximal enzyme rates ``V_max``. + + +.. code-block:: python + + from skimpy.utils.tabdict import TabDict + from skimpy.viz.controll_coefficients import plot_control_coefficients + + + # Compile mca with parameter elasticities with respect to Vmaxes + parameter_list = TabDict([(k, p.symbol) for k, p in + kmodel.parameters.items() if + p.name.startswith('vmax_forward')]) + + kmodel.compile_mca(sim_type=QSSA,ncpu=8, parameter_list=parameter_list) + + flux_control_coeff = kmodel.flux_control_fun(ref_fluxes, + ref_concentrations, + [parameter_values, ]) + + lac_control_coeff = flux_control_coeff.slice_by('sample',0).loc['LDH_D', :] + + lac_control_coeff.index = [v.replace('vmax_forward_','') + for v in lac_control_coeff.index ] + plot_control_coefficients(lac_control_coeff, + filename='lac_control_coeff.html', + backend='svg', + ) + diff --git a/doc/modal_analysis.rst b/doc/modal_analysis.rst index 8b13789..14f201e 100644 --- a/doc/modal_analysis.rst +++ b/doc/modal_analysis.rst @@ -1 +1,55 @@ +Modal analysis +=========================== +Modal analysis allows to get insight into the dynamics close the the steady state. +The modal matrix provides information on which eigenvalues contribute to the +dynamics of each concentration. Similar to MCA, Modal analysis is conducted around a steady-state +therefore modal analysis requieres to compute or import a valid set of steady-state concentrations: + +.. code-block:: python + + 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.analysis.modal import modal_matrix + from skimpy.core.parameters import ParameterValues + from skimpy.utils.namespace import * + from skimpy.utils.tabdict import TabDict + + from skimpy.viz.modal import plot_modal_matrix + + # 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) + + parameter_values = {p.symbol:p.value for p in kmodel.parameters.values()} + parameter_values = ParameterValues(parameter_values, kmodel) + + + +To compute the modal-matrix the model needs to have compiled Jacobian expressions, they are build by calling ``kmodel.prepare()`` and ``kmodel.compile_jacobian()``. + +.. code-block:: python + + kmodel.prepare() + kmodel.compile_jacobian(sim_type=QSSA,ncpu=8) + + M = modal_matrix(kmodel,ref_concentrations,parameter_values) + + plot_modal_matrix(M,filename='modal_matrix.html', + plot_width=800, plot_height=600, + clustered=True, + backend='svg', + )