Gray-Scott model
In this tutorial we are going to explore patterning arising through Reaction-Diffusion type mechanisms. We will use the classic Gray-Scott model, a Reaction-Diffusion model that can generate a wide range of static and dynamic patterns, based on combinations of just two parameters.
We will first write the model equations, the implementation, and then discuss how to simulate the model. Finally, we will explore some more advanced topics such as how to generate a Pearson plot (how to plot the model behaviour against two input parameters), and how to change the mesh used by our model.
Model Base
Navigate to the "models" tab from LimbNET's main (left) menu, and click on the "new model" button. A pop-up will appear with a text input where you can write the model name. Enter a suitable name and click "create"; you will be presented with the model editor view. By default, the model will be private until you specify otherwise.
Before starting the modelling, we recommended writing an informative description of the model. Of course this can also be updated at any point, as we add new features to the model. In this example, we may want to write something along the following lines:
model description
This is a basic implementation of the Gray-Scott Reaction-Diffusion model. Further information about the model can be found at, e.g., https://www.mrob.com/pub/comp/xmorphia/index.html.
Now we need to choose the base mesh/morphomovie for our model i.e., the mesh on which the simulation is going to run. You can find more information in morphomovies.
For this model we are going to use the square fixed size morphomovie, so select this in the top right-hand corner of the model editor:
With that set up, we are ready to move on to the model implementation!
The equations of the model can be written as
Where \(u\) and \(v\) are two interacting chemical species, and \(D_u\) and \(D_v\) are their diffusion coefficients, respectively.
The dynamics of the model are controlled by two key parameters: \(F\) - the feed rate, the rate at which new material is added to the system; and \(k\) - the kill rate, the rate at which \(v\) is removed from the system.
Since both rates are fixed global numerical values, we will define them as parameters, in the model's Global compartment.
As a starting point we can set up the parameters and corresponding values as k = 0.059 and F = 0.032, so that we will end up with nice stripes.
In the model editor, in the parameters section of the Global compartment, click on the add parameter button twice, and add the two parameters and their corresponding values:
We are now ready to input our equations into the model editor.
Since the model is running on the internal compartment of the mesh, we have to define the reactions in the Reactions (mesenchyme) section.
In LimbNET, the reaction and the diffusion are input in two different text fields, even though in the mathematical formulation they are part of the same overall equation.
Thus, for the \(u\) species, we have to input just the reaction component - F*(1.0 - u) - u*v*v - in the reaction input field, and then the value of \(D_u\) into the diffusion field.
For now, we can set diffusion = 100.
We can continue with the same process for \(v\); the reaction field will be u*v*v - (F + k)*v and the diffusion field, \(D_v\), will simply take a constant value of 50.
Tip
If you like, we recommend writing brief notes on the reactions, in the Notes field, for future reference and collaboration.
At this point the model structure should look like the following:
In LimbNET one can also specify the initial value for a variable, or in this case for a reactant.
Since we are modelling patterns in 2D, the initial conditions will also be in 2D.
They can be a constant value, a random value or a predefined pattern, selected from the inital value dropdown for a given reactant.
We would like to use a predefined pattern for the model's initial state - so from the initial value dropdown, select Pattern - and in this case we are going to use the square__anti-spots_multi pattern for \(u\) and square__spots_multi for \(v\).
For the time being we will not add noise to any of the species - a deterministic simulation - and thus we can leave the noise field empty.
Now that the model is fully defined, it is a good moment to save it - and now we are ready to run the simulation!
Tip
You can save your model at any time, so don't forget to save regularly.
If you want to adapt a different version of your current model you can use Save As to create a copy of the current model specification.
To run the simulation, click the simulation button in the top right corner of the editor.
The editor will prompt you to save if you did not do so already, and then you will be redirected to the simulation page.
To run the simulation click on simulate and wait for the backend response; a progress bar will inform you of the state of the simulation as results are being processed.
Once the simulation has finished, and the results have been loaded in the simulation view, we can visualise our results.
To do so, we need to select which patterns we want to associate with which colour channel in the visualisation; for example, select u in the green channel dropdown, and v in the blue channel dropdown to see the simulated pattern for u in green and the pattern for v in blue.
Then, by clicking on the play button or clicking/dragging the slider below the simulation view, we can move backwards and forwards in time through the simulated pattern and visualise the dynamic behaviour of our model:
In order to change the pattern, we need to modify the parameters k and F.
One way to do this could be to return to the editor view and change the model specification there, save it again, and re-run the simulation.
However if all we want is to make a quick parameter change, and then visualise the results, this is a bit long-winded so, instead, we can change the parameters directly without going back to the editor page.
To do so, just click on the parameters button in the simulation view and change the parameters as you like to the desired values.
Tip
Parameters that are changed while in the simulation view will not be saved in the stored model specification; changes here are temporary.
For instance, if we wanted to create a dot pattern we could use the values k = 0.062 and F = 0.028:
We re-run the simulation by closing the parameters dialog and once again clicking on simulate, and waiting for the results to arrive from the server.
As expected, these parameters produce a spotty pattern in v:
It is important to note that these changes - made here in the simulation view - do not affect the model specification itself and are not saved!
In order to modify the parameters permanently we would indeed have to return to the editor, change the values, and save the model itself.
"Pearson" plot
Ok!
We have now seen that we can simulate both stripes and dots, with the same equations but different parameters, but what if we wanted to explore the model as a whole and see which parameter combinations can generate some type of pattern?
Using the exact same set of equations, instead of treating X and Y as our spatial co-ordinate we will instead vary the value of the parameters k and F across the X and Y axes to explore how changing them affects the pattern.
We refer to this as a Pearson plot, in reference to1.
Let's begin by making a copy of the current version of our Gray-Scott model.
To do so we return to the editor view and click the save as button.
Remember to write a meaningful name (e.g., My Gray-Scott model - Pearson plot) and include a description.
For instance, we could write the following as a description to document this version of the model:
model description
This is a basic implementation of the Gray-Scott Reaction-Diffusion model. Further information about the model can be found at, e.g., https://www.mrob.com/pub/comp/xmorphia/index.html.
In this version we gradually change the parameters k and F over the X and Y axes to generate a Pearson Plot; a plot which will allow us to explore the model behaviour under various parameter combinations simultaneously.
Since we would now like our k and F values to vary over space, we will have to use 2D variables, instead of the Global Parameters that we are currently using. Moreover, we need each parameter to increase over the respective axis from some base value to a maximum, effectively forming a gradient over that axis.
To do so we can use Predefined patterns, in the mesenchyme.
So, in the mesenchyme section, by Predefined patterns, click twice on the new pattern button to add two new predefined pattern variables.
Let's firstly define our x-axis gradient variable xGrad and set its pattern to square__x_gradient_right; likewise we can define our y-axis gradient as yGrad and set its pattern to square__y_gradient_top:
To visually explore what these patterns look like, save the model, and open the simulation view by clicking on simulation.
Now, there is no need to simulate since the patterns we are interested in are predefined - simply select them in the colour channels, and you will see the X and Y gradients that we have defined!
Return to the editor view by clicking on the editor button.
We would like our k value to smoothly go from a minimal to a maximum value.
Since we may want to explore the values by changing the gradient on the fly while doing simulations, let's declare a minimum and a maximum as parameters, for example, k_min = 0.03 and k_max = 0.07.
Then, let's create a new mesenchymal variable by clicking on the new variable button in the Variables (mesenchyme) section of the editor to define our actual k value.
We can call the new variable kx and set its equation to k_min + xGrad * (k_max - k_min).
This will essentially scale our predefined X gradient, whose values go from 0 to 1, to the range k_min to k_max.
The same process can be repeated to create a gradient in the y-axis for F.
For instance, we could define two mode Global parameters, F_min = 0.0 and F_max = 0.08, and then create a second Variable in the mesenchyme called Fy with the equation F_min + yGrad * (F_max - F_min).
Finally, we need to modify our existing reaction equations for u and v in order to account for the fact that k and F are now spatial gradients with different names.
Let's change our equations to be consistent with these two new variables; change the reaction equations for u to du/dt=Fy*(1.0 - u) - u*v*v, and v to dv/dt=u*v*v - (Fy + kx)*v.
Finally, we can write the diffusion values as Global parameters for more efficient direct exploration in the simulation view.
We will have to add two more global parameters - let's call them diff_u and diff_v - and set their values to 40 and 20, respectively.
Now the mesenchyme compartment's Variables and Reactions sections should look like the following:
Make sure to save the model, and we are ready to simulate!
Return to the simulation view by clicking on the simulation button, and click on simulate.
Select the variable u under the red channel.
Next, select kx in the green channel dropdown, and Fy in the blue channel dropdown.
Finally, click on the play button or click/drag the slider below the simulation view to visualise the results.
In this way we can directly see the relationship between k and F and the resulting pattern.
Tip
Select any species - or combination of species - and then click on a triangle in the visualisation to see the value of those species at any point in space and time.
Having selected both gradients as separate colour channels, by clicking on the plot we can see the values of k and F that correspond to a particular intensity, or behaviour, of u:
As you can see in the visualisation, we have obtained a Pearson plot, but the different patterns are not as clear as they could be... One possible improvement could be to work with a finer mesh, rather than the coarse-grained model developed so far.
Fine mesh
In order to work with a finer-grained mesh, let's go back to the editor and change the bash mesh to square fixed size (fine), using the dropdown in the uppermost right corner of the editor.
You will see that once the base mesh has been changed to a different one, the patterns in the selectors (for the initial conditions, and for the predefined patterns) will disappear.
This is because any given predefined pattern exists explicitly for a given base mesh, so we need to select the equivalent patterns for our new mesh.
Tip
You can always see which predefined patterns are available for a given mesh in the Molecular Field Sequences tab of the main LimbNET interface.
First, let's fix the pattern used for the initial values of the reactions.
For the u reaction select square_fine__anti-spots_multi, and likewise for v select square_fine__spots_multi.
For the gradients, where we were using a predefined pattern variable, there is not a predefined pattern that is equivalent to the gradient we were using in the coarse mesh. But fear not - we can define this gradient ourselves as an explicit function of the spatial co-ordinates.
To do so, we will remove xGrad and yGrad from the Predefined patterns section, and instead write expressions for our x- and y-axis gradients directly in the equations for kx and Fy.
So, for kx, we will edit the equation field and replace the current equation by the expression
k_min + ((x + 1000) / 2000) * (k_max - k_min)
where x is the spatial co-ordinate in the x-axis.
Tip
Check the model specification section in the User Guide to see what other built-in variables are available.
Let's break this down a bit; we want kx to go from k_min to k_max as x goes from -1000 to 1000, the spatial extent of our square domain.
So, we will first shift x by 1000 so that it goes from 0 to 2000.
Then, we will divide by 2000 so that it goes from 0 to 1.
Finally, we will multiply by k_max - k_min so that it goes from 0 to k_max - k_min, and then we will add k_min shifting everything so that kx goes from k_min to k_max.
Likewise, for Fy we will replace the yGrad by ((y + 1000) / 2000), so we end up with our equation field for Fy containing
F_min + ((y + 1000) / 2000) * (F_max - F_min)
Don't forget to clear all the empty Predefined patterns or the model won't compile!
Save the model, and we are once again ready to simulate. The simulation will take a little longer than before, since the finer mesh requires more computations, however now you should be able to see the Pearson plot in an improved higher resolution.
-
JE Pearson. (1993) "Complex patterns in a simple system". Science, 1993; https://doi.org/10.1126/science.261.5118.189 ↩