Merge pull request #75 from yaml2sbml-dev/develop

Version 0.2.0

.github/workflows/ci_push.yml

@@ -7,7 +7,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.7, 3.8]
+        python-version: [3.6, 3.7, 3.8, 3.9]
 
     steps:
     - name: Check out repository

.gitignore

@@ -8,24 +8,39 @@ yaml2sbml/tests/test_yaml2sbml/sbml_test.xml
 ### Files generated when running the notebooks (AMICI + PEtab files)
 
 ## Lotka Volterra Python
-./doc/examples/Lotka_Volterra_python/amici_models
-./doc/examples/Lotka_Volterra_python/amici_models/*
+./doc/examples/Lotka_Volterra/Lotka_Volterra_python/amici_models
+./doc/examples/Lotka_Volterra/Lotka_Volterra_python/amici_models/*
 
-./doc/examples/Lotka_Volterra_python/Lotka_Volterra_AMICI
-./doc/examples/Lotka_Volterra_python/Lotka_Volterra_AMICI/*
+./doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_AMICI
+./doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_AMICI/*
 
-./doc/examples/Lotka_Volterra_python/Lotka_Volterra_basic.xml
+./doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_basic.xml
 
-./doc/examples/Lotka_Volterra_python/Lotka_Volterra_PEtab/*
-!./doc/examples/Lotka_Volterra_python/Lotka_Volterra_PEtab/measurement_table.tsv
+./doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_PEtab/*
+!./doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_PEtab/measurement_table.tsv
 
 
 ## Lotka Volterra CLI
-./doc/examplesLotka_Volterra_CLI/Lotka_Volterra_PEtab/*
-./doc/examplesLotka_Volterra_CLI/Lotka_Volterra_PEtab.yml
-./doc/examplesLotka_Volterra_CLI/Lotka_Volterra_PEtab/*.xml
+./doc/examples/Lotka_Volterra/Lotka_Volterra_CLI/Lotka_Volterra_PEtab/*
+./doc/examples/Lotka_Volterra/Lotka_Volterra_CLI/Lotka_Volterra_PEtab.yml
+./doc/examples/Lotka_Volterra/Lotka_Volterra_CLI/Lotka_Volterra_PEtab/*.xml
 
 
+## Lotka Volterra Model Editor
+./doc/examples/Lotka_Volterra/Lotka_Volterra_Model_Editor/*.yaml
+./doc/examples/Lotka_Volterra/Lotka_Volterra_Model_Editor/*.xml
+./doc/examples/Lotka_Volterra/Lotka_Volterra_Model_Editor/Lotka_Volterra_PEtab
+
+
+## Sorensen example
+./doc/examples/Sorensen/*.xml
+./doc/examples/Sorensen/amici_models
+
+
+## FSP examples
+./doc/examples/Finite_State_Projection/amici_models
+./doc/examples/Finite_State_Projection/*.xml
+
 ## Format Features
 
 ./doc/examples/Format_Features/__pycache__

doc/examples/Finite_State_Projection/utils.py

@@ -0,0 +1,121 @@
+import amici
+import amici.plotting
+
+import importlib
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+import sys
+
+
+def plot_AMICI(sbml_dir: str,
+               t: np.ndarray,
+               r_max: int,
+               p_max: int):
+    """
+    Compiles, simulates and plots the AMICI model for the FSP example.
+    """
+    n_t = len(t)
+    # compile and Simulate model
+    simulation, model = compile_and_simulate(sbml_dir,
+                                             t,
+                                             r_max,
+                                             p_max)
+
+    def get_obs_idx_by_id(obs_id: str):
+        try:
+            return model.getObservableIds().index(obs_id)
+        except ValueError:
+            raise IndexError(f'No observable with id {obs_id}')
+
+    for i in range(n_t):
+
+        r_marginal = [simulation.y[i, get_obs_idx_by_id(f'x_r{r}')] for r in range(r_max)]
+        p_marginal = [simulation.y[i, get_obs_idx_by_id(f'x_p{p}')] for p in range(p_max)]
+        # rna
+        plt.subplot(n_t, 2, 2*i+1)
+        plt.fill_between(np.arange(r_max), 0, r_marginal, facecolor='blue', alpha=0.5)
+        plt.plot(r_marginal)
+        plt.ylabel(f't={int(t[i])}')
+        plt.yticks([])
+        plt.xlim(0, r_max-1)
+        plt.ylim(0, 0.3)
+
+        if i < n_t-1:
+            plt.tick_params(labelbottom=False)
+        else:
+            plt.xlabel('mRNA abundance')
+
+        # protein
+        plt.subplot(n_t, 2, 2*(i+1))
+        plt.fill_between(np.arange(p_max), 0, p_marginal, facecolor='blue', alpha=0.5)
+        plt.plot(p_marginal)
+        plt.yticks([ ])
+        plt.xlim(0, p_max-1)
+        plt.ylim(0, 0.2)
+        if i < n_t-1:
+            plt.tick_params(labelbottom=False)
+        else:
+            plt.xlabel('protein abundance')
+
+    plt.subplots_adjust(wspace=0.07, hspace=0.1)
+    plt.show()
+
+
+def compile_and_simulate(sbml_dir: str,
+                         t: np.ndarray,
+                         r_max: int,
+                         p_max: int):
+    """
+    Utility function, that compiles and simulates the FSP model.
+    """
+    model_name = 'Gene_regulation_FSP'
+
+    # define marginal rna concentrations
+    observables_rna = {f'x_r{r}': get_marginal_rna_obs(r, p_max)
+                       for r in range(r_max)}
+
+    # define marginal protein concentrations
+    observables_protein = {f'x_p{p}': get_marginal_protein_obs(p, r_max)
+                           for p in range(p_max)}
+
+    observables = {**observables_rna, **observables_protein}
+
+    sbml_importer = amici.sbml_import.SbmlImporter(sbml_dir)
+    sbml_importer.sbml2amici(
+        model_name=model_name,
+        output_dir=os.path.join('amici_models', model_name),
+        observables=observables)
+
+    # Import the compiled model.
+    sys.path.insert(0,
+                    os.path.abspath(os.path.join('amici_models', model_name)))
+    model_module = importlib.import_module(model_name)
+    model = model_module.getModel()
+
+    # Simulate.
+    model.setTimepoints(t)
+    solver = model.getSolver()
+    simulation = amici.runAmiciSimulation(model, solver)
+
+    return simulation, model
+
+
+def get_marginal_rna_obs(r: int, p_max: int):
+    """Get the observable for a marginalized RNA abundance."""
+    marginal = ''
+    for p in range(p_max-1):
+        marginal += f'x_{r}_{p} + '
+    marginal += f'x_{r}_{p_max-1}'
+
+    return {'name': f'x_r{r}', 'formula': marginal}
+
+
+def get_marginal_protein_obs(p: int, r_max: int):
+    """Get the observable for a marginalized protein abundance."""
+    marginal = ''
+    for r in range(r_max-1):
+        marginal += f'x_{r}_{p} + '
+    marginal += f'x_{r_max-1}_{p}'
+
+    return {'name': f'x_p{p}', 'formula': marginal}

doc/examples/Lotka_Volterra/Lotka_Volterra_Model_Editor/Lotka_Volterra_Model_Editor.ipynb

@@ -0,0 +1,205 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Model Editor\n",
+    "\n",
+    "This notebook demonstrates functionality of the `yaml2sbml` model editor using the Lotka Volterra equations as an example. The \"Lotka-Volterra\" equations are given by \n",
+    "\n",
+    "\\begin{align*}\n",
+    "\\frac{d}{dt} x_1 &= \\alpha x_1 - \\beta x_1x_2, \\\\\n",
+    "\\frac{d}{dt} x_2 &= \\delta x_1x_2 - \\gamma x_2.\n",
+    "\\end{align*}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## ODE model\n",
+    "\n",
+    "We first generate a basic model, that only contains all necessary information for generating an SBML file for model simulation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from yaml2sbml import YamlModel\n",
+    "\n",
+    "# generate model\n",
+    "model = YamlModel()\n",
+    "\n",
+    "# add ODEs\n",
+    "model.add_ode(state_id='x_1', \n",
+    "              right_hand_side='alpha*x_1 - beta*x_1*x_2', \n",
+    "              initial_value=2)\n",
+    "model.add_ode(state_id='x_2', \n",
+    "              right_hand_side='delta*x_1*x_2 - gamma*x_2', \n",
+    "              initial_value=2)\n",
+    "\n",
+    "# add parameters\n",
+    "model.add_parameter(parameter_id='alpha', nominal_value=2)\n",
+    "model.add_parameter(parameter_id='beta', nominal_value=4)\n",
+    "model.add_parameter(parameter_id='gamma', nominal_value=3)\n",
+    "model.add_parameter(parameter_id='delta', nominal_value=3)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "`yaml2sbml` can export the `model` object either to YAML or to SBML directly, via  "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# write to YAML\n",
+    "model.write_to_yaml('Lotka_Volterra_basic.yaml', overwrite=True)\n",
+    "\n",
+    "# write to SBML\n",
+    "model.write_to_sbml('Lotka_Volterra_basic.xml', overwrite=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "There are further functions to: \n",
+    "* get all `parameter_id`s via `model.get_parameter_ids()` \n",
+    "* get a parameter by its id (`model.get_parameter_by_id('alpha')`) \n",
+    "* delete a parameter by its id (`model.delete_parameter('alpha')`)\n",
+    "\n",
+    "Similar functions also exist for the other model components."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Parameter Estimation Problem\n",
+    "\n",
+    "Now we want to extend the current `model` to include all the necessary information for parameter estimation in PEtab. Therefore we load the model from the `.yaml` file and modify the parameters, such that it also contains all information that is going to be written into the PEtab parameter table."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model = YamlModel.load_from_yaml('Lotka_Volterra_basic.yaml')\n",
+    "\n",
+    "# extend parameters \n",
+    "model.add_parameter(parameter_id='alpha',\n",
+    "                    nominal_value=2,\n",
+    "                    parameter_scale='log10',\n",
+    "                    lower_bound=0.1,\n",
+    "                    upper_bound=10,\n",
+    "                    estimate=1, \n",
+    "                    overwrite=True)\n",
+    "\n",
+    "model.add_parameter(parameter_id='beta',\n",
+    "                    nominal_value=4,\n",
+    "                    parameter_scale='log10',\n",
+    "                    lower_bound=0.1,\n",
+    "                    upper_bound=10,\n",
+    "                    estimate=1, \n",
+    "                    overwrite=True)\n",
+    "\n",
+    "model.add_parameter(parameter_id='gamma',\n",
+    "                    nominal_value=3,\n",
+    "                    parameter_scale='log10',\n",
+    "                    lower_bound=0.1,\n",
+    "                    upper_bound=10,\n",
+    "                    estimate=1,\n",
+    "                    overwrite=True)\n",
+    "\n",
+    "model.add_parameter(parameter_id='delta',\n",
+    "                    nominal_value=3,\n",
+    "                    parameter_scale='log10',\n",
+    "                    lower_bound=0.1,\n",
+    "                    upper_bound=10,\n",
+    "                    estimate=1,\n",
+    "                    overwrite=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A parameter fitting problem in PEtab allows for the specification of observables and experimental conditions:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# specify an observable:\n",
+    "model.add_observable(observable_id='prey_measured', \n",
+    "                     observable_formula='log10(x_1)',\n",
+    "                     noise_formula='noiseParameter1_prey_measured', \n",
+    "                     noise_distribution='normal',\n",
+    "                     observable_transformation='lin')\n",
+    "\n",
+    "# specify trivial condition\n",
+    "model.add_condition(condition_id='condition1', \n",
+    "                    condition_dict={})"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The modified model can either be exported to YAML or PEtab via"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# write to YAML\n",
+    "model.write_to_yaml('Lotka_Volterra_PEtab.yaml', overwrite=True)\n",
+    "\n",
+    "# write to PEtab\n",
+    "model.write_to_petab(output_dir='./Lotka_Volterra_PEtab',\n",
+    "                     model_name='Lotka_Volterra',\n",
+    "                     petab_yaml_name='Lotka_Volterra_problem', \n",
+    "                     measurement_table_name='measurement.tsv')"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}