Contact-mediated cellular communication supplements positional information to regulate spatial patterning during development
CContact-mediated cellular communication supplements positional information toregulate spatial patterning during development
Chandrashekar Kuyyamudi,
1, 2
Shakti N. Menon, and Sitabhra Sinha
1, 2 The Institute of Mathematical Sciences, CIT Campus, Taramani, Chennai 600113, India Homi Bhabha National Institute, Anushaktinagar, Mumbai 400 094, India (Dated: January 25, 2021)Development in multi-cellular organisms is marked by a high degree of spatial organization ofthe cells attaining distinct fates in the embryo. We show that receptor-ligand interaction betweencells in close physical proximity adaptively regulates the local process of selective gene expressionin the presence of a global field set up by a diffusing morphogen that provides positional cues. Thisallows information from the cellular neighborhood to be incorporated into the emergent thresholdsof morphogen concentration that dictate cell fate, consistent with recent experiments.
Spatial symmetry breaking is a fundamental prereq-uisite to morphogenesis, or the development of form, inliving organisms, such that an initially homogeneous do-main exhibits patterns in the concentrations of molecularspecies referred to as morphogens [1–4]. This can comeabout through either self-organizing reaction-diffusionprocesses [5–7] or from the anisotropy associated withthe concentration gradient of a morphogen produced bya localized source [8–10]. While in the simplest scenarioinvolving the latter mechanism, the morphogen diffusesthrough space subject to uniform linear degradation [11–15], more complex means of establishing a morphogengradient have been proposed [16–21]. Cells attain dif-ferent fates according to the positional information pro-vided by the local concentration of the morphogen vis-a-vis threshold values that emerge from the dynamics ofthe interpretation module of their genetic regulatory net-work [22–25]. However, the spatial pattern of cell fatesis not entirely determined by these local interactions asrecent experiments have highlighted the role of cell-cellcommunication in this process [26].Cells in the developing embryo are known to inter-act with other cells that are in close physical proxim-ity through contact-mediated signaling. This can occurthrough binding between membrane-bound receptors andligands on the surfaces of neighboring cells, a prominentexample being the evolutionarily conserved Notch sig-naling pathway [27]. Notch-mediated interactions, thatare believed to have arisen early in evolution, have beenshown to play a crucial role in the development of allmetazoans [27, 28]. It has been demonstrated to helpsharpen the boundaries between regions having differentcell fates in the presence of fluctuating morphogen con-centrations [29], providing an important mechanism forsystems to be robust with respect to noisy signals [30–32]. More importantly, Notch signaling is capable of self-regulation as the signaling between neighboring cells im-plements an effective feedback loop [33].In this paper we present a plausible mechanistic ba-sis for explaining how inter-cellular interactions influencecell fate determination, as indicated by recent experi- ments, e.g., on the mouse ventral spinal cord [26], by al-lowing Notch to alter the expression of genes in the mor-phogen interpretation module, which in turn control theproduction of Notch ligands. Using a three-gene interpre-tation module associated with the Sonic Hedgehog (Shh)morphogen gradient in vertebrate neural tubes [34–36],we show that specific types of Notch-mediated couplingallow the size of the domains corresponding to differentcell fates to be varied in a regulated manner. They re-tain the broad features of the reference pattern obtainedin the absence of any intercellular coupling, while avoid-ing phenotypes that do not preserve the number and thesequence of these domains [Fig. 1 (a)]. Our results sug-gest that the emergent thresholds for the morphogen con-centration that determine the localization of various cellfates are not only an outcome of the interaction betweenthe morphogen gradient manifested across an entire em-bryonic segment with the gene circuit dynamics at thecellular level, but also the intermediate-scale dynamicsof inter-cellular interactions.To investigate how the spatial patterning of cell fatesare affected by juxtacrine signaling, we consider a lineararray of cells responding to a morphogen whose concen-tration decays exponentially away from the source [11,37]. This spatial profile is reflected in the response ofthe cells in terms of the concentration of the downstreamsignaling molecules released as a result of binding of mor-phogen molecules to receptors on the cell membrane, viz., S M ( x ) = S M (0)exp( − x/λ M ), where x is the distance ofa cell from the source of the morphogen. The externalsignal concentration sensed by each cell through its recep-tors affects the expression of a set of genes that functionsas the morphogen interpretation module. We choose onethat has been proposed in the context of vertebrate neu-ral tube patterning, comprising the genes Pax6, Olig2and Nkx2.2, in the presence of a Sonic hedgehog (Shh)morphogen gradient [36]. Fig. 1 (b) shows the modulewith the regulatory motif of three patterning genes B, Wand R, that mutually repress each other, with the soleexception of W by B. The gene having the highest ex-pression level in each cell determines its fate, indicated a r X i v : . [ q - b i o . T O ] J a n FIG. 1.
Contact-mediated signaling regulates the dif-ferential expression of cell fates dictated by mor-phogen concentration profiles. (a) Schematic diagramsillustrating the French Flag problem, namely, how positionalinformation provided by spatial gradients of morphogen con-centration specify patterns of cell fates in embryonic tissue.Equally sized domains of cells exhibiting one of three differ-ent fates, viz., blue, white and red, characterize the ideal-ized situation (left), shown for the case of patterning in thevertebrate neural tube by a gradient of Sonic hedgehog mor-phogen (whose concentration profile is displayed). Under dif-ferent conditions, variations preserving the chromatic orderand number of fate boundaries of the idealized situation canarise (a: right, top row); however, other variations may violatethese (a: right, bottom row). (b) Schematic diagram of a pairof cells coupled via Notch signaling in the presence of an ex-ternal morphogen. Each cell contains a morphogen interpre-tation module comprising a regulatory circuit of fate-inducinggenes B, W and R. Notch intra-cellular domains (NICD), re-leased upon successful binding of Notch receptors to ligandsfrom the neighboring cell, affect expression of B, W and Rwith strengths θ , , , respectively. This in turn regulates theproduction of Notch ligand with strengths θ , , . (c) Spa-tial variation of the response S M to the morphogen across aone-dimensional domain comprising 30 cells. The three insetsdisplay the time evolution of gene expression levels Y (= B, Wor R, in arbitrary units) for cells that are subject to low, in-termediate and high morphogen concentrations, respectively.(d) The resulting final expression levels Y of the patterninggenes. The maximally expressed gene at each cell determinesits fate, as shown in the schematic representation of the 1Ddomain displayed at the bottom. by blue, white or red, which correspond to genes B, Wand R, respectively [Fig. 1 (c-d)]. As Pax6 is the onlygene whose expression occurs even in the absence of theShh morphogen, we consider this pre-patterning gene (B)to be expressed at very high levels initially, in contrastto the other two. The time-evolution of the expression ofthe three genes are described by: dBdt = α + ϕ N b K N (cid:0) RK (cid:1) h + (cid:0) WK (cid:1) h + ξ N b K N − k B , (1) dWdt = βS M + ϕ N b K N S M + ξ N b K N
11 + (cid:0) RK (cid:1) h − k W , (2) dRdt = γS M + ϕ N b K N S M + ξ N b K N
11 + (cid:0) BK (cid:1) h + (cid:0) WK (cid:1) h − k R , (3)where α, β, γ are the maximum growth rates and k , , are the decay rates of expression for the three genes, while K , K N and h , , , , specify the nature of the responsefunctions. The parameters ϕ , , and ξ , , are associatedwith the juxtacrine coupling of adjacent cells through thecanonical Notch signalling pathway [27, 28]. To describethe dynamics resulting from the coupling, Eqs. (1)-(3)are augmented with the time-evolution equations of theconcentrations L and N b of the Notch ligand and theNotch intra-cellular domain (NICD), respectively: dLdt = β L + φ BK + φ WK + φ RK ζ BK + ζ WK + ζ RK − Lτ L , (4) dN b dt = β N b L trans K + L trans − N b τ N b . (5)Here, the parameters β L,N b and τ L,N b correspond to themaximum growth rates and mean lifetimes for the lig-and and NICD, respectively. The binding of Notch re-ceptors of a cell to corresponding ligands of neighboringcells ( L trans ) causes the receptor’s intracellular domainto be released and translocated to the nucleus [38]. Weconsider Notch and the patterning genes to regulate theexpression of each other [see Fig. 1 (b)]. Specifically, weconsider four classes of inter-cellular interactions basedon whether NICD up or downregulates the expression ofB, W and R genes, and whether Notch ligand produc-tion is promoted or repressed by the patterning genes(mirroring the response of the ligands Jagged and Delta,respectively [39–41]). For simplicity, the ligand is as-sumed to be either activated by all the genes or inhibitedby each of them, while the genes themselves are regu-lated by NICD in a qualitatively identical manner. Thus,the four classes of inter-cellular coupling, defined by up( G + ) or downregulation ( G − ) of the patterning genes,and promotion ( L + ) or repression ( L − ) of the ligand,and specified by the parameter set ( ϕ i , ξ i , φ j , ζ j ), corre-spond to type I: G − , L − (0 , θ i , , θ j ); type II: G − , L + (0 , θ i , θ j , G + , L − ( θ i , , , θ j ); and type IV: G + , L + ( θ i , , θ j , i = 1 , , j = 4 , , ϕ i = 0, ξ i = 0, ∀ i ) yields an ide-alized flag with three chromatic regions of equal length,each corresponding to distinct cell fates (see SI). To seehow Notch signalling between adjacent cells can alter theordered pattern of cells having different fates, even whenthe morphogen gradient and the parameters of the inter-pretation module are kept unchanged, we systematicallyinvestigate the six-dimensional parameter space spannedby Θ = { θ , θ , θ , θ , θ , θ } . For each of the four typesof coupling described above, we consider 10 realizationsof the model obtained by randomly sampling Θ. Each ofthe parameters θ ,..., is sampled from the interval [1 , G + and L + ) or [0 . , .
0] (for G − and L − ). Alteringthe nature and strength of inter-cellular interactions, weobserve a diversity of resulting patterns of distinct cellfates that differ from the flag obtained in the uncoupledcase not only in terms of the lengths of the individualchromatic regions, but also in terms of their number andsequential order. To quantify the variation in the flagsobtained from the different realizations, we characterizethem by (i) the number n B of fate boundaries, which aredefined by adjacent cells having different fates, and (ii)a Hamming distance d H to the idealized flag (obtainedin the absence of coupling), determined by enumeratingthe number of cells whose fates are different in the twoflags. Depending on whether NICD up or downregulatesthe expression of the patterning genes, we obtain twoqualitatively different outcomes. While repression of B,W, R almost always results in flags having two bound-aries [Fig. 2 (a-b)], promoting their expression yields amuch wider range of n B [Fig. 2 (c-d)]. Furthermore, theflags generated for coupling types I and II are typicallycloser (in terms of d H ) to the idealized flag as comparedto types III and IV.In contrast to the parameters governing the regulationof B, W and R by NICD, those associated with modulat-ing the effect of the patterning genes on ligand productionappear to have little or no effect on the resulting flags.We use a variance-based sensitivity analysis technique toquantify the contribution of each of these parameters indetermining the cell fates [1]. We consider the final stateof each cell i comprising the domain to be represented bya discrete scalar variable F i ∈ { , , } corresponding toblue, white and red. Prior to quantifying the role playedby the parameters Θ at each cell, we quantify the vari-ance ( σ ) in the fate F i across the different realizations[Fig. 3 (a-d), upper panels]. For coupling types I andII, we note that σ is negligible throughout the array, ex-cept around the location of the two fate boundaries in theidealized flag. In contrast, σ has a finite value at all lo-cations in coupling types III and IV. The contribution ofthe different parameters θ ,..., to the observed variationin the fate of each cell is measured by the respective first- FIG. 2.
The diversity in the spatial patterns of cellfates is controlled by the nature of interactions un-derlying Notch-mediated inter-cellular coupling.
Theinter-cellular interactions can be classified into four types, de-termined by whether NICD up or downregulates the pattern-ing genes, and in turn, the genes up or downregulate ligandproduction, represented by the four motifs in the upper leftcorners in (a)-(d) [arrows representing up/downregulation areas indicated in Fig. 1 (b)]. For each type, the frequency dis-tributions of different flags, i.e., patterns representing the se-quential arrangement of distinct cell fates, are obtained byrandomly sampling θ ,..., , are shown in (a)-(d) for a 1D do-main comprising 30 cells subject to the morphogen gradientshown in Fig. 1 (c). In the absence of inter-cellular coupling,the domain is divided into three equal segments of cells hav-ing different fates [as in Fig. 1 (d)]. The flags obtained uponcoupling the cells are characterized by the number of fateboundaries n B and the difference d H with the pattern in theuncoupled system (which has equal chromatic divisions). Fortypes I, II (where NICD downregulates B, W and R) almostall flags have the same chromatic order and n B as the ideal-ized flag shown in Fig. 1 (a, left), with d H limited to very lowvalues [a and b, cf. Fig. 1 (a, right, top row)]. In contrast,the flags seen for types III, IV (where NICD upregulates B,W and R) exhibit large variation from the uncoupled case interms of both d H and n B [c and d, cf. Fig. 1 (a, right, bottomrow)]. For each type, sample flags are displayed in ascendingorder of d H along the corresponding axis. order sensitivity indices S
1, expressed as the variance of (cid:104) F i | θ j (cid:105) θ k ( (cid:54) = j ) normalized by σ (see SI). Fig. 3 (a-d, lowerpanels) show that only θ and θ contribute significantlyin all coupling types, while for coupling types III and IV, θ also plays an important role. We note that the bulk ofthe variation in F i can be explained by S FIG. 3.
Sensitivity of the flags to inter-cellular cou-pling parameters. (a-d) Dependence of the variation in cellfates on the spatial location of each cell in a 1D domain, aswell as, the differential contributions of the coupling parame-ters Θ to the variation, for the four types of Notch-mediatedinteractions. The top half of each panel shows the variance σ of the discrete variable representing the three possible fates(blue/white/red) that a cell can attain. The bottom halvesdisplay the fraction of the variance that can be accountedfor by independently varying each of the parameters (coloredaccording to the legend), as quantified by the first order sen-sitivity index S
1. When NICD downregulates the pattern-ing genes (a-b), most of the variation is localized around thetwo fate boundaries of the uncoupled case and is sensitive tochanges in θ and θ . In contrast, variation is seen across thedomain when NICD upregulates the genes (c-d), with most ofthe contribution from θ , θ and θ . (e-f) Focusing on types I,II for which chromatic order and n B of the flags are invariant,we observe that the lengths of the red and white segments ( l R and l W , respectively) are narrowly distributed around thosein the uncoupled case ( l ∗ R = 10, l ∗ W = 10). The sensitivityof the segment lengths to the parameters Θ are shown in theconcentric piecharts (outer: l R , middle: l W , inner: l B ). (g-h)The dependence of the location i b of the two fate boundaries(shown in blue and red, respectively) on the parameters θ and θ , with contour lines shown at the top. ing that the observed diversity can be largely explainedin terms of the independent actions of each parameter.As flags that do not conserve n B or the chromatic or-der of the idealized flag represent pronounced aberrationsthat are undesirable in the context of development, wefocus on coupling types I and II that are extremely un-likely to generate such flags. Indeed, the localization ofvariation in cell fates for these coupling types is consis-tent with the resulting flags typically having low d H (seeFig. 2). Moreover, almost all of them have n B = 2, whichallows the flags to be uniquely specified by the lengths ofany two out of the three chromatic regions. Fig. 3 (e-f)shows that the joint distribution of the lengths l R , l W of the regions having red and white fates, respectively, isconcentrated around that of the flag obtained in absenceof coupling (viz., l R = l W = 10 for an array of 30 cells)for both coupling types. The outer, middle and innerrings in the adjoining piecharts represent the contribu-tion of each parameter θ ,..., to the variation observedin l R , l W and l B , respectively. This is quantified by thecorresponding first-order sensitivity indices, expressed interms of the angles subtended by each of the colored seg-ments representing the different parameters. Note thatthe bulk of the observed variance in the lengths can beattributed to changes in each of the parameters, inde-pendent of the others. As θ and θ appear to be almostexclusively responsible for the observed variation in theflags, in Fig. 3 (g-h) we explicitly show how the locationsof the two boundaries i b between R, W (red) and W,B (blue) change on varying these two parameters. Forboth coupling types, increasing θ is observed to expandboth the red and blue region at the expense of the whiteregion in the middle, while increasing θ results in reduc-tion of the red region but with little impact on the W-Bboundary. Thus, the variation in the flags resulting fromdown-regulation by NICD of the genes forming the mor-phogen interpretation module can be explained by usingonly a pair of parameters controlling the repression of Wand R genes, respectively. The predicted alterations inthe resulting flag upon changing Notch expression canbe tested experimentally to validate the role of inter-cellular interactions in determining the spatial patternof cell fates outlined here.To conclude, our work reveals that juxtacrine signalingbetween cells could play a key role in adaptively regulat-ing cellular differentiation that results in morphogenesis.While the diffusing morphogen, by setting up a globalfield acts as the signal triggering the breaking of the in-trinsic symmetry, and the gene regulatory circuit form-ing the interpretation module translates the local mor-phogen concentration into the eventual cell fates, inter-cellular interactions allow the information from the en-vironment of each cell to be incorporated into the pro-cess. Apart from their utility in correcting for fluctua-tions in the signal in the presence of noise [29], such anintermediate-scale process can increase the robustness ofthe system in generating the desired flag by compensat-ing for mutations affecting the production and/or inter-pretation of the morphogen. We show this by using amodeling approach that integrates two apparently dis-parate paradigms for investigating biological pattern for-mation [31], namely, that of boundary-organized mecha-nisms involving a pre-pattern such as a morphogen con-centration gradient, and self-organized mechanisms in-volving interactions between constituents [43].We would like to thank James P. Sethna for help-ful discussions. SNM has been supported by the IMScComplex Systems Project (12th Plan), and the Center ofExcellence in Complex Systems and Data Science, bothfunded by the Department of Atomic Energy, Govern-ment of India. The simulations required for this workwere supported by IMSc High Performance Computingfacility (hpc.imsc.res.in) [Nandadevi]. [1] M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. ,851 (1993). doi:10.1103/RevModPhys.65.851[2] A. J. Koch and H. Meinhardt, Rev. Mod. Phys. , 1481(1994). doi:10.1103/RevModPhys.66.1481[3] P. Ball, The Self-made Tapestry: Pattern Formation inNature (Oxford University Press, Oxford, 1999). doi:10.1119/1.880339[4] H. T. Zhang and T. Hiiragi, Annu. Rev. Cell Dev. Biol. ,405 (2018). doi: 10.1146/annurev-cellbio-100617-062616[5] A. M. Turing, Philos. Trans. R. Soc. B , 37 (1952).doi:10.1098/rstb.1952.0012[6] H. Meinhardt, Models of Biological Pattern Formation (Academic Press, London, 1982).[7] S. Werner, T. St¨uckemann, M. B. Amigo, J. C. Rink,F. J¨ulicher, and B. M. Friedrich, Phys. Rev. Lett. ,138101 (2015). doi:10.1103/PhysRevLett.114.138101[8] L. Wolpert, J. Theor. Biol. , 1 (1969). doi:10.1016/s0022-5193(69)80016-0[9] L. Wolpert, Development , 3 (1989). doi:10.1016/j.jtbi.2010.10.034[10] J. Sharpe, Development, , dev185967 (2019). doi:10.1242/dev.185967[11] F. Crick, Nature (Lond.) , 420 (1970). doi:10.1038/225420a0[12] W. Driever and C. N¨usslein-Volhard, Cell , 95 (1988).doi:10.1016/0092-8674(88)90183-3[13] A. A. Teleman, M. Strigini, and S. M. Cohen, Cell ,559 (2001). doi:10.1016/s0092-8674(01)00377-4[14] A. D. Lander, Q. Nie, and F. Y. M. Wan, Dev. Cell ,785 (2002). doi:10.1016/s1534-5807(02)00179-x[15] A. D. Lander, Cell , 245 (2007). doi:10.1016/j.cell.2007.01.004[16] T. Bollenbach, K. Kruse, P. Pantazis, M. Gonz´alez-Gait´an, and F. J¨ulicher, Phys. Rev. Lett. , 018103(2005). doi:10.1103/physrevlett.94.018103[17] J. L. England and J. Cardy, Phys. Rev. Lett. , 078101(2005). doi:10.1103/physrevlett.94.078101[18] G. Hornung, B. Berkowitz, and N. Barkai, Phys. Rev. E , 041916 (2005). doi:10.1103/physreve.72.041916[19] D. Ben-Zvi and N. Barkai, Proc. Natl. Acad. Sci. USA , 6924 (2010). doi:10.1073/pnas.0912734107[20] S. B. Yuste, E. Abad, and K. Lindenberg, Phys. Rev. E , 061123 (2010). doi:10.1103/physreve.82.061123[21] C. B. Muratov, P. V. Gordon, and S. Y. Shvarts-man, Phys. Rev. E , 041916 (2011). doi:doi.org/10.1103/physreve.84.041916[22] J. B. Gurdon and P-Y. Bourillot, Nature (Lond.) ,797 (2001). doi:10.1038/35101500[23] H. L. Ashe and J. Briscoe, Development , 385 (2006).doi:10.1242/dev.02238[24] K. W. Rogers and A. F. Schier, Annu. Rev. Cell Dev.Biol. , 377 (2011). doi:10.1146/annurev-cellbio-092910-154148[25] S. F. Gilbert, Developmental Biology (Sinauer, Sunder-land, MA, 2013).[26] J. H. Kong, L. Yang, E. Dessaud, K. Chuang,D. M. Moore, R. Rohatgi, J. Briscoe, and B. G. Novitch,Dev. Cell , 373 (2015). doi:10.1016/j.devcel.2015.03.005[27] S. Artavanis-Tsakonas, M. D. Rand, and R. J. Lake, Sci-ence , 770 (1999). doi:10.1126/science.284.5415.770[28] R. Kopan and M. X. G. Ilagan, Cell , 216 (2009).doi:10.1016/j.cell.2009.03.045[29] D. Sprinzak, A. Lakhanpal, L. LeBon, J. Garcia-Ojalvo,and M. B. Elowitz, PLoS Comput. Biol. , e1002069(2011). doi:10.1371/journal.pcbi.1002069[30] T. Erdmann, M. Howard, and P. R. Ten Wolde,Phys. Rev. Lett. , 258101 (2009). doi:10.1103/physrevlett.103.258101[31] A. D. Lander, Cell , 955 (2011). doi:10.1016/j.cell.2011.03.009[32] A. D. Lander, Science , 923 (2013). doi:10.1126/science.1224186[33] The feedback loop arises from the Notch intra-cellulardomain (NICD) of a cell i controlling the production ofNICD in the neighboring cell j by affecting the productionof the corresponding ligand in cell i , which in turn resultsin NICD of j controlling that of i , thereby completing theloop.[34] E. Dessaud, A. P. McMahon, and J. Briscoe, Develop-ment , 2489 (2008). doi:doi.org/10.1242/dev.009324[35] E. Dessaud, V. Ribes, N. Balaskas, L. L. Yang,A. Pierani, A. Kicheva, B. G. Novitch, J. Briscoe,and N. Sasai, PLoS Biol. , e1000382 (2010). doi:10.1371/journal.pbio.1000382[36] N. Balaskas, A. Ribeiro, J. Panovska, E. Dessaud, N.Sasai, K. M. Page, J. Briscoe, and V. Ribes, Cell , 273(2012). doi:doi.org/10.1016/j.cell.2011.10.047[37] We assume that the concentration of morphogens acrossthe domain is implicitly set by a spatially uniform lin-ear synthesis-diffusion-degradation (SDD) model, whichyields an exponential gradient in the morphogen concen-tration M at its steady state. We assume that S M is aproxy for the concentration of the signaling molecule thatis triggered when morphogen molecules bind successfullywith the cell membrane receptors.[38] We assume the number of receptors in each cell to besufficiently high such that saturation is not reached.[39] H. Shimojo, T. Ohtsuka, and R. Kageyama, Front. Neu-rosci. , 78 (2011). doi:10.3389/fnins.2011.00078[40] L. J. Manderfield, F. A. High, K. A. Engleka, F. Liu, L.Li, S. Rentschler, and J. A. Epstein, Circulation , 323(2012). doi:10.1161/circulationaha.111.047159[41] M. Boareto, M. K. Jolly, M. Lu, J. N. Onuchic, C.Clementi, and E. Ben-Jacob, Proc. Natl. Acad. Sci. USA , E402 (2015). doi:10.1073/pnas.1416287112[42] A. Saltelli and I. M. Sobol, Matem. Mod. , 16 (1995).[43] J. B. A. Green and J. Sharpe, Development , 1203(2015). doi:10.1242/dev.114991 SUPPLEMENTARY INFORMATIONContact-mediated cellular communication supplements positional information toregulate spatial patterning during development
Chandrashekar Kuyyamudi, Shakti N. Menon and Sitabhra Sinha
LIST OF SUPPLEMENTARY FIGURES
1. Fig S1: The final expression levels B , W and R of the patterning genes, as well as, the concentrations of theNotch ligand, L , and NICD, N b , for each of the coupling types..2. Fig S2: Spatio-temporal evolution of the expression levels of the patterning genes B, W, and R at different cellsin a 1D array comprising 30 cells.3. Fig S3: The effect of differential expression of Notch on cell fates.4. Fig S4: “Sloppy parameter sensitivity” of the flags to inter-cellular coupling. MODEL PARAMETER VALUES parameter S M (0) λ M α β γ β L β N b K K N k k k τ L τ N b h h h h h value 100.0 0.3 4.0 6.3 5.0 5.0 5.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 6 2 5 1 1TABLE S1. The values for the model parameters used for all simulation results reported (unless specified otherwise). TEMPORAL EVOLUTION AND PARAMETER DEPENDENCE OF THE PATTERNS REPRESENTINGTHE SEQUENTIAL ARRANGEMENT OF DISTINCT CELL FATES
The fate of each cell in the linear array we consider in our simulation is determined by the expression levels ofthe three patterning genes B,W and R at t max = 100. Gene B is chosen to be the pre-patterning gene, such thatits expression level is high (= 5 arb. units) initially while the initial expression level of the genes W and R is 0.Fig. S1 shows for each of the coupling types considered in the main text, the time evolution of the expression levelsfor B, W and R with a representative set of values chosen for the coupling parameters. For coupling types I and II,B continues to be the maximally expressed gene in cells that have relatively low exposure to the morphogen, whilein cells subject to intermediate and high morphogen concentrations W and R genes (respectively) are the maximallyexpressed ones. In case of types III and IV, more than two fate boundaries can emerge and the chromatic order seenin the uncoupled system is not conserved. This implies that the coupling types III and IV give rise to a large diversityof flags, compared to types I and II.Fig. S2 shows the final expression levels of the three pattering genes B, W, R in all the cells of the linear array, alongwith the concentrations of the Notch ligand, L , and the Notch intra-cellular domain NICD, N b , for four representativesets of parameter values for each of the four coupling types.We have considered the maximal production rate of NICD, β N b to be 5 for all the results described in the maintext. To understand the implications of the over or under-expression of Notch for the results of our model, we havealso performed simulations by considering a a wide range of values of β N b . We consider only coupling types I and IIwhich conserve the number of boundaries and the chromatic order observed in the uncoupled system, ensuring thatpathological patterns of cell fates do not arise for β N b = 5. Fig. S3 shows the variation in length of the chromaticregions as a function of β N b for four different sets of coupling parameter values for each of the two coupling types. Itis evident that, depending on the values of the coupling parameters chosen, the higher expression of Notch can lead toeither increase or reduction of the red region, while the width of the blue region remains largely invariant. This resultsuggests that experimental observation of the effect of under or over-expression of Notch on the length of the regionswith different cell fates can provide us with information about the type of coupling and the strength of interactionsin different biological systems. FIG. S1.
Spatio-temporal evolution of the expression levels of the patterning genes B, W, and R at differentcells in a 1D array comprising cells. The resulting cell fates are indicated on the top surface of each panel correspondingto the coupling types I-IV. For each coupling type, a representative parameter set Θ is used for the simulation. Note that, allcells initially exhibit high expression levels for B.
FIG. S2.
The final expression levels B , W and R of the patterning genes, as well as, the concentrations of theNotch ligand, L , and NICD, N b , for each of the coupling types. The maximally expressed gene at each cell in the 1Darray determines its fate, indicated by the colors blue, white or red. For each coupling type, results obtained using differentchoices of values for the parameter set Θ are shown. FIG. S3.
The effect of differential expression of Notch on cell fates.
The variation of the spatial extent of the threechromatic regions (indicated using the colors blue, white and red) with the NICD maximal production rate, β N b , for couplingtypes I and II which preserve the chromatic order and number of boundaries seen in the uncoupled system. For all resultsshown in the main text, we have chosen β N b = 5 (indicated using a broken line). The result of under or over-expression ofNotch relative to this value is shown for different choices of values for the parameter set Θ in the case of each coupling type. PARAMETER SENSITIVITYVariance-based sensitivity analysis
To quantify the contribution of the parameters governing the regulation of the patterning genes by NICD indetermining the cell fates, we have used a variance-based sensitivity analysis often referred to as the Sobol method [1].This method can be used to investigate the effect of varying any one of the parameters, or a pair of them at a time,or any other higher-order combinations. We consider a system that yields a scalar output as a function of parameters θ j , Y = f ( θ , θ , .., θ k ) . (6)The first-order sensitivity index, which corresponds to the fraction of the total variation in Y that can be attributedto varying only θ j , keeping the other parameters fixed, is defined as S j ) = V θ i ( E θ − j ( Y | θ j )) V ( Y ) . (7)Here, − j refers to all other parameters except j . For each cell this sensitivity index is, by definition, bound withinthe range 0 ≤ S j ) ≤
1. Note that the fraction of the total variance in the output variable that can be explained bychanging one parameter at a time is given by the sum over all the first-order sensitivity indices, (cid:80) j S j ) ≤ i ( i = 1 , . . . ,
30) represented by a discrete scalar variable, F i ∈ { , , } . The contribution of eachof θ k ( k = 1 , . . . ,
6) to the observed variation in cell fates is measured by the respective first-order sensitivity indices S
1, expressed as the the variance of (cid:104) F i | θ j (cid:105) θ k ( (cid:54) = j ) normalized by σ ( i = 1 , . . . , l R , l W and l B , which correspond to the lengths of the red, white and blue regions,respectively, and which can take discrete values in the range [0 , S R , S W and S B , using quasi Monte Carlo methods.1 FIG. S4. “Sloppy parameter sensitivity” of the flags to inter-cellular coupling.
The concentric pie charts (a: type I,b: type II) show the fractional contribution of the components of the parameter set Θ, represented by colors indicated in thekey at the top, to the sensitivity of the system to parameter variation. This is expressed in terms of the spectral characteristicsof the Hessian matrix, whose eigenvalues λ obtained for a specific parameter set for each coupling type are shown at the bottom.For both coupling types, the two largest eigenvalues are comparable in magnitude while there is a large gap between these andthe subsequent eigenvalues. The bar plots (c: type I, d: type II) show the probability distribution [in logarithmic scale] of theexponential function of entropy S for the eigenvector components of the Hessian calculated for an ensemble of 10 randomlychosen parameter sets Θ. Note that, as exp( S ) provides a measure for the number of dominant components in an eigenvector,most eigenvectors are dominated by a single component. This is in agreement with the observation that almost all of theconcentric shells in the pie charts show the predominance of one color. The aggregated contribution w of each parameter to theeigenvectors (e: type I, f: type II) indicates that, consistent with the variance-based sensitivity analysis reported in the maintext, θ and θ are almost exclusively responsible for the observed variation in the flags. “Sloppy parameter sensitivity” analysis As a supplement to the variance-based sensitivity analysis, we have characterized the sensitivity of the model outputto variation of the parameters Θ = { θ , . . . , θ } using “sloppy model analysis” [2–4].This is done by varying each parameter over a relevant range and calculating a Jacobian matrix that capturesthe variation of output variables of interest. The Jacobian is then used to obtain a Hessian matrix whose spectrumindicates the sensitivity of the system to each of the parameters.The output variables that we use to characterize the sloppiness of the model are the lengths of the blue,white and red regions. These are specified using the 3-tuple ( B, W, R ), each of which can take integer valuesbetween 0 and 30, subject to the constraint that B + W + R = 30. In conventional sloppy analysis of models, aspecific set of parameter values is chosen as the reference set Θ ∗ and the results of variations from this set arethen investigated. As in our model, there is no such privileged parameter set, we have carried out the analysisusing several different Θ ∗ obtained by randomly choosing values of each of the parameters from their respective ranges.To quantify the sloppiness of our model system, we compare the lengths of the regions ( B ∗ , W ∗ , R ∗ ) obtainedusing a given set Θ ∗ with those obtained using perturbed parameter sets. Each perturbed set Θ i,j is obtained byindependently varying the value of the parameter θ j in the set Θ ∗ over the relevant range in small steps of ∆ θ , whilekeeping values of the other five parameters fixed. Thus, for a given choice of Θ ∗ and Θ i,j , the residue is R i,j = (cid:113) ( B i,j − B ∗ ) + ( W i,j − W ∗ ) + ( R i,j − R ∗ ) , (8)2where ( B i,j , W i,j , R i,j ) denotes the lengths of the regions obtained using the parameter set Θ i,j . This is then used toobtain the Jacobian matrix J as J i,j = R i +1 ,j − R i,j ∆ θ , (9)which comprises ( k −
1) rows and 6 columns, where k is the number of perturbed parameter sets considered. Thissubsequently yields the 6 × H = J T J. (10)We calculate the eigen spectrum of H and normalize the eigenvalues by the largest eigenvalue. The number of“sloppy” directions corresponds to the number of eigenvalues whose normalized magnitude (cid:28)
1. The eigenvectorcorresponding to an eigenvalue whose magnitude is very small represents an axis in the 6-dimensional parameterspace along which any variation has relatively little impact on the output. Our results indicate that most eigenvectorshave a single large component.We obtain the Hessian spectra for each of 10 different random choices of Θ ∗ (Fig. S4). Using these, we establisha hierarchy of the six parameters θ j in terms of their “weights” defined as w θ j = 1 N Θ ∗ (cid:88) Θ ∗ (cid:88) k =1 λ k v kθ j ∗ v kθ j , (11)where the first summation is over all 10 parameter sets Θ ∗ and the second is over the six components of the spectraof H . Further, λ k and v kθ j are the eigenvalue and eigenvector of the k -th component of the Hessian spectra obtainedupon varying the parameter θ j . [1] A. Saltelli and I. M. Sobol, Matem. Mod. , 16 (1995).[2] K. S. Brown and J. P. Sethna, Phys. Rev. E , 021904 (2003). doi:10.1103/physreve.68.021904[3] J. J. Waterfall,F. P. Casey, R. N. Gutenkunst,K. S. Brown, C. R. Myers, P. W. Brouwer,V. Elser and J. P. Sethna Phys.Rev. Lett. , 150601 (2006). doi:10.1103/physrevlett.97.150601[4] R. N. Gutenkunst,J. J. Waterfall,F. P. Casey,K. S. Brown,C. R. Myers and J. P. Sethna PLoS Comput. Biol.3