DOC: Expand documentation for ORACLE/MCA/MA and dynamic simulations

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)

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'
+                   )
 

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',
+                              )
+
 

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',
+                      )