Demo

Here we show how to implement a genome-scale biosimulation using our Matlab API.

At a glance

In this brief document, we will describe the following (complete, with some variations) workflow:

import Netabolics.*

% Select a model
Cell = Biology.Template.Model("generic");

% Select some phenotypes 
Healthy = Biology.Template.Phenotype(Cell,"normal");
Diseased = Biology.Template.Phenotype(Cell,"cancer");

% Simulate healthy-disease transition
Scenario = Biology.Simulation(Cell,[Healthy,Diseased]);
Scenario.solve();

% Display the 10 most affected species
Scenario.Statistics.Species(1:10)

% Deconvolve the disease-healthy transition
Task = AI.Agent(Diseased,Healthy);
Task.search();

% Display the best target combination
Task.Statistics.Combination(1)

A structure as simple as the above can be the starting point for more advanced analyses.

1. Initialize the toolbox

IMPORTANT: This will clear the current workspace, so it is recommended to include the following operations in your startup.m file:

cd("/home/user/Netabolics/");
init__();

Let's add the Netabolics namespace to the current import list:

import Netabolics.*

and we are ready to go.

2. Select a model

For our purposes, a model is a cellular network reconstruction. In principle, a model could also represent a multicellular system, but for now we work at the single-cell level. Often a model is not to be interpreted as a real individual entity (one single cell), but rather as a population of cells of the same type (homogeneous population), or in some cases as a population of cells of different type (heterogeneous population) from a specific tissue (e.g., something like an "average" cell). The level of detail depends on the datasets at hand, for example bulk RNAseq vs scRNAseq.

While it is possible to use third-party SBML models, we have developed our proprietary and manually curated models. The toolbox is guaranteed to work correctly only with our curated models.

Let's use our latest human cell template model:

Cell = Biology.Template.Model("generic-10.22.89");

The model version (10.22.89) indicates the approximate number of genes (10K), molecular species including proteins and metabolites (22K), and gene-associated chemical reactions (89K).

The biological network models we use can be visualized as a directed hypergraph using:

Cell.plot(EdgeAlpha=0.2,Layout="force",UseGravity=true);

which produces the following figure:

Here species nodes are colored orange and reaction nodes/edges are colored blue.

3. Define phenotypes

3.1 Using experimental data

The above template model comes with a default parametrization (e.g., species concentrations, reaction rates, etc). In order to be of any use, we need to associate a phenotype (i.e., a state) to the model. This step could be seen as a "differentiation" into a cell type and/or a given condition.

To start with, we define a healthy generic phenotype for the model Cell:

Healthy = Biology.Template.Phenotype(Cell);

In the absence of data, this default template phenotype can already be used (e.g., for testing purposes). Of course, it is possible, and advisable, to load experimental datasets (e.g., RNAseq):

Healthy.loadDataset("HealthySample.dat");

Similarly, we can load experimental datasets from diseased and treated patients (here using the Biology.Phenotype constructor directly, instead of loading the dataset into the template phenotype as done above):

Diseased = Biology.Phenotype(Cell,"DiseasedSample.dat");
Treated = Biology.Phenotype(Cell,"TreatedSample.dat");

3.2 Using prebuilt templates

In the absence of data, and also to see some capabilities of the toolbox that are used for target deconvolution, let's see an example of defining phenotypes using prebuilt templates.

In this example, we use prebuilt phenotypes from the OncoDB public database (https://oncodb.org/) for a specific tissue as the healthy and diseased (here colon adenocarcinoma, COAD) states:

Healthy = Biology.Template.Phenotype(Cell,"Tissue/OncoDB/DGE/COAD/normal");
Diseased = Biology.Template.Phenotype(Cell,"Tissue/OncoDB/DGE/COAD/cancer");

To implement a treatment, we can start from prebuilt drugs (here as a combination) and then define a treated phenotype by converting the therapy to a new phenotype:

Therapy = Biology.Template.Drug(["aspirin","genistein"]);
Treated = Therapy.phenotype(Cell);

Overall, at the end of the above sequence of steps, we have the following:

• Model (Cell), created using a prebuilt template model.
• Healthy phenotype (Healthy), created using experimental data from normal tissue.
• Diseased phenotype (Diseased), created using experimental data from cancer tissue.
• Treated phenotype (Treated), created using drug combination.

With these in place, we can now proceed to the simulation phase. Notice, however, that depending on the specific goals, not all of the above phenotypes might be necessary (e.g., if we just want to see how a specific state looks like).

3.3 Using manual perturbations

Although very inefficient for large-scale models, manual perturbations can be useful for quick tests (e.g., during hypothesis generation).

The toolbox allows for the following perturbations to be handled manually:

• Absolute or relative changes of gene expression (e.g., mRNA levels). These are typically derived from transcriptomics data.
• Absolute or relative changes in species levels (e.g., protein abundance or metabolite concentration). These are typically derived from proteomics or metabolomics data.
• Relative changes in gene accessibility. These are typically derived from epigenomics data.
• Relative changes in protein function (e.g., gain-of-function or loss-of-function mutations). These are typically derived from genomics data together with a variety of experimental functional studies.

In addition, the following perturbations are available but extremely difficult to handle (not recommended):

• Absolute changes of protein activity rates (e.g., enzyme velocity). These are typically derived from fluxomics data.
• Absolute changes of gene expression rates (e.g., velocity of gene transcription). In principle, these could be derived from time-resolved transcriptomics data.

Detailing manual perturbations is beyond the scope of the present demonstration. Here is a simple illustration:

% Initialize an empty phenotype
Custom = Biology.Phenotype(Cell);

% Manually include a perturbation
Custom.addPerturbation("PTGS1","LOFM",0.5);

The above will result in a 50% loss-of-function mutation (LOFM) for PTGS1-encoded protein prostaglandin-endoperoxide synthase 1.

4. Simulate cellular state transitions

We can simulate the temporal evolution across phenotypes using our proprietary methods for creating and solving a dynamical system of ordinary differential equation (ODE) starting from the static model reconstruction and data. Of course, in real conditions the disease phenotype does not occur at-once as we do here (a more realistic picture would be possible if we had the temporal evolution of the diseased phenotype, which would require data to be acquired across the disease evolution time-span, something that however is seldom available in the practice).

It is possible to simulate as many phenotypes as we need in an ordered sequential manner. Here we simulate our three phenotypes, which translates into two transitions (healthy→diseased and diseased→treated):

% Prepare the ODE system (with default settings)
Scenario = Biology.Simulation(Cell,[Healthy,Diseased,Treated]);

% Ensemble solving 
Scenario.solve(); 

We can now plot the results, for example the top 40 species with largest responses (these are automatically sorted for variance):

figure(Color="w");

for i = 1:40
    subplot(4,10,i);
    Scenario.plot(i);
end

This will produce the following figure (notice that by default each phenotype is simulated for 24 hours):

We note that the first transition at 24 h (healthy→diseased) is substantially more impactful than the second transition at 48 h (diseased→treated), which implies that the treatment we designed has almost no effect.

We can obtain information on individual elements by visiting relevant online resources. For instance, to visit online resources for FABP1 (top left in the above figure) we do:

FABP1 = Cell.SBML.species(Scenario.Statistics.species(1).Index);
URL = FABP1.cvterms(1).resources{:}; % 'https://identifiers.org/uniprot:P07148'

% Open URL in a web browser window
web(URL); 

The simulation results can be quickly examined by using predefined cellular processes of physiopathological significance. As an illustration, here we generate a table with the largest changes from healthy to diseased phenotype (first transition):

Table = splitvars(struct2table(Scenario.Statistics.process));
Table.Change = Table.Mean_2-Table.Mean_1;

sortrows(Table,"Change","descend")

which will yield the following table (first 5 rows; not all columns are shown):

| Index | Symbol | Change | ... | Units | | - | - | - | - | - | | 2129 | FATTY ACID OXIDATION | 140.67 | ... | "AU"
| 2001 | ACTIN CYTOSKELETON REORGANIZATION | 130.45 | ... | "AU"
| 2159 | GRANULOCYTE DIFFERENTIATION | 81.482 | ...| "AU"
| 2389 | TRANSLATIONAL REGULATION | 70.218 | ... | "AU"
| 2107 | DNA REPAIR | 66.312 | ... | "AU"
| ... | ... | ... | ... | ...

Notably, analysing changes in many thousands molecular species requires advanced statistical methods, which is beyond the scope of the present demonstration.

5. Run target deconvolution

Now let's finally use AI to find the best treatment for a disease given a set of therapeutic op