Toward Synthetic Spatial Patterns in Engineered Cell Populations with

Subscriber access provided by ORTA DOGU TEKNIK UNIVERSITESI KUTUPHANESI

Article

Towards synthetic spatial patterns in engineered cell populations with chemotaxis Salva Duran-Nebreda, and Ricard Sole ACS Synth. Biol., Just Accepted Manuscript • DOI: 10.1021/acssynbio.5b00254 • Publication Date (Web): 23 Mar 2016 Downloaded from http://pubs.acs.org on March 28, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

ACS Synthetic Biology is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 25

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

ACS Synthetic Biology

Towards synthetic spatial patterns in engineered cell populations with chemotaxis Salva Duran-Nebreda†,‡ and Ricard V. Solé∗,†,‡,¶ †ICREA-Complex Systems Lab, Universitat Pompeu Fabra, 08003 Barcelona ‡Institute of Evolutionary Biology, UPF-CSIC, 08003 Barcelona ¶Santa Fe Institute, 1399 Hyde Park Road, 87501 Santa Fe, New Mexico, USA E-mail: [email protected]

Abstract A major force shaping form and patterns in biology is based in the presence of amplification mechanisms able to generate ordered, large-scale spatial structures out of local interactions and random initial conditions. Turing patterns are one of the best known candidates for such ordering dynamics, and their existence has been proved in both chemical and physical systems. Their relevance in biology, although strongly supported by indirect evidence, is still under discussion. Extensive modeling approaches have stemmed from Turing’s pioneering ideas, but further confirmation from experimental biology is required. An alternative possibility is to engineer cells so that self-organized patterns emerge from local communication. Here we propose a potential synthetic design based on the interaction between population density and a diffusing signal, including also directed motion in the form of chemotaxis. The feasibility of engineering such a system and its implications for developmental biology are also assessed.

keywords: Turing patterns, Pattern formation, Morphogenesis, Chemotaxis.

1

ACS Paragon Plus Environment


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

Introduction Living systems are spatially-extended, out of equilibrium systems with highly organized structures. In order to achieve such complex organization and especially during development, molecules must get distributed in ordered ways within cells and tissues. At the organ level, spatially differentiated groups involving different cell types need to exchange signals capable of modifying the molecular responses of other distant cells, ensuring proper coherent responses to the environment. 1 This requires information storage and processing, connecting different scales and thus needs the involvement of long-range dynamical processes. 2,3 A especially relevant problem in synthetic biology is the understanding of how spatially organized patterns emerge from local interactions among neighboring cells. 4 This problem was approached by the great British mathematician Alan Turing in his classical 1952 work "On the chemical basis of morphogenesis". 5,6 In that paper, Turing presented a groundbreaking approach to the problem by using a model involving two interacting molecules (the morphogens) which could also diffuse over space. Under an appropriate conditions, 7 these reaction-diffusion systems can be shown to display spatial instabilities: the combination between local amplification due to reactions together with the propagation of such amplifications due to diffusion -which plays here the role of a transmission system- is able to display periodic arrangements of morphogens concentration. Reaction-diffusion (RD) systems are usually described in terms of a system of coupled partial differential equations. In its simplest form, two morphogens M1 and M2 would interact through a nonlinear model: ∂M1 = f1 (M1 , M2 ) + D1 ∇2 M1 ∂t ∂M2 = f2 (M1 , M2 ) + D2 ∇2 M2 ∂t where the functions fn (M1 , M2 ), n = 1, 2 define the specific form of the molecular interactions taking place between M1 and M2 and that can involve different kinds of chemical 2


Page 2 of 25

Page 3 of 25

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


reactions: synthesis, degradation, catalysis and multimerization among other possibilities. 8 Turing implicitly considered morphogens as Brownian particles, which allows the use of a diffusion term Dk ∇2 Mn where Dn is the standard diffusion coefficient, weighting how fast the particles move on average. One successful way of finding Turing-generating instabilities is considering that one of the morphogens acts as an activator (hereafter A) and the other as an inhibitor (I). 9–11 Specifically, the most common interaction map includes autocatalytic activity and cross activation for the activator, and it is often assumed that the I inhibits A by either repressing its synthesis or accelerating its degradation (Figure 1A). Moreover, a widely accepted condition for instability is provided by the assumption that DI ≫ DA . In other words, that the inhibitor diffuses much faster than the activator. The RD approach has been widely successful in reproducing the patterning in different scales and systems: from skin pigmentation in animals 12,13 to organism distribution in ecological scales. 14–17 However, the identification of the molecular mechanisms underpinning some of these processes have only been recently identified, 18–20 and are still under discussion. Another perspective can be taken by asking if one or both of the morphogens can be discrete, active agents such as cells. 21–23 Moreover, cells are not constrained to a passive role, they can react to fields of chemical substances, promoting reactions or displaying active movement. While not included in Turing’s initial approach, this possibility appears as an interesting line of inquiry to consider, especially under the light of recent experimental results. 24,25 In them, it was reported that pigment cell rearrangement during development is required for proper pattern formation in zebrafish. This raises important questions regarding the role of cell movement in naturally occurring Turing patterns, whether its just a transmission system as mentioned before or can have a more relevant activity, 26,27 but it also opens up new opportunities in the synthetic pattern formation domain. This perspective has been explored in the context of natural structures by some authors, 28,29 including general models of chemotactic interactions between simple agents and chemotactic fields. The solutions to this kind of systems contain relevant biological struc-

3



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

tures and order, from travelling 29 to standing waves, 30 and are thought to underpin colony pattern formation in bacteria. 31,32 Here we will explore the possibility of creating periodic structures with finite wavelength arising from a new interaction map (figure 1B), containing also two entities: a population of cells capable of chemotactic movement (x) and a freely diffusing signaling molecule and chemoattractant (Q). Essentially, we will consider that cells constitutively synthesize Q and are drawn to it following the steepest gradient, while at the same time, the chemoattractant drives an apoptotic response in x and diffuses as a brownian particle. Given the properties assumed for each morphogen, we think that a suitable pair of candidates would be fibroblasts and TNF-α (Figure 1C). It has been established that fibroblasts are capable of chemotactic movement towards several key hormones 33,34 including our proposed signaling molecule. Additionally, these eukaryotic cells move in a way that closely matches our model, especially when compared to the swim and tumble movement of model prokaryotic organisms. 35,36 It has also been described that fibroblasts do naturally synthesize and respond to several hormones including this tumor necrosis factor. 37–39 Accordingly, since all the relevant proteins for synthesis and transport of this hormone are in place, a simple integrative or episomal vector harboring a synthetic gene that constitutively expresses TNFα would introduce the desired behavior. On the other hand, TNF-α activates antagonistic signaling pathways that can either increase survival or force apoptotic responses depending on the genetic background of cells and other signaling processes. 40–42 Fortunately, extensive efforts have been put on identifying the key effectors of both pathways and blocking the proliferative effect of TNF-α with a simple knock out is possible.

Results and discussion In order to test our hypothesis regarding the possibility of building periodic arrangements of cells using a chemotactic process, we analyzed a continuous mathematical model incor-

4


Page 4 of 25

Page 5 of 25


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 ACS Paragon Plus Environment


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

Page 6 of 25

interactions: 7 ∂x ∂ ∂Q ∂ 2x = ρ − δx xΦ(Q) − G x + Dx 2 ∂t ∂r ∂r ∂r ∂ 2Q ∂Q = αx − δq Q + Dq 2 ∂t ∂r

(1) (2)

Here both variables are space- and time-dependent, i.e. x = x(r, t) and Q = Q(r, t), where r is defined within a given spatial domain Γ. The first two terms of the right hand side of both equations include growth and decay terms. The rate ρ introduces the constant increase of cell numbers whereas −δx xΦ(Q) stands for mortality, which is affected by the concentration of Q through some functional form Φ(Q). A Hill function response is used here, i.e. we assume Φ(Q) = Q2 /(Θ2 + Q2 ). Similarly, we have an increase in Q due to the presence of cells which we assume linear (αx) and a linear decay (−δq Q) is also considered. The last terms in the right hand side of both equations accounts for random movements of both cells and molecules. These are standard diffusion terms that would lead to a uniform distribution of both components. The third term in the right hand side of the first equation incorporates a different class of phenomenon: the coupling between both populations through chemotactic interactions. Here ∂Q/∂r weights the local gradient of Q across the linear dimension. The product x∂Q/∂r measures the interaction between the local cell population and the gradient: if any of them is zero (no cells or no local gradient are present) this term is zero. The full term introduces the overall effect of these chemotactic interactions, weighted by a constant coefficient G. In order to test the presence of Turing-like instabilities, we first determine the fixed points of our model for a spatially homogeneous case, where the spatially-dependent terms can be ignored. We thus use the so called mean-field model: dx = ρ − δx xΦ(Q) dt dQ = αx − δq Q dt 6


(3) (4)

Page 7 of 25

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


First, we determine the equilibrium (fixed) points of these equations, namely the pairs S = (x∗ , Q∗ ) such that (dx/dt)P = (dQ/dt)P = 0. This leads to one positive fixed point P = (x∗ , Q∗ ) which is obtained as the intersection of two curves in the (x, Q) plane, namely a straight line x = δq Q/α and a smooth decaying function x = ρ(Θ2 /Q2 + 1)/δx . The stability of this point is obtained by analysis of this system, using the Jacobian matrix J defined as 

 J (P) = 

∂ P˙ ∂P

∂ P˙ ∂M

˙ ∂M ∂P

˙ ∂M ∂M

  

(5) P

where P˙ and M˙ indicate the time derivatives (that are replaced by the right hand sides of equations 3-4). The stability of the fixed point P is determined by the eigenvalues of this matrix, evaluated at P . Specifically, if both have a negative real part, the equilibrium state will be stable. This is a required condition while searching for spatially-induced instabilities. For our system, the matrix reads: 

∗

 −δx Φ(Q ) J (P) =  α

−δx x ∂ΦP ∂Q −δq

  

(6) P

It can be shown that the stability of the system can be determined by the signs of the trace of the Jacobian (i.e. T rJ = ∂ P˙ /∂P + ∂ M˙ /∂M ) and the determinant Det(P). Here we have the following functional forms:

T r(P) = − (δq + δx Φ(Q∗ )) < 0 ∂Φ(Q∗ ) ∗ Det(P) = δx δq Φ(Q ) + αδx x >0 ∂Q

(7) (8)

whose signs determine that the point is a stable spiral: any perturbation falls asymptotically to (x∗ , Q∗ ) with a damped oscillation (figure 2). Following the standard analysis of stability of reaction-diffusion models, we need to determine if, under the presence of diffusion dynamics, the previous fixed point becomes unstable, 7



Dq = 1.0 0.04

a

b

1.5

Dx = 0.1 0.02

Dx = 0.05

H(k 2 )

1

Q

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

Page 8 of 25

0.00

km

0.5

Dx = 0.01

-0.02

km

0 0

0.5

1

2

1.5

0

1

2

3

x

4

5

k2

Figure 2: Diffusion-induced instability. In our model without space (equations 3-4) a single fixed point P = (x∗ , Q∗ ) is present, as shown in the example (a) where we represent the isoclines (red and fray lines) whose intersection define the location of P . Several trajectories are also plotted. Here we use α = 0.075, δx = δ1 = 0.05, Θ = 1.1 and ρ = 0.01. This solution is then used in the extended model incorporating space. In (b) we display the H(k 2 ) curves as defined by equation (17), now including the additional parameters Dq = 1.0, G = 1.0 and several values of Dx . Once the asymmetry between diffusion rates is high enough, we obtain negative values of H(k 2 ) indicating wave lengths where small spatial perturbations are amplified.

leading to ordered spatial structures of a characteristic size. This linear stability starts from a spatially homogeneous state where x(r, t) = x∗ and Q(r, t) = Q∗ for all spatial points r ∈ Γ. It is easy to see that this state will be a stable solution of the spatial system (12), since all spatially-dependent derivatives are zero and thus the model collapses into the spatially-independent system (3-4). Now we consider a deviation from this homogeneous state, namely we consider a new set of values:

x(r, t) = x∗ + y(r, t)

Q(r, t) = Q∗ + q(r, t)

(9)

where y(r, t) and q(r, t) are small deviations from the stationary state, i.e. we have |y(r, t)/x∗ | ≪ 1 and |q(r, t)/Q∗ | ≪ 1. Close to this homogeneous state, the previous equations adopt a linear form, which is obtained by expanding the original equations (1-2) around the perturbed

8


Page 9 of 25

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


states (equations 9) and keeping the linear (first-order) terms. 7 This leads to: ∂y ∂Φ(Q∗ ) ∂ 2y ∂ 2x = −δx Φ(Q∗ )y − δx x∗ q − Gx∗ 2 + Dx 2 ∂t ∂Q ∂r ∂r ∂ 2q ∂q = αy − δq q + Dq 2 ∂t ∂r

(10) (11)

In order to determine if the spatial coupling (through both chemotaxis and passive diffusion) is capable of amplifying these initially small perturbations, we consider a Fourier decomposition of the initial perturbations. Specifically, we assume that the perturbation terms are a sum of Fourier modes modulated by a time-dependent exponential factor. Formally, this is written as linear superpositions of (independent) periodic waves, namely:

yk (r, t) =

X

Ak sin(kr)eΓk t

qk (r, t) =

k

X

Bk sin(kr)eΓk t

(12)

k

In order to make our analysis simpler, we consider each of these periodic functions separately, and just indicate a given Fourier component as the functions y(r, t) = A sin(kr)eΓt

q(r, t) = B sin(kr)eΓt

(13)

where k is the so called wave number, which scales as an inverse of wave length (the smaller k, the longer the associated wavelength). The exponential factor contains a parameter Γ whose sign determines the stability of a given periodic solution with wavenumber k: if Γ > 0 the perturbation will grow in time whereas Γ < 0 tell us that it will die out and no structure with a wavelength k will be present. The underlying idea under the Turing approach is that we can determine the presence of structures of a given size (characterized by k) by finding the conditions under which the time-dependent term grows in time. Using the previous periodic functions, we calculate their time- and space-dependent partial derivatives and replace the resulting functions in the linearized equations (10-11). It can

9



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

Page 10 of 25

be shown that this leads to a pair of algebraic equations: (Γδx φ1 + Dx k 2 )A + (δx φ2 − γk 2 )B = 0

(14)

−αA + (Γ + δq + Dq k 2 )B = 0

(15)

where we use φ1 = Φ(Q∗ ), φ2 = x∗ (∂Φ(Q∗ )/∂Q) and γ = Gx∗ . This equation has nontrivial solutions provided that the associated determinant is zero, i.e. (Γδx φ1 + Dx k 2 )(Γ + δq + Dq k 2 ) + α(δx φ2 − γk 2 ) = 0

(16)

The stability of the spatial perturbations is determined now by the two roots of the previous p quadratic equation for Γ. We have now Γ± = − 21 X1 ± X1 − 4H(k 2 ), where we use (δx φ1 + δq +(Dx +Dq )k 2 ) and H(k 2 ) = (δx φ1 +Dx k 2 )(δq +Dq k 2 )+α(δx φ2 −γk 2 ). These eigenvalues will

have negative real part provided that the condition H(k 2 ) > 0 and this defines the required condition for no spatial pattern. The critical boundary separating stable (homogeneous) from unstable (heterogeneous, patterned) spatial distributions is thus obtained form the marginal stability condition: H(k 2 ) = Dx Dq k 4 + (δx φ1 Dq + Dx δq − αγ)k 2 + Det(P)

(17)

2 This corresponds to a parabola in the H − k 2 plane. There is an extreme value km associated

to the condition ∂H(k 2 ) =0 ∂k 2

(18)

which is given by 2 km =

γα − δx φ1 Dq − δq Dx 2Dx Dq

10


(19)

Page 11 of 25

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


which corresponds to a minimum, since the condition ∂ 2 H(k 2 ) = 2Dx Dq > 0 ∂(k 2 )2

(20)

holds. By using this value and introducing it in the H(k 2 ) equation we obtain the value of the H function at the minimum. This result is well known within the mathematical treatment of reaction-diffusion systems displaying Turing instabilities. 7 Since the condition for stability is H > 0, those values of the wave number such that H < 0 will correspond to wave lengths that are amplified by the nonlinearities of our system. Using η = γα − δx φ1 Dq − δq Dx , this value is 2 H(km )=

η η2 + Det(P) − >0 2 2Dx Dq

(21)

The parameter space of our model can be described in terms of two main phases, corresponding to the absence H > 0 or the presence H < 0 of spatial instabilities. The critical boundary separating these two phases is given by the value kc2 such that H(kc2 ) = 0. From this value, it can be shown that the two diffusions are connected by the following condition separating uniform from patterned:

Dxc =

η2 >0 ηDq + 12 Det(P)

(22)

In figure 2b we show the parabolas associated to H(k 2 ) (a) and an example of the two phases exhibited by the model. From the previous formula we can see that the higher the diffusion rate of the chemical signal Dq the lower is the required value for cell spreading rates Dx . The previous mathematical analysis is useful because it predicts that spatial instabilities are expected to occur provided that the standard asymmetry of diffusions associated to our "activator" (x) and our "inhibitor" (Q) is present. We used a one-dimensional model with periodic boundary conditions that is simple and easy to manage. Reality however is usually

11



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

two-dimensional concerning pattern formation in natural and synthetic systems. Below we present a version of the previous model where individual cells are introduced by means of an individual-based modeling approach, whereas the chemical signal Q is again defined in terms of coupled partial differential equations. As shown below, the previous results hold and Turing-like instabilities are obtained.

Two-dimensional computational model We also created an embodied computational Netlogo model for the two-dimensional problem, incorporating the two entities previously described: cells (x) and the chemoattractant (Q). The first one is introduced as a collection of point objects, placed in a bidimensional domain Ω, each one characterized by its position and velocity. The second one is modeled as a regular square lattice with periodic boundary conditions. Following the interactions displayed in Figure 1b, we assume that cells are introduced at a fixed rate -or alternatively could be considered to differentiate from a subjacent stem cell population- with essentially randomized positions. At the same time, cells are killed -or detach from the culture surface- if the signaling molecule Q surpasses a given threshold (Θ). Moreover, cells react to the local values of signal concentration -i.e. they have access to the concentration of Q in their eight nearest neighbors- and move uphill following the chemical gradient. More precisely, our cells decide their velocity by integrating the gradient in the local neighborhood, accounting also for other cells that might be in their vicinity (in order to avoid overlaps). The first part is computed by increasing the chemotactic pull in a particular direction following a Michaelis-Menten like formulation, namely:

~vi,j =

X

(i′ ,j ′ )∈Γi,j

~u(ij→i′ j ′ )

β(Qi′ ,j ′ − Qi,j ) , (kQ + (Qi′ ,j ′ − Qi,j ))

where ~vi,j is the velocity a cell would experience if it stood in the position (i, j), Γi,j is a

12


Page 12 of 25

Page 13 of 25

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


standard Moore’s neighborhood of size one centered at (i, j) and ~u(ij→i′ j ′ ) is the unitary vector from (i, j) to each of the neighbors (i′ , j ′ ). The second part, is incorporated by checking if there is any other cell at a radius of 1.5 cell widths, then a repulsive component to the velocity is applied to both cells which keeps them separated. On the other hand, we consider that Q is locally synthesized at a constant rate by cells, decays exponentially and diffuses into a Moore’s neighborhood of size one, computed using the standard discretization: DQ ▽2 Qi,j = DQ

X

(i′ ,j ′ )∈Γ

[Qi′ ,j ′ − Qi,j ], i,j

Taking all these processes into account, a simple ODE formulation can be obtained for the finite element (i, j) in Ω, which reads as: dxi,j = B(x, Q) + ρ − δx xi,j Φ(Θ, Q) dt dQi,j = DQ ▽2 Qi,j + αx − δQ Qi,j , dt where B(x, Q) is the cellular flux associated to the chemotactic and repulsive processes, which depend on the fields of cells and signaling molecule as described before, ρ is the constant input of cells, Φ(Θ, Q) is a step function that equals 1 if Q ≥ Θ and implements the signal mediated death, α is the parameter regulating cell-dependent synthesis of Q and δQ Qi,j is a standard exponential decay for the chemotactic signal. First we wanted to assess the existence of domains in parameter space that yielded distinguishable phases in terms of the formed patterns. More precisely, under what conditions periodic structures where formed in the distributions of cells and morphogens. This periodicity has been commonly analyzed by the existence of a peak in the frequency domain of the Fourier Transform of spatial data, corresponding to dominant wavelength of the pattern. In order to examine this possibility we run the model with varying parameter sets and found

13




Page 14 of 25

Page 15 of 25

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


These are higher order properties of the distributions of morphogens, i.e. how the activated regions spatially connect and extend. Figure 4A shows how parameter variation can force the creation of the aforementioned types of structures in the distribution of cells. In particular, for low chemotactic pull and high input of cells, increasingly longer structures are formed, eventually reaching an almost completely interconnected phase. On the other hand, when cells are added less frequently to the simulation space and move slower, spots tend to be created, arranging themselves in quasi-hexagonal lattice (Figure 4B). Moreover, the rules upon which this model is based readily enable intercalation and fusion phenomena described both in other models and natural systems (arrows in Figure 4B). These refer to the apparition of new structures between distant activated regions and the collapse of two activated regions into a single spot when they are too close. An important issue regarding the pattern specification is the dominant wavelength of the spatial distribution. This typically depends on the various parameters directing the nonlinear processes of interaction among morphogens and the ratio between diffusivities. For our proposed model involving chemotactic movement, we found that three parameters are mainly responsible for defining the dominant wavelength, the diffusion coefficient (DQ ), the chemotactic pull parameter (β) and the cellular input (ρ). Figure 4C shows for a fixed value of cellular input how β and DQ affect the dominant wavelength (λ). One interesting property observed in natural systems is the formation of not only periodic structures but linear arrangements of high morphogen concentration, a well documented outcome of fish pigmentation patterning across different species. 44–47 This effect can be obtained in canonical reaction-diffusion systems with the introduction of particular boundary conditions, which create a pre-pattern and force the organization of the morphogen structures in accordance with that pre-existing information. However, another possibility arises in our model when considering that chemotactic movement is necessarily mediated by signaling pathways and that receptors can become saturated. In particular, we tested the impact of such effect by modifying the semi saturation con-

15




Page 16 of 25

Page 17 of 25




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

opmental theories that do not include self-organizing phenomena. 3,47,48 Moreover, although there are examples of artificially induced self-organization in the synthetic pattern formation domain, no clear examples of synthetic Turing-type mechanisms have been built, adding to the general impression that self-organizing systems are difficult to engineer. Here we have proposed a robust alternative, with patterns appearing over a wide range of parameter values with only alterations in the characteristic scale of the domains. Additionally, the model introduced here offers interesting predictions regarding the straightness of the activated regions, should the biochemical saturation of the signaling processes come to play a relevant role. This is in accordance to more recent trends 8 where researchers explore the effect of new interaction logics between the morphogens in order to expand the pattern forming region, thus facilitating its translation to a synthetic biology device.

Methods For the analytical treatment presented in the first section a periodic boundary condition was used to explore the formation of regular structures. In the second part, a hybrid agent-based ODE model was custom built in Netlogo (available upon request) with space discretized through a 200 × 200 grid and agent population ranging between 103 and 2 × 104 (depending on the parameter set). The boundary conditions were either periodic (Figures 3, 4c and 5c,d) or zero-flux (Figures 4a and 5b). The initial conditions were set close to the unstable solution with randomised positions for the agents, which introduces the necessary noise to promote the formation of regular patterns.

Acknowledgement We thank the members of the Complex Systems Lab for useful discussions. This work was supported by MINECO BES-2010-038940 fellowship (SDN), by the Fundación Botín through its global Universities program and the Santa Fe Institute. 18


Page 18 of 25

Page 19 of 25

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


References (1) Goodwin, B. C. (2000) The life of form. Emergent patterns of morphological transformation. C. R. Acad. Sci., Ser. III 323, 15–21. (2) Kondo, S., and Miura, T. (2010) Reaction-diffusion model as a framework for understanding biological pattern formation. Science 329, 1616–20. (3) Green, J. B. a., and Sharpe, J. (2015) Positional information and reaction-diffusion: two big ideas in developmental biology combine. Development 142, 1203–1211. (4) Basu, S., Gerchman, Y., Collins, C. H., Arnold, F. H., and Weiss, R. (2005) A synthetic multicellular system for programmed pattern formation. Nature 434, 1130–4. (5) Turing, A. (1952) The Chemical Basis of Morphogenesis. Philos. Trans. R. Soc., B 237 . (6) Ball, P. (2015) Forging patterns and making waves from biology to geology: a commentary on Turing (1952)ÔThe chemical basis of morphogenesisÕ. Philos. Trans. R. Soc., B 370, 20140218. (7) Murray, J. D. (1982) Parameter space for turing instability in reaction diffusion mechanisms: a comparison of models. J. Theor. Biol. 98, 143–163. (8) Diambra, L., Senthivel, V. R., Barcena Menendez, D., and Isalan, M. (2014) Cooperativity to increase Turing pattern space for synthetic biology. ACS Synth. Biol. (9) Koch, A., and Meinhardt, H. (1994) Biological pattern formation: from basic mechanisms to complex structures. Rev. Mod. Phys. 66, 1481. (10) Gierer, A., and Meinhardt, H. (1972) A theory of biological pattern formation. Kybernetik 12, 30–39.

19



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

(11) Meinhardt, H. (2012) Turing’s theory of morphogenesis of 1952 and the subsequent discovery of the crucial role of local self-enhancement and long-range inhibition. Interface Focus 2(4), 407–416. (12) Kondo, S., and Asai, R. (1995) A reaction-diffusion wave on the skin of the marine angelfish Pomacanthus. Nature 376, 765–768. (13) Nakamasu, A. (2009) Interactions between zebrafish pigment cells responsible for the generation of Turing patterns. Proc. Natl. Acad. Sci. U. S. A. 106, 8429–8434. (14) Kawasaki, K., Mochizuki, A., Matsushita, M., Umeda, T., and Shigesada, N. (1997) Modeling spatio-temporal patterns generated by Bacillus subtilis. J. Theor. Biol. 188, 177–185. (15) Lejeune, O., and Tlidi, M. (1999) A model for the explanation of vegetation stripes. J. Veg. Sci. 10, 201–208. (16) Thar, R., and Kühl, M. (2005) Complex pattern formation of marine gradient bacteria explained by a simple computer model. FEMS Microbiol. Lett. 246, 75–79. (17) Deng, P., de Vargas Roditi, L., van Ditmarsch, D., and Xavier, J. B. (2014) The ecological basis of morphogenesis: branching patterns in swarming colonies of bacteria. New J. Phys. 16, 015006–15006. (18) Economou, A. D., Ohazama, A., Porntaveetus, T., Sharpe, P. T., Kondo, S., Basson, M. A., Gritli-Linde, A., Cobourne, M. T., and Green, J. B. (2012) Periodic stripe formation by a Turing mechanism operating at growth zones in the mammalian palate. Nat. Genet. 44, 348–351. (19) Sheth, R., Marcon, L., Bastida, M. F., Junco, M., Quintana, L., Dahn, R., Kmita, M., Sharpe, J., and Ros, M. a. (2012) Hox genes regulate digit patterning by controlling the wavelength of a Turing-type mechanism. Science 338, 1476–80. 20


Page 20 of 25

Page 21 of 25

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


(20) Raspopovic, J., Marcon, L., Russo, L., and Sharpe, J. (2014) Digit patterning is controlled by a Bmp-Sox9-Wnt Turing network modulated by morphogen gradients. Science 345, 566–570. (21) Cocho, G., Pérez-Pascual, R., and Rius, J. (1987) Discrete Systems , Cell-Cell Interactions and Color Pattern of Animals . I . Conflicting Dynamics and Pattern Formation. J. Theor. Biol. 419–435. (22) Bonabeau, E. (1997) From classical models of morphogenesis to agent-based models of pattern formation. Artificial life 3, 191–211. (23) Adamatzky, A., Martínez, G. J., and Mora, J. C. S. T. (2006) Phenomenology of reaction–diffusion binary-state cellular automata. Int. J. Bifurcation Chaos Appl. Sci. Eng. 16, 2985–3005. (24) Frohnhöfer, H. G., Krauss, J., Maischein, H.-m., and Nüsslein-volhard, C. (2013) Iridophores and their interactions with other chromatophores are required for stripe formation in zebrafish. Development 3007, 2997–3007. (25) Mahalwar, P., Walderich, B., and Singh, A. P. (2014) Local reorganization of xanthophores fine-tunes and colors the striped pattern of zebrafish. Science 345, 1362– 1364. (26) Edelstein-Keshet, L., and Ermentrout, G. (1990) Models for contact-mediated pattern formation : cells that form parallel arrays. J. math. biol. 33–58. (27) Chen, T.-H., Guo, C., Zhao, X., Yao, Y., Boström, K. I., Wong, M. N., Tintut, Y., Demer, L. L., Ho, C.-M., and Garfinkel, A. (2012) Patterns of periodic holes created by increased cell motility. Interface focus 2, 457–64. (28) Myerscough, M., Maini, P., and Painter, K. (1998) Pattern formation in a generalized chemotactic model. Bull. Math. Biol. 60, 1–26. 21



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

(29) Wang, Z., and Hillen, T. (2007) Classical solutions and pattern formation for a volume filling chemotaxis model. Chaos 17, 037108. (30) Myerscough, M., and Murray, J. (1992) Analysis of propagating pattern in a chemotaxis system. Bull. Math. Biol. 54, 77–94. (31) Budrene, E. O., and Berg, H. C. (1991) Complex patterns formed by motile cells of Escherichia coli. Nature 349, 630. (32) Ben-Jacob, E., Cohen, I., and Levine, H. (2000) Cooperative self-organization of microorganisms. Adv. Phys. 49, 395–554. (33) Seppa, H., Grotendorst, G., Seppa, S., Schiffmann, E., and Martin, G. R. (1982) Platelet-derived Growth Factor Is Chemotactic for Fibroblasts Sources of Growth Factors Synthesis in Fibroblast Cultures. J. Cell. Biol. 92, 584–588. (34) Postlethwaite, a. E., and Seyer, J. M. (1990) Stimulation of fibroblast chemotaxis by human recombinant tumor necrosis factor alpha (TNF-alpha) and a synthetic TNFalpha 31-68 peptide. J. Exp. Med. 172, 1749–1756. (35) Kuo, S. C., and Koshland, D. (1987) Roles of cheY and cheZ gene products in controlling flagellar rotation in bacterial chemotaxis of Escherichia coli. J. Bacteriol. 169, 1307– 1314. (36) Alon, U., Surette, M. G., Barkai, N., and Leibler, S. (1999) Robustness in bacterial chemotaxis. Nature 397, 168–171. (37) Vilcek, J., Palombella, V. J., Henriksen-DeStefano, D., Swenson, C., Feinman, R., Hirai, M., and Tsujimoto, M. (1986) Fibroblast growth enhancing activity of tumor necrosis factor and its relationship to other polypeptide growth factors. J. Exp. Med. 163, 632–643.

22


Page 22 of 25

Page 23 of 25

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


(38) Battegay, E. J., Raines, E. W., Colbert, T., and Ross, R. (1995) TNF-alpha stimulation of fibroblast proliferation. Dependence on platelet-derived growth factor (PDGF) secretion and alteration of PDGF receptor expression. J. Immunol. 154, 6040–6047. (39) Mansbridge, J. N., Liu, K., Pinney, R. E., Patch, R., Ratcliffe, A., and Naughton, G. K. (1999) Growth factors secreted by fibroblasts: role in healing diabetic foot ulcers. Diabetes, Obes. Metab. 1, 265–279. (40) Miura, M., Friedlander, R. M., and Yuan, J. (1995) Tumor necrosis factor-induced apoptosis is mediated by a CrmA-sensitive cell death pathway. Proc. Natl. Acad. Sci. U. S. A. 92, 8318–8322. (41) Mohan, R. R., Mohan, R. R., Kim, W.-J., and Wilson, S. E. (2000) Modulation of TNF-a-induced apoptosis in corneal fibroblasts by transcription factor NF-kB. Invest. Ophthalmol. Vis. Sci. 41, 1327–1335. (42) Wajant, H., Pfizenmaier, K., and Scheurich, P. (2003) Tumor necrosis factor signaling. Cell Death Differ. 10, 45–65. (43) Nakamasu, A., Takahashi, G., Kanbe, A., and Kondo, S. (2009) Interactions between zebrafish pigment cells responsible for the generation of Turing patterns. Proc. Natl. Acad. Sci. U. S. A. 106, 8429–8434. (44) Shoji, H., Mochizuki, A., Iwasa, Y., Hirata, M., Watanabe, T., Hioki, S., and Kondo, S. (2003) Origin of directionality in the fish stripe pattern. Dev. Dyn. 226, 627–633. (45) Shoji, H., and Iwasa, Y. (2005) Labyrinthine versus straight-striped patterns generated by two-dimensional Turing systems. J. Theor. Biol. 237, 104–116. (46) Painter, K., Maini, P., and Othmer, H. (1999) Stripe formation in juvenile Pomacanthus explained by a generalized Turing mechanism with chemotaxis. Proc. Natl. Acad. Sci. U. S. A. 96, 5549–5554. 23



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

(47) Ball, P., and Borley, N. R. The self-made tapestry: pattern formation in nature; Oxford University Press Oxford, 1999. (48) Forgacs, G., and Newman, S. A. Biological physics of the developing embryo; Cambridge University Press, 2005.

24


Page 24 of 25

Page 25 of 25

agent-based ACS Synthetic Biology

population up gradient motion

1 2 3 4 5 6 7 8 9 10 11 12

cell death synthesis


chemoattractant field

Toward Synthetic Spatial Patterns in Engineered Cell Populations with

Recommend Documents