Tutorials

Basic tutorials

This is a short introduction to using the BMA based on a set of demos. Please add new examples of how to do common tasks and interpret results

Building models

The BMA is designed to allow users to rapidly build networks of gene interactions based on the available literature using a standard model type, the qualitative network. The most simple way to develop models is simply to draw the networks as you understand them, adding proteins, genes and phenotypes and their relationships. By default Boolean networks are built. This means that each variable (that is, protein or gene) can occupy two states and the default target function (which determines how variables change over time) is the average of the activators - the average of the inhibitors. Once the gene interactions are built in the GUI, the user can increase the granularity of individual variables in the network (the number of states which a variable can occupy) and the target function.

Simulation and stability analysis

The BMA presently allows the system to be analysed using two different approaches:

1. Simulation- where the changes from a single state are observed as a function of (discrete) time

2. Stability analysis- where the end points of all possible simulations are queried

These two type of analysis treat the QN as synchronous. This means that all variables update with each step of the simulation. One common alternative is asynchronous updates, where a single variable is updated (non-deterministically) at each timestep.

Stability analyses ask the question- do all initial states of a simulation end up in the same final state, or do they end up in multiple states, or in a loop (i.e. infinitely alternating through a series of states)? If it ends in a single state, the system is said to be stable- if not, the "Further testing" button can be used to show how the system is unstable (either by finding the multiple end states or the loop). For biologists, stability can be thought of as a measure for robustness- if you expect a system to reliably end up in a single state then this is a powerful way to test for this. Examples of this could be metabolic networks, or signalling systems which respond reliably, or chemical networks at a robust equilibrium. The stability progression graph shows how the algorithm has proved stability- it works by iteratively reducing the range of individual variables until they can no longer be reduced.

The demos

Modelling Dicty cAMP signalling

Dictyostelium Discoidum us a single celled organism which undergoes a transition under starvation conditions to become a multicellular "slug". Starvation initially causes aggregation by cells starting synchronous oscillations in cAMP levels. Other cells detect this over a distance, join the cycles and chemotax towards the cAMP concentrations. In non-starving cells, externally applied cAMP oscillations are sufficient to cause this aggregation behaviour.

This network has been modelled previously using a synchronous chemical reaction network and found to generate stable oscillations. However, it was subsequently suggested that this property was highly dependent on the choice of parameters, and that an asynchrous model showed more robust cycles.

We can build the network as above to test the proposal that a synchronous network generates robust oscillations. Starting as a Boolean network, the model can be built as drawn above. To test the behaviour we can use stability analysis. The behaviour which the model must reproduce (the "specification") is the single cell must show an oscillation in the cAMP levels, and have a single fix point (representing the non-starvation state).

As can be seen from the model, it can be built almost identically to the image above. The key difference is that we explicitly include a constant +1 input to RegA as a variable. Testing the model, we find that the system has a single fix point (by simulating from an all zero state) and an oscillation in cAMP.

We can further modify the model to include multiple cells, which share a pool of external cAMP. Testing shows that these cells oscillate synchronously, demonstrating the biological relevance

Which network is correct?>

These three diagrams represent three possible networks of a phosphorylation cascade. Biologically it is known that if the cell is in the presence of a low concentration of a ligand I, the level of phosphorylation of X is low. If I increases suddenly, X phosphorylation temporarily increases before dropping down to the same level as when the concentration of I was low. You know from the literature one of these is probably right, but want to use modelling to find out which.

Step 1: create a specification

1. When I is low, X is low. This is a robust result, so the system is also stable.

2. When I is high, X is low. This system is also stable

3. When I increases, X transiently increases

Step 2: Build each model and test the specifications

1. Build each model as above, with the default target functions

2. Test each model in turn according to the three specifications above. To mimic the high and low I states, the target function of I can be set to 1 or 0 respectively

3. The first model passes the first specification (tested using stability analysis), but fails the second (tested using stability analysis)- when I is high, X is high

4. The second model also passes the first specification, but when the second spec is tested a cycle is observed and so the model fails (shown using further testing in stability analysis)

5. When tested using stability analysis both the first and second specification are passed. The third specification can be tested using simulation, starting from the low I stable state (with the target function of I set to 1 i.e. high). This passes as X transiently becomes 1

Nitrate and Nitrite sensing in E. coli- Building and understanding models with a tight specification

Nitrite and nitrate sensing in E. coli is achieved through receptors and two downstream regulators, NarQ and NarX, and NarP and NarL. Both NarQ and NarX respond to both ligands, and both activate the downstream effectors NarP and NarL. However, the receptors respond differently to the ligands and have complex relationships with the effectors. The network of interactions is known

Based on some of the early work (http://www.ncbi.nlm.nih.gov/pubmed/8501030), I aimed to build a model showing how the protein signals combine to control the expression of a single operon, aeg-46.5. The table above is taken from the paper, and translated into a formal specification as below

Genotype		NarL+NarP+			NarL NarP+			NarL+NarP			NarL NarP
		None	NO3	NO2	None	NO3	NO2	None	NO3	NO2	None	NO3	NO2
NarX+	NarQ+	0	1	2	2	4	4	0	0	0	0	0	0
NarX	NarQ+	0	1	0	2	4	4	0	0	0	0	0	0
NarX+	NarQ	0	0	0	2	4	3	0	0	0	0	0	0
NarX	NarQ	0	0	0	0	0	0	0	0	0	0	0	0

The protein interaction network was taken from more recent work (http://www.ncbi.nlm.nih.gov/pubmed/25873583), giving a model which looks like the below.

In building the model a few arbitrary decisions were made initially- all proteins had a granularity of 5 (i.e. a range of 0-4) and didn't interact with themselves.

Several things can be noted from looking at the specification and the paper.

1. NarP is required for any aeg activity, whilst loss of NarL increases aeg activity uniformly- hence NarP activates aeg whilst NarL inhibits aeg.
   

2. NarX and NarQ appear to respond with different sensitivities to Nitrite and Nitrate, with NarQ responding more strongly to Nitrite than NarX
   

3. Nitrate/Nitrite=2 represents the presence of the ligand
   

4. The model is expected to be stable- the responses are implied to be reliable
   

To build the model, the "NarP+ NarL " genotype was considered first. In the absence of nitrite or nitrate, and if at least one of NarX or NarQ is expressed, basal activity from the receptors activates NarP. The specific function for NarP to NarX and NarQ was then built up from the single mutant data, then tested against the NarX+ NarQ+ genotype to give the functions for different proteins

NarP:

min(1,var(NarX)+var(NarQ))*2 + max(0,(var(NarQ)-1))*2 + max(0,(var(NarX)-1))

Interpreted as:

Basal activity from receptors + activation from NarQ + activation from NarX

NarX:

1 + var(Nitrate) + max(0,(var(Nitrite)-1))

Interpreted as:

Basal activity + Nitrate + Weaker response to Nitrite

NarQ:

1 + var(Nitrite) + max(0,(var(Nitrate)-1))

Interpreted as:

Basal activity + Nitrite + Weaker response to Nitrate

The next step was to consider the wild type. Given the model describes the activity of the system without the activity of NarL, the contribution of NarL to aeg activity can be considered a specification in isolation and a target function proposed which can generate that specification. This is a complex formula and should be considered carefully in terms of how it can be interpreted- this is how I broke it down. n means that the value can be anything

`

	X	Q	P	L	min(1,var(NarX)) + min(1,var(NarQ))	max((var(NarX)-var(NarQ)),(var(NarQ)-var(NarX)))	-min(1,var(NarX))*min(1,var(NarQ))*min(1,max(0,(var(NarQ)-var(NarX))))	Total
None	1	1	2	2	2	0	0	2
NO3	3	2	4	3	2	1	0	3
NO2	2	3	4	2	2	1	1	2
KO+None	0	1	2	2	1	1	0	2
KO+None	1	0	2	2	1	1	0	2
KO+None	0	0	0	n	0	0	0	0
KO+NO3	0	2	4	3	1	2	0	3
KO+NO3	3	0	4	4	1	3	0	4
KO+NO3	0	0	0	n	0	0	0	0
KO+NO2	0	3	4	4	1	3	0	4
KO+NO2	2	0	3	3	1	2	0	3
KO+NO2	0	0	0	n	0	0	0	0
`

This function is ultimately therefore

min(1,var(NarX)) + min(1,var(NarQ)) + max((var(NarX)-var(NarQ)),(var(NarQ)-var(NarX))) - min(1,var(NarX))*min(1,var(NarQ))*min(1,max(0,(var(NarQ)-var(NarX))))

Interpreted as

The basal activity contribution of NarX and NarQ + antagonistic activation by NarQ/NarX (the difference determines the signal strength) - modifier to reflect stronger response to NarX

I found this by building a table as above and testing out different function components, whilst speculating on their meaning. The more complex a function is, the more important it is to be able to justify its complexity. It may be worth noting that if this was rendered as a truth table, the break down of the function into its interpreted behavior may not have been possible and thus the function harder to justify.

This model could be further built out to include the responsiveness to other operons- its worth noting that this set of functions gives no guarantees that the others will work trivially. This may require refinement of the either functions or variable granularity, or even replacement of the functions with ones which work more generically. This would need to be found by experimentation.

Finding missing interactions in skin model signalling

Early models of epithelial development were based around the concept that Wnt and Notch activity varied across 5 layers of skin cells, creating a gradient in gene activities and defining the cell development process. Note that this model was built on the best available knowledge at the time, and may have subsequently been superseded by recent findings about the stem cell populations of the epithelia- this is not relevant for understanding this model and its behavior. It is further worth noting that the subsequent prediction from working with this model was verified experimentally.

The skin exists at homeostasis, and the layers of the skin are expected to be stable with a gradient in Jagged activity across the layers- specifically we expect to see more activity in the bottom layer. The model was built based on known gene interactions, with a granularity of 5. A 1D model of the 5 layers was initially made, with one cell representing each layer.

Despite the model including all known gene interactions, the model was unstable. Subsequent reading of the literature suggested that in other systems, Jagged could be generated by B-Catenin, via its gene target (labelled GT1). Including this interaction (Skin1D.json) gives a stable model where the level of Jagged activity is highest in the basal cells. Alternative single interactions also stabilise the system, but crucially they don't generate the gradients expected.

Later work extended the models into 2D and 3D, identifying other bugs in the model (http://research.microsoft.com/pubs/139258/VMCAI_11.pdf).