There are times that you discover that someone has actually named the phenomenon you have been thinking about for a significant amount
of time. This name opens the doors of an infinite set of new books, names and ideas associated with the
original one, and maybe, most importantly, you feel that you are not the only one thinking about it.
I guess, in this context, the first word I got stuck with was "asymmetry", since I was interested in particle physics in high school more than I never was. Of course I am not talking about the literal meaning of these words, and I don't even remember how old I was when I first heard the word "asymmetry". What I mean is more of a discovery of what these words actually mean in nature. As expected, then came the word "randomness", which also bothered me for a long time. These words are quite powerful, they shape the path of your curiosity, and make you discover more and more words as you dig into them. In my case, I ended up with the word "evolution" when I began college. And when I started my MSc, then came the word "diffusion". Since then, these four words are frequently present in all the books I read, as well as all the things I write. And I am incredibly happy that they made me discover my last phrase, which is hopefully my future-research to be, "mathematical biology".
So here comes the inspirational talk of this post, which contains all the words I have mentioned above. The talk is given by Philip Maini, a mathematical biologist, whom I enjoy to listen to his talks a lot. I found the chance to listen to this talk in a summer school I have recently attended, and listening to him for the second time gave me the motivation to investigate Alan Turing's paper on morphogenesis and pattern formation more carefully, and simulate it in order to understand it.
It was quite impressive to find out that Alan Turing, who is known as the father of computer by most people, had also spent quite amount of time thinking on symmetry breaking and pattern formation in nature. If you look around carefully, it is incredibly easy to see millions of examples for asymmetry. Take the branches on a tree, for instance. Their point of branching out of the trunk breaks the circular symmetry of the trunk itself. Black and white patterns of a cow is also a quite good example for spatial asymmetry. Same phenomenon is also present in the early stages of embryonic development. It all begins with a single cell, but you end up having a nervous system, muscles, bones, etc. At some point, symmetry is broken, and your initial cells begin to differentiate from each other. You can see an increadibly cool video of the embryonic development of Fruit fly Drosophila, where at the beginning, all the cells at the surface are the same, and suddenly, they begin to migrate to the places they need to be, and they begin to form a nervous system (you can see the vertical stripes developing a pattern on the surface). This whole process is called morphogenesis, and what Turing was trying to do was to model these patterns that occur through this process, mathematically.
What may be the cause of such differentiation? Recall that all the cells were daughters of the initial cell, which means that their genetic material is identically the same. In order to differentiate, these cells need an external stimulus, some kind of a chemical for instance (growth hormones can be a good example), which will trigger their process of differentiation. Turing named these chemicals morphogens, and assumed that cells begin to behave in a different way when the concentration of this morphogen is high, and vice versa. The only way for such a chemical to reach these cells, is of course, via diffusion.
Indubitably, this system is full of different morphogens interacting with each other, which also must be considered in the system model. Introducing reactions between morphogens plus diffusion gives us a reaction-diffusion system. Simplest model containing two morphogens can be described with the two following equations \begin{align} \frac{\partial u}{\partial t} &= D_{u}\nabla^2 u+f(u,v),\\ \frac{\partial v}{\partial t} &= D_{v}\nabla^2 v+g(u,v), \end{align} where \(u\) and \(v\) denote the concentrations of these two morphogens, \(D_{u}\) and \(D_{v}\) denote the diffusion coefficents of these morphogenes, and \(f(u,v)\) and \(g(u,v)\) denote the reaction kinetics between these morphogens, respectively. As you can see, without the reaction kinetics, this would be a system considering only the good old diffusion equation.
So how does this system, governed by these two equations, becomes instable, and produce these patterns? In order to answer that, we need to think of what stability is in terms of a chemical system. We need to investigate the cases where we have only reaction, only diffusion, and both, in terms of stability.
To begin with, let us assume we have a reaction-only system, which means that \(D_{u} = D_{v} = 0\), and the morphogens interact in such a way that if you mix these chemicals to be homogeneous enough, there would be a spatially uniform steady state for both \(u\) and \(v\), which is stable to perturbations [1]. This means that if I add random noise to this well-mixed system, the system will come back to its steady state. As a result, a reaction-only system can not produce any patterns, since it is stable, and comes back to a homogenous state in case of a perturbation. In other words, if the cow was grey, and I draw little black and white dots on it, it will become grey again, not due to the diffusion of these black and white dots (remember that for this case we have removed diffusion completely), but the interaction between them (you can assume such a scenario that they simply produce a grey dot, or destruct each other).
Now let us assume a diffusion-only system, where \(f(u,v)=g(u,v)=0\). This is a bit more intuitive, since heuristics tells us that if we drop a black ink into the water, after a while, the water turns grey. You can refer to this post for more information on diffusion, and how it stabilizes the system and tends to make everything homogenous as time elapses. So in this case, if the cow was grey, and I draw little black and white dots on it, it will become grey again, not due to the interaction between these dots, but due to their diffusion on the surface of the cow.
So none of these cases by themselves, reaction- or diffusion-only, gives us the big dotted cow we want. What is more counterintuitve though, is the linear combination of these systems, a reaction-diffusion system, is actually capable of doing it. If you have a stable reaction system, and you add diffusion onto it, which is also stabilizing, you will have an instable system, i.e., a real cow. In other words, stability plus stability, leads to instability.
In Turing's own words,
"It is suggested that a system of chemical substances, called morphogens, reacting together and diffusing through a tissue, is adequate to account for the main phenomena of morphogenesis. Such a system, although it may originally be quite homogeneous, may later develop a pattern or structure due to an instability of the homogeneous equilibrium, which is triggered off by random disturbances."[2]
So how does these chemicals should interact such that they can create patterns? In order to have a pattern visually, absence and presence of these two chemicals must be competitive. Turing hypothesized the existence of two molecules, an activator and an inhibitor, and in order to create patterns, these molecules must interact in a way such that,
If all the conditions above are met, in case of a small perturbation in the activator concentration (the random noise introduced to the initial homogeneous concentration), this perturbation will amplify because the activator will enhance its own production. But at the same time, since activator also promotes the production of the inhibitor, and since the inhibitor propagates faster than the activator, it will run and conquer its own areas where the activator can no longer be present. As a result, a pattern will form where there are inhibitor and activator peaks, spatially separated, which is also called "short range activation, long range inhibition".
So do the morphogen pairs really exist as Turing envisioned? They exist in chemistry, that's for sure. These simulated patterns can also be created in a petri dish. But do they exist in biology? It is still an open question.
References
[1] P. K. Maini, T. E. Woolley, R. E. Baker, E. A. Gaffney and S. Seirin Lee (2012), "Turing’s model for biological pattern formation and the robustness problem." Roy. Soc. Interface Focus "Computability and the Turing Centenary" 2(4):487-496.
[2] Alan Turing, "The chemical basis of morphogenesis".
I guess, in this context, the first word I got stuck with was "asymmetry", since I was interested in particle physics in high school more than I never was. Of course I am not talking about the literal meaning of these words, and I don't even remember how old I was when I first heard the word "asymmetry". What I mean is more of a discovery of what these words actually mean in nature. As expected, then came the word "randomness", which also bothered me for a long time. These words are quite powerful, they shape the path of your curiosity, and make you discover more and more words as you dig into them. In my case, I ended up with the word "evolution" when I began college. And when I started my MSc, then came the word "diffusion". Since then, these four words are frequently present in all the books I read, as well as all the things I write. And I am incredibly happy that they made me discover my last phrase, which is hopefully my future-research to be, "mathematical biology".
So here comes the inspirational talk of this post, which contains all the words I have mentioned above. The talk is given by Philip Maini, a mathematical biologist, whom I enjoy to listen to his talks a lot. I found the chance to listen to this talk in a summer school I have recently attended, and listening to him for the second time gave me the motivation to investigate Alan Turing's paper on morphogenesis and pattern formation more carefully, and simulate it in order to understand it.
It was quite impressive to find out that Alan Turing, who is known as the father of computer by most people, had also spent quite amount of time thinking on symmetry breaking and pattern formation in nature. If you look around carefully, it is incredibly easy to see millions of examples for asymmetry. Take the branches on a tree, for instance. Their point of branching out of the trunk breaks the circular symmetry of the trunk itself. Black and white patterns of a cow is also a quite good example for spatial asymmetry. Same phenomenon is also present in the early stages of embryonic development. It all begins with a single cell, but you end up having a nervous system, muscles, bones, etc. At some point, symmetry is broken, and your initial cells begin to differentiate from each other. You can see an increadibly cool video of the embryonic development of Fruit fly Drosophila, where at the beginning, all the cells at the surface are the same, and suddenly, they begin to migrate to the places they need to be, and they begin to form a nervous system (you can see the vertical stripes developing a pattern on the surface). This whole process is called morphogenesis, and what Turing was trying to do was to model these patterns that occur through this process, mathematically.
What may be the cause of such differentiation? Recall that all the cells were daughters of the initial cell, which means that their genetic material is identically the same. In order to differentiate, these cells need an external stimulus, some kind of a chemical for instance (growth hormones can be a good example), which will trigger their process of differentiation. Turing named these chemicals morphogens, and assumed that cells begin to behave in a different way when the concentration of this morphogen is high, and vice versa. The only way for such a chemical to reach these cells, is of course, via diffusion.
Indubitably, this system is full of different morphogens interacting with each other, which also must be considered in the system model. Introducing reactions between morphogens plus diffusion gives us a reaction-diffusion system. Simplest model containing two morphogens can be described with the two following equations \begin{align} \frac{\partial u}{\partial t} &= D_{u}\nabla^2 u+f(u,v),\\ \frac{\partial v}{\partial t} &= D_{v}\nabla^2 v+g(u,v), \end{align} where \(u\) and \(v\) denote the concentrations of these two morphogens, \(D_{u}\) and \(D_{v}\) denote the diffusion coefficents of these morphogenes, and \(f(u,v)\) and \(g(u,v)\) denote the reaction kinetics between these morphogens, respectively. As you can see, without the reaction kinetics, this would be a system considering only the good old diffusion equation.
So how does this system, governed by these two equations, becomes instable, and produce these patterns? In order to answer that, we need to think of what stability is in terms of a chemical system. We need to investigate the cases where we have only reaction, only diffusion, and both, in terms of stability.
To begin with, let us assume we have a reaction-only system, which means that \(D_{u} = D_{v} = 0\), and the morphogens interact in such a way that if you mix these chemicals to be homogeneous enough, there would be a spatially uniform steady state for both \(u\) and \(v\), which is stable to perturbations [1]. This means that if I add random noise to this well-mixed system, the system will come back to its steady state. As a result, a reaction-only system can not produce any patterns, since it is stable, and comes back to a homogenous state in case of a perturbation. In other words, if the cow was grey, and I draw little black and white dots on it, it will become grey again, not due to the diffusion of these black and white dots (remember that for this case we have removed diffusion completely), but the interaction between them (you can assume such a scenario that they simply produce a grey dot, or destruct each other).
Now let us assume a diffusion-only system, where \(f(u,v)=g(u,v)=0\). This is a bit more intuitive, since heuristics tells us that if we drop a black ink into the water, after a while, the water turns grey. You can refer to this post for more information on diffusion, and how it stabilizes the system and tends to make everything homogenous as time elapses. So in this case, if the cow was grey, and I draw little black and white dots on it, it will become grey again, not due to the interaction between these dots, but due to their diffusion on the surface of the cow.
So none of these cases by themselves, reaction- or diffusion-only, gives us the big dotted cow we want. What is more counterintuitve though, is the linear combination of these systems, a reaction-diffusion system, is actually capable of doing it. If you have a stable reaction system, and you add diffusion onto it, which is also stabilizing, you will have an instable system, i.e., a real cow. In other words, stability plus stability, leads to instability.
In Turing's own words,
"It is suggested that a system of chemical substances, called morphogens, reacting together and diffusing through a tissue, is adequate to account for the main phenomena of morphogenesis. Such a system, although it may originally be quite homogeneous, may later develop a pattern or structure due to an instability of the homogeneous equilibrium, which is triggered off by random disturbances."[2]
So how does these chemicals should interact such that they can create patterns? In order to have a pattern visually, absence and presence of these two chemicals must be competitive. Turing hypothesized the existence of two molecules, an activator and an inhibitor, and in order to create patterns, these molecules must interact in a way such that,
- Activator has an ability to enhance its own production,
- Activator promotes the production of the inhibitor,
- Inhibitor impedes the production of the activator,
- Activator propagates slower than the inhibitor, which corresponds to the difference between diffusion coefficients, \(D_{u}\) and \(D_{v}\), in the system model.
If all the conditions above are met, in case of a small perturbation in the activator concentration (the random noise introduced to the initial homogeneous concentration), this perturbation will amplify because the activator will enhance its own production. But at the same time, since activator also promotes the production of the inhibitor, and since the inhibitor propagates faster than the activator, it will run and conquer its own areas where the activator can no longer be present. As a result, a pattern will form where there are inhibitor and activator peaks, spatially separated, which is also called "short range activation, long range inhibition".
Numerical Simulation
In order to simulate this kind of a behaviour in 2D, we will discretize the system, and make use of finite element methods to find approximate solutions to the partial differential equations described in Equations (1) and (2). Since the system is discretized, we will approximate the Laplacian ( the \(\nabla^2\) operator ) using the five-point stencil finite difference method. We will also impose zero-flux boundary conditions, where none of the chemicals can diffuse out of the domain. The initial distribution of both \(u\) and \(v\) will be homogeneous, with a random noise added in order to give a perturbation for patterns to occur. Explicit equations employed in the simulations are as follows \begin{align} \frac{\partial u}{\partial t} &= 2.8 \times 10^{-4}\nabla^2 u+ 0.6 u- u^{3}-v,\\ \frac{\partial v}{\partial t} &= 5 \times 10^{-2}\nabla^2 v+1.5u-2v. \end{align} You can see the simulation results of such a system (with slightly different parameters) in Figure 1, where initially the system is homogeneous, and then random noise comes into play. You can also simulate reaction- and diffusion-only systems too, and end up with a red monochrome heatmap at the end. |
| Figure 1: Pattern formation in a reaction-diffusion system. |
So do the morphogen pairs really exist as Turing envisioned? They exist in chemistry, that's for sure. These simulated patterns can also be created in a petri dish. But do they exist in biology? It is still an open question.
Code
Data generation is done in C++ using XCode, and visualizations are done using MATLAB. I am giving both codes below, where you can generate the data, and use the MATLAB script to read it from the directory you have specified in order to create heatmaps.C++ Script (XCode)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 | // // main.cpp // turing patterns // // Created by Burcu Tepekule on 02/12/14. // Copyright (c) 2015 Burcu Tepekule. All rights reserved. // #include <iostream> #include <random> #include <fstream> #include <cmath> #include <tgmath.h> #include <cstdlib> #include <ctime> #include <stack> #include <random> #include <math.h> using namespace std; int main(int argc, const char * argv[]) { // SAY WHERE TO SAVE (Specify your own directory) string directory; directory = "/Users/burcu/Dropbox/TuringSims/"; stringstream myRoot; myRoot << directory; stringstream stream; string fileName; FILE *result; // INITIALIZE int i,k,j,t,r,m; srand((unsigned)time(NULL)); int size = 50; // Grid will be 50 x 50. Don't make it more than 100 for a reasonable runtime. double dx = 2/(double)size; int totalTime = 10; double dt = .9 *(dx*dx)/2; int timeSteps = (totalTime/dt); double **u_old; double **v_old; double **u_new; double **v_new; // Allocate memory u_old = new double*[size]; v_old = new double*[size]; u_new = new double*[size]; v_new = new double*[size]; for (i=0;i<size;i++){ u_old[i] = new double[size]; v_old[i] = new double[size]; u_new[i] = new double[size]; v_new[i] = new double[size]; } // Assign values - initial matrix is random for (i=0;i<size;i++){ for (j=0;j<size;j++){ u_old[i][j] = 3+((double)rand()/(RAND_MAX)); v_old[i][j] = 3+((double)rand()/(RAND_MAX)); } } // Save the initial matrices stream.str(string()); stream << myRoot.str() << "u_initial.dat"; fileName = stream.str(); result = fopen(fileName.c_str(),"w+"); for (i=0;i<size;i++){ for (j=0;j<size;j++){ fprintf(result, "%f\n", u_old[i][j]); }} fclose(result); stream.str(string()); stream << myRoot.str() << "v_initial.dat"; fileName = stream.str(); result = fopen(fileName.c_str(),"w+"); for (i=0;i<size;i++){ for (j=0;j<size;j++){ fprintf(result, "%f\n", v_old[i][j]); }} fclose(result); double Du = 2.8e-4; double Dv = 5e-2; for(t=0;t<timeSteps-1;t++){ // calculate 2D laplacian double **deltaU = new double*[size-1]; double **deltaV = new double*[size-1]; for(m=0;m<size-1;m++){ deltaU[m] = new double[size-1]; deltaV[m] = new double[size-1]; } for(k=1;k<size-1;k++){ for(r=1;r<size-1;r++){ deltaU[k-1][r-1] = (u_old[k-1][r]+u_old[k+1][r]+u_old[k][r-1]+u_old[k][r+1]-4*u_old[k][r])/(dx*dx); deltaV[k-1][r-1] = (v_old[k-1][r]+v_old[k+1][r]+v_old[k][r-1]+v_old[k][r+1]-4*v_old[k][r])/(dx*dx); } } for(k=1;k<size-1;k++){ for(r=1;r<size-1;r++){ // ONLY DIFFUSION SYSTEM // u_new[k][r] = u_old[k][r] + dt*(Du*deltaU[k-1][r-1]); // v_new[k][r] = v_old[k][r] + dt*(Dv*deltaV[k-1][r-1]); // ONLY REACTION SYSTEM // u_new[k][r] = u_old[k][r] + dt*(0.6*u_old[k][r]- pow(u_old[k][r],3)-v_old[k][r]); // v_new[k][r] = v_old[k][r] + dt*(1.5*u_old[k][r]-2*v_old[k][r]); // REACTION-DIFFUSION SYSTEM u_new[k][r] = u_old[k][r] + dt*(Du*deltaU[k-1][r-1] + 0.6*u_old[k][r]- pow(u_old[k][r],3)-v_old[k][r]); v_new[k][r] = v_old[k][r] + dt*(Dv*deltaV[k-1][r-1] + 1.5*u_old[k][r]-2*v_old[k][r]); } } for(k=0;k<size;k++){ // ZERO FLUX BOUNDARY CONDITION u_new[0][k] = u_old[1][k]; u_new[k][0] = u_old[k][1]; v_new[0][k] = v_old[1][k]; v_new[k][0] = v_old[k][1]; u_new[k][size-1] = u_old[k][size-2]; u_new[size-1][k] = u_old[size-2][k]; v_new[k][size-1] = v_old[k][size-2]; v_new[size-1][k] = v_old[size-2][k]; } delete deltaU; delete deltaV; for(k=0;k<size;k++){ for(j=0;j<size;j++){ u_old[k][j]=u_new[k][j]; v_old[k][j]=v_new[k][j]; } } }; // Save the last matrices stream.str(string()); stream << myRoot.str() << "u.dat"; fileName = stream.str(); result = fopen(fileName.c_str(),"w+"); for (i=0;i<size;i++){ for (j=0;j<size;j++){ fprintf(result, "%f\n", u_new[i][j]); }} fclose(result); stream.str(string()); stream << myRoot.str() << "v.dat"; fileName = stream.str(); result = fopen(fileName.c_str(),"w+"); for (i=0;i<size;i++){ for (j=0;j<size;j++){ fprintf(result, "%f\n", v_new[i][j]); }} fclose(result); return 0; } |
MATLAB Script
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | clear all;clc;close all; %Specify the same directory with the C++ code. directory = '/Users/burcu/Dropbox/TuringSims/'; u = readFile('u',directory); v = readFile('v',directory); u_initial = readFile('u_initial',directory); v_initial = readFile('v_initial',directory); size = sqrt(length(v)); u=reshape(u,size,size); v=reshape(v,size,size); u_initial=reshape(u_initial,size,size); v_initial=reshape(v_initial,size,size); HeatMap(v,'Colormap','hsv'); HeatMap(u,'Colormap','hsv'); |
1 2 3 4 5 6 7 | function [ variable_out ] = readFile(varname,directory) longname = [directory varname '.dat']; fid = fopen(longname,'rt'); variable = textscan(fid, '%f', 'HeaderLines',0); variable_out = variable{1}; fclose(fid); end |
References
[1] P. K. Maini, T. E. Woolley, R. E. Baker, E. A. Gaffney and S. Seirin Lee (2012), "Turing’s model for biological pattern formation and the robustness problem." Roy. Soc. Interface Focus "Computability and the Turing Centenary" 2(4):487-496.
[2] Alan Turing, "The chemical basis of morphogenesis".
