An Iterative Genetic and Dynamical Modelling Approach Identifies Novel Features of the Gene Regulatory Network Underlying Melanocyte Development

The mechanisms generating stably differentiated cell-types from multipotent precursors are key to understanding normal development and have implications for treatment of cancer and the therapeutic use of stem cells. Pigment cells are a major derivative of neural crest stem cells and a key model cell-type for our understanding of the genetics of cell differentiation. Several factors driving melanocyte fate specification have been identified, including the transcription factor and master regulator of melanocyte development, Mitf, and Wnt signalling and the multipotency and fate specification factor, Sox10, which drive mitf expression. While these factors together drive multipotent neural crest cells to become specified melanoblasts, the mechanisms stabilising melanocyte differentiation remain unclear. Furthermore, there is controversy over whether Sox10 has an ongoing role in melanocyte differentiation. Here we use zebrafish to explore in vivo the gene regulatory network (GRN) underlying melanocyte specification and differentiation. We use an iterative process of mathematical modelling and experimental observation to explore methodically the core melanocyte GRN we have defined. We show that Sox10 is not required for ongoing differentiation and expression is downregulated in differentiating cells, in response to Mitfa and Hdac1. Unexpectedly, we find that Sox10 represses Mitf-dependent expression of melanocyte differentiation genes. Our systems biology approach allowed us to predict two novel features of the melanocyte GRN, which we then validate experimentally. Specifically, we show that maintenance of mitfa expression is Mitfa-dependent, and identify Sox9b as providing an Mitfa-independent input to melanocyte differentiation. Our data supports our previous suggestion that Sox10 only functions transiently in regulation of mitfa and cannot be responsible for long-term maintenance of mitfa expression; indeed, Sox10 is likely to slow melanocyte differentiation in the zebrafish embryo. More generally, this novel approach to understanding melanocyte differentiation provides a basis for systematic modelling of differentiation in this and other cell-types.

Let us consider the Transcription Factor (TF) A, which regulates the two copies of gene G, G 1 and G 2 , one on each chromosome. Regulation proceeds by binding of A to the appropriate regulatory elements upstream of G 1 and/or G 2 . General TF kinetics, transcription itself , mRNA export and translation are not explicitly modelled in this approach.
Upon binding of A, the regulatory element makes a transition from the "empty" state into the "occupied" (or "regulated") state, that is , where p(A i ) (respectively p(∼ A i )) is the probability that A is bound (respectively unbound) to its binding site.
When the gene G i gets into the active state, it synthesizes the protein P with a rate g. When mRNA levels are not explicitly modelled, the rate g can be thought of as an effective rate which corresponds to the maximal expression level of each of the genes G i , and incorporates binding of general TFs, as well as mRNA dynamics [1][2][3]: Degradation of proteins is also implemented, d P being the degradation rate.
The equations representing processes (1), (2), and (3) read as follows: Let us assume that binding and unbinding of A to the regulatory element is a much faster process than protein production and degradation: Then the equilibration of r Ai will be faster than the equilibration of P , and we can solve equations (4) and (5) with the quasi steady state assumption d[r Ai ]/dt = 0: whence: Summing over i we finally obtain: where we set g P = 2g. The dynamics is driven by a function of Michaelis-Menten type, or more generally a Hill function with Hill coefficient equal to 1. This is appropriate if the TF binds to DNA as a monomer. Notice the two limits for A very small and very large: A small implies that the production rate goes to zero (no TF activating the gene), while A large reproduces the maximal expression level of the gene. Since we do not model polymerase dynamics, in particular translocation, the present approach is valid in the case of weak activators, or otherwise A small [3].

Complex formation
Let us assume that the efficiency of the regulatory TF is enhanced by complex formation. In other words, let us assume that one or more molecules of A first bind together to form a regulatory complex, and then this complex binds to DNA onto a single binding site. Process (1) is then modified as follows The symbol A n represents the complex of n monomers of A. The corresponding reaction reads: and accordingly We assume again time scale separation among the different processes involved, namely: In this way, the association reaction equilibrates first, and we can safely replace [A n ] = (k 1 /k 2 )[A] n in equation (10). The same procedure used in the previous section then gives: The production rate results in the well-known Hill function. The larger n, the more sigmoidal the shape of the Hill function.

Two binding sites: The AND gate
Let us consider now the case when two different activators A and B bind independently to two different regulatory elements. By AND gate, we mean the situation when binding of both A and B is needed in order for gene G to be activated.
The processes involved are: which at steady state give and .
The average number of activated genes is where p(A, B) is the joint probability that both A and B are bound. Therefore: The generalization of more than two TF's binding each to its specific binding site is straightforward.

Two binding sites: The OR gate
In this case gene G is activated by the binding of at least one of either A or B. Activation of G implies where [r A ] and [r B ] are as given by equations (17). This leads to the following equation for protein production: Notice that the maximal expression rate g P is reached in the limit [A] → ∞, or [B] → ∞, or both. Only for both [A] → 0 and [B] → 0, expression drops to 0.

COMPETITIVE ACTIVATION
Consider now two activators A and B competing for binding to the same binding site. In order to calculate the concentration of active genes, we proceed as follows: where we used Bayes' Theorem, and the fact that binding of A and binding of B are mutually exclusive events. Notice that the conditional probabilities p(A| ∼ B) and p(B| ∼ A) correspond to [r A ] and [r B ] respectively: Also: Combining (23) and (27), we obtain: and finally

REGULATION BY REPRESSIVE TRANSCRIPTION FACTORS
Let us consider gene G in the active state G * . When a repressor R binds to the proper binding site it downregulates expression of G. The dynamics is the same as the dynamics for activation with the only difference that the active state G * is the one with the repressor unbound. Therefore where we used The quasi steady state solution reads: where again we assumed α 0 , α 1 g P , d P . The expression level of the gene is now maximal when the concentration of repressor is zero. On the other hand, for large repressor concentrations the protein production drops to zero.

NON-COMPETITIVE REPRESSION
Let us consider now a regulatory process that involves an activator and a repressor acting on two different binding sites. This case is formally very similar to activation through two activatory proteins as discussed in section 3, the only difference being the computation of the number of active genes.
The relevant processes read: with the concentrations [r A ] and [r R ] at steady state given by: and .
Since now activation depends upon the regulatory element for the repressor to stay unbound, we have now , and therefore:

REGULATION BY COMPETITIVE ACTIVATION AND REPRESSION
Consider now an activator A and a repressor R competing for binding to the same binding site. The concentration of active genes reads: By using (28) this turns out to be and finally

MELANOCYTE DIFFERENTIATION IN ZEBRAFISH -MODEL A
Let us consider the gene regulatory network depicted in Fig. 9, Model A. The transcription factor A is provided externally so as to activate Sox10. The gene X is an unknown putative gene, which plays the same role as Phox2A in the GRN for sympathetic neurons [4]. For simplicity of notation, we define S to be the concentration of Sox10 and M the concentration of Mitfa. We make the following assumptions: • Sox10 is activated by a complex regulatory process which involves a number of transcription factors. We model this as an effective TF, which binds as a monomer to the appropriate binding site, and for which we assume the following step-like time dependency: .
Parameter A 0 fixes the maximum value of A, β fixes the rising time of the A signal, and t A defines the position in time of the step.
• Binding of Sox10 occurs in both monomer and dimer forms [5]. For simplicity we assume that only one binding site is present, and Sox10 monomers and dimers can bind in a competitive fashion.
• All other proteins produced in the circuit, acting as regulatory TF's, bind as monomers to the respective regulatory elements.
• Activators and repressors of the same gene bind non-competitively at different binding sites.
• Binding and unbinding of TFs is fast, so that a quasi steady state approximation is valid.
The regulatory processes and respective equations read: S : M : i.e. dr where by e (G) F and r we indicate the average number of empty (respectively occupied) binding sites for TF F in the promoter region of gene G.
Let us consider now the rate equation for S: The quasi steady state approximation means to set to zero the derivatives in equations (40), (41), and (42), but to keep the derivative in (43), and accordingly to get rid of the terms in the second line of eq. (43). Therefore (43) becomes: Similarly, equations for M and X read: The set of equations (44), (45), and (46) represent the full dynamics of the network depicted in Fig. 9 (Model A).

MELANOCYTE DIFFERENTIATION IN ZEBRAFISH -MODEL B
Let us consider the gene regulatory network depicted in Fig. 9 (Model B). Activation of Sox10 is followed by activation of Mitfa and Factor Y, which define a positive feedback loop. Also, repression of Mitfa upon Sox10 is now further detailed by introducing Mitfa activation of Sox10 inhibited by Hdac1, and this is effectively modelled as a competitive activation/repression process (as from section 6). On top of the assumptions made already for Model A, we also assume the OR gate (see section 2) for the activation of Mitfa by Sox10 and Factor Y. Again for simplicity of notation, we define H as the concentration of Hdac1, T as the concentration of Tyrp1, and D as the concentration of Dct. Because of the topology of the network, we report equations and numerical data only for Tyrp1 and Dct. The regulatory processes and respective equations then read: Under the quasi steady state approximation, the relative equations read: The set of equations (53), (54), (55), (56), (57), and (58) represents the full dynamics of the network depicted in Fig. 9 (Model B). At steady state the dynamics is determined by S, M , and Y only, and can be determined by setting all derivatives to zero. In particular it is interesting to consider the Sox10 mutant, which can be selected by imposing where For D M D Y δ 0 σ 0 − δ 1 σ 1 > 0 the non zero steady state is stable and the null steady state is unstable. In contrast, for D M D Y δ 0 σ 0 − δ 1 σ 1 < 0 the null steady state is stable, while the non zero one is unstable.

MELANOCYTE DIFFERENTIATION IN ZEBRAFISH -MODEL C
Model B shows three drawbacks. The first one relates to the fact that Mitfa repression of Sox10 mediated by Hdac1 does not appear sufficient to drive Sox10 expression below a detection threshold after its initial high expression phase. We cure this by imposing Hdac1 inhibition onto Factor A as well, as shown in Fig. 9 (Model C).
The second drawback relates to the unstable features of the null (M, Y ) steady state as from (65), (66) for D M D Y δ 0 σ 0 − δ 1 σ 1 > 0 in the Sox10 mutant. In this case the basal expression of either Mitfa or Factor Y would be sufficient to start the feedback loop, and get all target genes normally expressed. This is at odds with experimental observation. To avoid this circumstance we insert a threshold mechanism based on Mitfa and Factor Y concentration, which effectively mimicks a regulatory dynamics not characterized in the present model.
The third drawback is related to the behaviour of both Mitfa and Sox10 mutants. The experimentally observed existence of a transient Dct signal in both mutants implies the existence of a parallel pathway that provides only weak activation and is overwhelmed in the wild type situation, but becomes the predominant one in mutants. We assume this pathway to be mediated by Factor Z, and driven by a transient activator, Factor B.