Proximo-Distal patterning

In this tutorial we will go through the steps to derive the equations representing a Gene Regulatory Network (GRN) model. We will create a model that simulates a Proximal-Distal (PD) patterning process, using Meis, Hoxa11 and Hoxa13 as markers of stylopod, zeugopod, and autopod, respectively.

The model will be based on the published network described as "Model F" in the following publication: Uzkudun et al., 2015.

Our model can be represented using the network illustrated below:

Proximo-Distal patterning model network scheme. Retinoic Acid (R) and FGF (F) are the proximal and distal signals, respectively. Meis (M), Hoxa11 (H11) and Hoxa13 (H13) are the predicted patterns of the model.

In this tutorial we will first introduce some general concepts of the modelling process, and of the type of reaction-diffusion model that we will develop, and then we will go through the process of implementing the model in LimbNET.

General GRN modelling

Dynamics

Before writing the set of equations to describe the system, let's refresh some basics of gene regulatory network systems. In general, we can describe the dynamics of a given gene \(g\) according to the equation

\[ \frac{\partial g}{\partial t} = p_g - \lambda_g g + D_g\nabla^2g \]

where \(p_g\) is the production of \(g\), \(\lambda_g\) is the degradation rate of \(g\), and \(D_g\) is its diffusion coefficient. \(\frac{\partial g}{\partial t}\) simply means the rate of change of \(g\) with respect to time.

Modulation

Other genes, transcription factors, can modulate this overall production of \(g\) by upregulating (activating) or downregulating (inhibiting) it.

Upregulation (activation)

For instance, let's imagine that there is another gene \(a\) that activates \(g\), with a strength \(k_a\). For now, we will assume that all gene interactions can be modelled via sigmoidal (Hill) functions. In the case of activation, we introduce a function \(H_{up}(x, k, \mu)\), where \(x\) is the input gene, \(k\) is the strength coefficient of the link, and \(\mu\) is the steepness of the Hill function.

We can modulate our overall production rate according to this upregulation equation, proportional to \(a\), resulting in a final reaction-diffusion PDE for \(g\) expressed as

\[ \frac{\partial g}{\partial t} = p_gH_{up}(a, k_a, \mu) - \lambda_g g + D_g\nabla^2g \]

Downregulation (inhibition)

Now let's imagine that our activator \(a\) itself has an inhibitor, \(b\), with strength \(k_b\) and steepness \(\mu\). The resulting downregulation can be modelled using a \(H_{down}(b, k_b, \mu)\) function:

\[ \frac{\partial a}{\partial t} = p_aH_{down}(b, k_b, \mu) - \lambda_a a + D_a\nabla^2a \]

Implementation

Our final system of equations for our three gene products, \(g\), \(a\) and \(b\) would be

\[ \begin{align} \frac{\partial g}{\partial t} &= p_gH_{up}(a, k_a, \mu) - \lambda_g g + D_g\nabla^2g \\ \frac{\partial a}{\partial t} &= p_aH_{down}(b, k_b, \mu) - \lambda_a a + D_a\nabla^2a \\ \frac{\partial b}{\partial t} &= p_b - \lambda_bb + D_b\nabla^2b \end{align} \]

In LimbNET, we can use the built-in functions hillup for \(H_{up}(x, k, \mu)\) and hilldn for \(H_{down}(x, k, \mu)\) (see built-in functions for further details). Therefore, we could write our system of equations for the reactions only (note that the diffusion uses another input field) with the following notation:

dg/dt = p_g*hillup(a, ka, mu) - lam_g*g

da/dt = p_a*hilldn(b, kb, mu) - lam_a*a

db/dt = p_b - lam_b*b

Writing down the system of reactions

Ok! Now we are ready to tackle our model. We will start with only the reaction component (we will add the diffusion later).

For simplicity we will assume that the production and degradation rates of Cyp, Meis, Hoxa11, Hoxa13 are the same, so, we will declare the production parameters as \(p_{all}\) and degradation as \(\lambda_{all}\). All the Hill functions shall have the same steepness \(\mu\), except one - the Hoxa13 to Hoxa11 interaction.

Genes

Let's start with the stylopod marker Meis. This signal has one activator RA, with strength \(k_2\), so its equation would be

\[ \frac{\partial \mathrm{Meis}}{\partial t} = p_{all} H_{up}(\mathrm{RA}, k_2, \mu) - \lambda_{all}\mathrm{Meis} \]

Following the same logic, we can derive the equation for the autopod marker Hoxa13 as

\[ \frac{\partial \mathrm{Hoxa13}}{\partial t} = p_{all}H_{down}(\mathrm{RA}, k_5, \mu) - \lambda_{all}\mathrm{Hoxa13} \]

The zeugopod marker Hoxa11 has two inputs: FGF as activator, and Hoxa13 as inhibitor, with strengths \(k_6\) and \(k_4\), respectively. We can combine the effect of both regulators by multiplying their Hill functions - in this case, the combined effect is somewhat akin to a logical AND - ending up with the following equation:

\[ \frac{\partial \mathrm{Hoxa11}}{\partial t} = p_{all}H_{up}(\mathrm{FGF}, k_6, \mu)H_{down}(\mathrm{Hoxa13}, k_4, \mu^\prime) - \lambda_{all}\mathrm{Hoxa11} \]

Similarly, we can write the equation for Cyp as

\[ \frac{\partial \mathrm{Cyp}}{\partial t} = p_{all}H_{up}(\mathrm{FGF}, k_1, \mu)H_{down}(\mathrm{Meis}, k_9, \mu) - \lambda_{all}\mathrm{Cyp} \]

Retinoic Acid (proximal input)

Retinoic Acid (RA) is a molecule with different behaviour than the previous reactants, which represent genes.

For this model, we assume that RA is produced in the flank and passively diffuses. To express this, we can add a binary variable \(\mathrm{flank}_{in}\), dependent on our spatial co-ordinates \(\mathbf{x}\):

\[ \mathrm{flank}_{in}(\mathbf{x})= \begin{cases} 1, & \text{if } \mathbf{x} \in \text{flank} \\ 0, & \text{otherwise} \end{cases} \]

Thus, production of RA is 0 everywhere except within the flank. Moreover, RA is linearly inhibited by Cyp with strength \(c_1\), and has a degradation rate \(\lambda_{ra}\). And so the full equation for RA can be written as

\[ \frac{\partial \mathrm{RA}}{\partial t} = p_{RA} \mathrm{flank}_{in}(x) - c_1 \mathrm{Cyp} \cdot \mathrm{RA} - \lambda_{RA} \mathrm{RA} \]

FGF (distal input)

Finally, we have FGF, the distal input.

The degradation rate parameter for both FGFs is the same, and we denote it as \(\lambda_{FGF}\). The total production of FGF is the sum of \(\mathrm{FGF}_4\) and \(\mathrm{FGF}_8\), both of which have their own production rate: \(p_{fgf_4}\) and \(p_{fgf_8}\) respectively. Since both FGFs act as feed-forward inputs in this model, we define two new spatial variables to specify where the molecules are produced, \(\mathrm{FGF4}_{in}(x)\) and \(\mathrm{FGF8}_{in}(x)\). The values of these variables are derived from experimental data, and are pre-defined.

The final reaction equation for FGF is

\[ \frac{\partial \mathrm{FGF}}{\partial t} = p_{FGF_4} \mathrm{FGF_4}_{in}(x) + p_{FGF_8}\mathrm{FGF_8}_{in}(x) - \lambda_{FGF}\mathrm{FGF} \]

The complete PDE system

With the final input, FGF, now defined our system of six coupled reaction equations is complete. We can summarize the system as

\[ \begin{align} \frac{\partial \mathrm{FGF}}{\partial t} &= p_{FGF_4} \mathrm{FGF_4}_{in}(x) + p_{FGF_8}\mathrm{FGF_8}_{in}(x) - \lambda_{FGF}\mathrm{FGF} \\ \frac{\partial \mathrm{RA}}{\partial t} &= p_{RA} \mathrm{flank}_{in}(x) - c_1 \mathrm{Cyp} \cdot \mathrm{RA} - \lambda_{RA} \mathrm{RA} \\ \frac{\partial \mathrm{Cyp}}{\partial t} &= p_{all}H_{up}(\mathrm{FGF}, k_1, \mu)H_{down}(\mathrm{Meis}, k_9, \mu) - \lambda_{all}\mathrm{Cyp} \\ \frac{\partial \mathrm{Meis}}{\partial t} &= p_{all} H_{up}(RA, k_2, \mu) - \lambda_{all}\mathrm{Meis} \\ \frac{\partial \mathrm{Hoxa11}}{\partial t} &= p_{all}H_{up}(\mathrm{FGF}, k_6, \mu)H_{down}(\mathrm{Hoxa13}, k_4, \mu^\prime) - \lambda_{all}\mathrm{Hoxa11} \\ \frac{\partial \mathrm{Hoxa13}}{\partial t} &= p_{all}H_{down}(RA, k_5, \mu) - \lambda_{all}\mathrm{Hoxa13} \\ \end{align} \]

Now, that we've got the set of reactions, let's add them to our model...

Creating our model in LimbNET

Now that we have defined our full set of equations, we can start creating our model in LimbNET. The first step is to create a new model specification: select the "models" tab in the left-hand panel, and click the "new model" button in the top-right corner. A dialog will appear, prompting you to enter a name for your model.

New model creation dialog.

Once you have entered a name, click "create" to create the model.

You will be taken to the model specification view, where you can add parameters, variables and reactions to your model specification:

Model specification view for the new model.

For the time being we can leave the base mesh defined as-is - the default mouse limb morphomovie. Note that the model specification is not shared by default, it is only visible to you; should you wish to later share it you can toggle the switch labelled "public".

We will now move on to defining the parameters, variables and reactions for our model.

Since we have already written down our full system of equations, we know in advance which parameters and pre-defined patterns we will need in order for our equations to be valid. So we can proceed step-by-step in the order presented in the model editor interface:

1. Define our parameters,
2. Define our predefined patterns (system inputs),
3. Finally, define our reactions, which are dependent on the rest.

You may notice that as you are typing (or pasting in) the equations below, that the equation editor will sometimes indicate a validation error. Don't worry about this for now - the model is being validated against the equations as you type, and during this time you may be writing equations that depend on parameters or variables that you have not yet defined. Once the model is fully defined, it should be free of errors and will validate successfully.

Setting parameters

The first step in defining our model is to set the parameters. Analogous to writing computer code, parameters and variables must be defined somewhere in order for them to be used later in equations.

In our model, we have a number of parameters that we will need to define - 16 in total. So, under global parameters let's add all of them.

To add a parameter, simply click the "new parameter" button, and then enter a name for the parameter. The names of our parameters are as follows:

Name	Description	Value
p_all	Generic production rate (genes)	0.05
p_fgf4	Production rate of FGF4	0.0335
p_fgf8	Production rate of FGF8	0.005
p_ra	Production rate of RA	0.01
lam_all	Generic degradation rate (genes)	0.05
lam_fgf	Degradation rate of FGF	0.00704
lam_ra	Degradation rate of RA	0.0001
mu	Hill function steepness	1.25
mu_prime	Hill function steeepness for Hoxa13->Hoxa11 interaction	3.76
c1	Inhibition strength of RA by Cyp	2.46
k1	Midpoint of FGF activation of Cyp	0.353
k2	Midpoint of RA activation of Meis	0.145
k4	Midpoint of Hoxa13 inhibition of Hoxa11	0.622
k5	Midpoint of RA activation of Hoxa13	0.000699
k6	Midpoint of FGF activation of Hoxa11	0.06
k9	Midpoint of Meis inhibition of Cyp	1

Add all parameters in the above table to the model, and set their values as above.

parameter values

Why do the parameters take the values in the table above? To find out optimal values, the authors of the papers used an optimisation algorithm, to fit the model output to experimental data; here we will just use their published values directly.

At this point, your model specification should look something like this:

Model specification view with parameters added.

Save the model, by clicking "save" at the top of the model editor and let's move on to adding the predefined input patterns.

Adding predefined patterns

In our model, we will be using three predefined patterns, flank_in, fgf4_in and fgf8_in. In LimbNET, we can use predefined patterns to create a value that changes in time and space, but is not dependent on the model equations.

These patterns will all be defined in the mesenchyme. Therefore, in the model editor, scroll down to the "mesenchyme" section, and click the "new pattern" button. For each of the three patterns, you will need to enter a name, and select the pattern name from the dropdown.

For the first one type flank_in as the name, and select flank_region from the dropdown; the dropdown will also autocomplete, you can start typing in the name of the desired pattern and the list will be filtered for you.

Now, repeat this for the remaining two patterns:

Name	Pattern name
fgf4_in	fgf4exp__manu_presim
fgf8_in	fgf8exp__manu_presim

Finally, your list of predefined patterns should be complete:

Model specification view with predefined patterns added.

Save the model, by clicking "save" at the top of the model editor and let's move on to the final step: adding the reactions!

Adding the reaction equations to the model

As with the predefined patterns, all reactions that we are interested in will also take place in the mesenchyme, so make sure that you have scrolled down to the "mesenchyme" section of the model editor.

We will begin by adding the reactions for our two input molecules, RA and FGF.

Retinoic Acid

Let's add the RA.

In the "mesenchyme" section of the editor, click on the "new reaction" button to add a new reaction. We have already defined our flank region, \(\mathrm{flank}_{in}(\mathbf{x})\), as a named predefined pattern, which is available in our equation editor as the variable flank_in. Thus, we can just add the reaction equation for RA directly, as:

Name	Equation
ra	p_ra*flank_in - c1*cyp*ra - lam_ra*ra

Upon entering the equation, you will immediately see an error message appear, that the variable cyp is not defined. As you may be thinking already this is because we have not yet defined the variable cyp in our model; we will do this shortly, and this equation will then be valid.

FGF

Secondly, we add the FGF.

As with RA, the total FGF input is dependent on some predefined patterns, which we have already defined as fgf4_in and fgf8_in. We can click on the "new reaction" button again, and add the equation for FGF as:

Name	Equation
fgf	p_fgf4*fgf4_in + p_fgf8*fgf8_in - lam_fgf*fgf

Patterning genes (Cyp, Meis, Hoxa11, Hoxa13)

Finally, let's add the Cyp, Meis, Hoxa11 and Hoxa13 reactions.

You will have to add four new reactions, then input the name, and equation for each.

As stated before, we can write down the Hill functions using the keyword hillup or hilldn. While entering the equations, you may notice that the equation editor will sometimes indicate a validation error. Again, don't worry about this for the time being - you will see that the errors are generally related to incomplete parts of the model, or missing reactions, and will disappear as you add the remaining reactions and the model is revalidated.

The four reactions are outlined in the table below:

Name	Equation
cyp	p_all*hillup(fgf, k1, mu)*hilldn(meis, k9, mu) - lam_all*cyp
meis	p_all*hillup(ra, k2, mu) - lam_all*meis
hoxa11	p_all*hillup(fgf, k6, mu)*hilldn(hoxa13, k4, mu_prime) - lam_all*hoxa11
hoxa13	p_all*hilldn(ra, k5, mu) - lam_all*hoxa13

In summary, we should now have six reactions in the mesenchyme compartment: RA, FGF, Cyp, Meis, Hoxa11 and Hoxa13.

Model specification view with reactions added.

Note

Remember to regularly save your model as you add and change components!

Adding diffusion to the model

In this model, only FGF and RA diffuse, with coefficient values of 100, and 600 respectively.

We can add these values in the diffusion field in the section for each reactant.

Model specification view with diffusion added.

At this point, the model is complete, and ready for simulation. There should not be any more validation errors.

Simulation

After all of our hard work we can save the model and go to the simulation view by clicking on "simulation", and then "save and simulate" in the dialog that appears (if we had made no changes, and had already saved everything, we could instead just click "simulate only" to go directly to the simulation view without saving).

Once the simulation view opens, we can directly click "simulate" to run the simulation; an indicator will appear in the centre of the view, indicating simulation progress, and will disappear once the simulation is complete, allowing us to view the results.

Once the backend has returned a simulation result, you can select a reactant in one of the three colour channels' dropdowns and visualise it by clicking the play button or moving the frame slider:

Simulation view, results of the final model.

And with that, we have completed our model.

Congratulations, and happy modelling!