A Multigenerational Turing Model Reproduces Transgressive Petal Spot Phenotypes in Hybrid Mimulus

The origin of phenotypic novelty is a perennial question of genetics and evolution. To date, few studies of biological pattern formation specifically address multi-generational aspects of inheritance and phenotypic novelty. For quantitative traits influenced by many segregating alleles, offspring phenotypes are often intermediate to parental values. In other cases, offspring phenotypes can be transgressive to parental values. For example, in the model organism Mimulus (monkeyflower), the offspring of parents with solid-colored petals exhibit novel spotted petal phenotypes. These patterns are controlled by an activator-inhibitor gene regulatory network with a small number of loci. Here we develop and analyze a model of hybridization and pattern formation that accounts for the inheritance of a diploid gene regulatory network composed of either homozygous or heterozygous alleles. We find that the resulting model of multi-generational Turing-type pattern formation can reproduce transgressive petal phenotypes similar to those observed in Mimulus. The model gives insight into how non-patterned parent phenotypes can yield phenotypically transgressive, patterned offspring, aiding in the development of empirically testable hypotheses.


Introduction
The origin of phenotypic novelty is a persistent question in evolutionary genetics.In answering this question, evolutionary biologists often focus on intergenerational B Gregory D. Conradi Smith greg@wm.educhanges over long time scales.However, it is well known that a dramatic amount of phenotypic variability can arise within a single generation.In plant biology, for example, breeding divergent varieties produces offspring whose yield is significantly greater than the yield of either parent, a phenomenon known as heterosis, (i.e., hybrid vigor).In heterosis, the interaction between alleles from different parents somehow produces an emergent phenotype.This paper is a theoretical investigation of this source of phenotypic novelty.
Previous work has demonstrated that crosses between different Mimulus species with solid-colored petals can produce F1 offspring with spotted petals (P A ×P B → F1).In particular, crosses of M. luteus var.variegatus (red) with M. cupreus (yellow or orange) produce F1 hybrids that are yellow with red spots (Fig. 1A, bottom two rows).The distribution of phenotypes in an F2 population generated from the selfing of F1 hybrid plants (F1 × F1 → F2) yields additional novelty (Cooley and Willis 2009).We observe a complex distribution of phenotypes with petals having a range of spot number, size, location, intensity, and color (Fig. 1B).We seek to investigate the explanation for phenotypes changing so dramatically over just two generations.A working hypothesis must begin with the current understanding of the gene regulatory network controlling anthocyanin production, the red-purple pigment found in Mimulus petals.
In monkeyflowers, the biochemical pathway that produces anthocyanin pigment is controlled by MYBs, a well-known transcription factor (TF) family (see Stracke et al. 2001).These include two R2R3-MYB TFs, NEGAN/MYB5 in the nectar guide and PELAN in the petal lobes, and one R3-MYB TF, red tongue (RTO) (Ding et al. 2020).Two accessory proteins, a basic helix-loop-helix and a WD40 (MlANbHLH1 and MlWD40a, respectively), form a complex with the activator MYBs (called the MBW complex).The MBW complex has been widely hypothesized to operate as an activator within a reaction-diffusion mechanism, while the R3-MYB operates as an inhibitor (Bouyer et al. 2008;Ding et al. 2020;Ishida et al. 2008;Larkin et al. 1996;Pesch and Hülskamp 2004).Research on various MBW complexes has shown that the R3-MYBs (the inhibitors) move intercellularly, expanding away from the cells in which they are synthesized; R2R3-MYBs (the activators) do not (Albert et al. 2014;Ding et al. 2020;Kurata et al. 2005).These results are consistent with the requirements for long-range inhibitor and short-range activation in Turing-type pattern formation (Turing 1990;Meinhardt 2012;Kondo and Miura 2010;Meinhardt 2009).
To confirm that a reaction-diffusion-mediated pattern formation mechanism is plausible in monkeyflowers, researchers performed genetic manipulation on a species whose wild type always produces spotted nectar guides.These results show that elimination of activator via NEGAN RNAi knockouts exhibit nectar guides lacking anthocyanin.Inhibitor was suppressed in homozygous and heterozygous RTO mutants exhibit, showing very high and intermediate amounts of anthocyanin compared to the wild type, respectively.Finally, activator was eliminated throughout the flower with WD40a knockout mutants which resulted in flowers devoid of anthocyanin (Ding et al. 2020).
Based on these empirical results, Ding et al. (2020) hypothesized that a Turing-type reaction-diffusion mechanism mediates Mimulus petal patterning.They demonstrated that a simple two-variable (activator and inhibitor) reaction-diffusion system can mimic experimentally observed phenotypes (e.g., variation of spot size, elimination of (1) These PDEs account for the diffusion of activator and inhibitor and include nonlinear reaction terms for three of the four activator-inhibitor interactions typically found in Turing systems (no auto-inhibition) (Ding et al. 2020).
The Ding et al. model shows that activator-inhibitor reaction-diffusion equations can mimic the wild-type spotted nectar guide phenotype and experimental perturbations.This work is similar in spirit, but our focus is the appearance of novel (i.e., transgressive) patterned phenotypes in F1 hybrids and the variety of patterns found in the F2 generation.Section 2 shows how to model inheritance within the reactiondiffusion framework, beginning with reaction terms that instantiate a hypothesized mechanism of transcription factor binding and regulation.This enables the derivation of a reaction-diffusion system that accounts for diploidy-in particular, the effect of heterozygosity-on gene regulatory network function and pattern formation.Because genetic inheritance occurs by propagating parameters across simulated generations, we refer to our framework for modeling emerging petal patterns in hybrid Mimulus as a multi-generational Turing model.(haploid).The horizontal lines represent two unlinked genetic loci that code for transcription factors: an activator u and an inhibitor v. Filled circles represent transcription factor binding sites with affinities κ uu , κ uv , κ vu , and κ vv .The parameter α u is the baseline production rate of activator u.The parameter γ u is the increase in production rate that is regulated by the binding of the activator and the inhibitor (magenta curves).Similarly, α v and γ v are the baseline and regulatable production rate of inhibitor v (see Eq. 3) (Color figure online)

Model Formulation
The biochemical pathway that produces anthocyanin pigment in monkeyflowers involves activator and inhibitor transcription factor complexes.Our model formulation begins with a system of reaction-diffusion equations for an activator (u) and an inhibitor (v), where D u and D v are diffusion coefficients.The reaction terms, f (u, v) and g (u, v), are As illustrated in Fig. 2, these expressions represent an idealized gene regulatory network with the following properties.The activator u is produced at the baseline rate, α u , and degrades at a rate proportional to u (with a first-order rate constant β u ).In addition, the production rate of u may be increased (by as much as γ u ) through a Hill function that represents the competitive binding of activator u and inhibitor v to a pair of TF binding sites.For simplicity, these binding sites are assumed to be identical and independent.The inhibitor is modeled similarly, using the parameters α v , β v , and γ v .
The Hill functions occurring in Eq. 3 involve four equilibrium association constants: κ uu , κ uv , κ vu , and κ vv .In each case, the first subscript denotes the TF, u or v, whose production rate is being regulated; the second subscript denotes the TF that is binding.For example, κ uv is the association constant for the inhibitor (v) binding to a site that regulates the activator production rate (see Fig. 2).
Fig. 3 Schematic diagram of the diploid activator-inhibitor gene regulatory network.In essence, the diploid regulatory network is obtained by duplicating Fig. 2.There are four gene products, i.e., one activator and one inhibitor allele from each parent.The parental haplotype is distinguished by the presence or absence of a caret (u and v from parent 1, û and v from parent 2).Filled circles represent transcription factor binding sites with affinities κ u * , κ û * , κ v * , and κ v * (where * is a placeholder for u, v, û and v).Thick and thin curves indicate an interaction between binding sites and transcription factors that are inherited from the same and opposite haplotypes, respectively (Color figure online)

Accounting for Diploidy
The reaction-diffusion system presented above (Eq.2) represents two genetic loci (one for activator, one for inhibitor) known to exist on different chromosomes.To account for multi-generational aspects of Mimulus genetics, it is necessary to duplicate the model to represent two copies of each locus (one for each haploid genome).The resulting four-variable reaction-diffusion system takes the form where u and û represent the (possibly distinct) gene products associated with the activator locus, and similarly for v and v (see Fig. 3).The reaction terms appropriate for diploid Mimulus are an elaboration of those in Eq. 3, (5) The parameters α * , β * , and γ * are the production rate, degradation rate, and maximum regulation rate of each of the four gene products (the subscript u, û, v, v).The four Hill functions in Eq. 5 include sixteen equilibrium association constants κ * * .As in Eq. 3, the first subscript of these binding affinities denotes the TF (u, û, v, or v) whose production rate is being regulated.The second subscript denotes which TF is binding (also u, û, v, or v, for a total of 16 parameters).For simplicity, our simulations assume that the rates of diffusion for both activator and inhibitor are independent of allele type (D û = D u , D v = D v ).Consistent with Turing-type pattern formation, the diffusion coefficient for the inhibitor is assumed to be greater than activator (D v > D u ).

Parameter Assignment: Homozygous Parents and Doubly Heterozygous F1 Hybrid
Consider heterozygous F1 offspring from the cross of two parents with distinct alleles for both activator and inhibitor (P A × P B → F1).Table 1 shows the assignment of rate constants (α * u , β * u , . .., γ * v ) for the three diploid reaction-diffusion systems (Eq.4) that model this cross.Superscripts indicate allele type, e.g., α A u is the baseline production rate for the activator derived from parent A. The allele type chosen for each parameter follows Mendelian inheritance of the diploid genotypes for the activator (U A A , U AB , and U B B ) and inhibitor (V A A , V AB , and V B B ). Table 2 shows the inheritance of the association rate constants (κ * * * * ).Here, the interaction between the activator and inhibitor, each with two potentially different allele types, leads to the 16 possible values.The first subscript-superscript pair denotes the allele type of the TF (κ ) whose production rate is being regulated.The second subscript- ).The F1 hybrid model uses all 16 binding constants.Each of the homozygous parents utilizes 4 binding constants (parent A highlighted pink, parent B highlighted cyan).That is, eight binding constants are only relevant for heterozygous Mimulus.Note that although the genotype U B A is not distinguishable from U AB , and similarly for V B A and V AB , the binding constants κ B A * * are distinct from κ AB * * .For example, the parameter κ B A vv refers to the binding of the allele type A inhibitory transcription factor to the binding site associated with the regulation of type B inhibitor.The binding constant κ AB vv has a different interpretation and need not have the same value as κ B A vv .Using both Tables 1 and 2 and Eq. 4 we derive the reaction-diffusion system for the heterozygous F1 offspring of the P A × P B → F1 cross: These equations can be compared and contrasted with the equations for parent A: Because parent A is homozygous, the equations for u and û are identical, as are the equations for v and v.It is a simple matter to derive an equivalent reactiondiffusion system with one equation for each distinct gene product.Defining the total concentration of activator and inhibitor as ū = u + û and v = v + v, we arrive at the equivalent two-variable system, As expected, the zeroth-order rate constants (α and γ ) are scaled by a factor of 2, while the first-order rate constants (β) are not.That is, a homozygous diploid model is equivalent to an appropriately scaled haploid model (compare Eq. 8 to Eqs. 2 and 3).The reaction-diffusion system for parent B (also homozygous) is identical to Eq. 7 with the replacement of B for A.

Inheritance of Parameters in a Simulated Population of F2 Hybrids
Hybrids of M. cupreus and M. l. variegatus show a striking distribution of phenotypes in the F2 generation (recall Fig. 1).In these experiments, the F2 hybrids are produced by selfing an F1 flower and growing the progeny (F1 × F1 → F2).In our model, the distribution of genotypes in the simulated F2 generation is achieved by assigning parameters in Eq. 4 according to the Mendelian logic of the previous section.
Because the F1 hybrid is heterozygous for both activator and inhibitor (U AB V AB ), the F2 hybrid population exhibits 9 distinct genotypes, including the doubly homozy-

F2 Hybrid
The highlighted parameters differ from the doubly heterozygous F1 hybrid U AB V AB (cf.Tables 1 and 2) gous parental genotypes U A A V A A and U B B V B B , the doubly heterozygous F1 hybrid genotype U AB V AB , and 6 genotypes that are unique to the F2 generation (

and U B B V AB
).The parental and F1 cases were described in the previous section.For an example that is unique to the F2 hybrids, consider the U AB V B B genotype, which is homozygous for activator and heterozygous for inhibitor.Simulations for this genotype use the rate and binding constants in Table 3.Because this genotype is homozygous for the inhibitor (V B B ) but not the activator (U AB ), the corresponding contracted reaction-diffusion system has three equations (u, û and v = v + v): The reaction-diffusion equations for the remaining genotypes are derived similarly (see Appendix A).
Note that the parameters used to simulate the doubly heterozygous F1 hybrid (U AB V AB ) effectively define a given multi-generational simulation.Each of the other 8 genotypes (two parents, 6 unique F2 genotypes) use a subset of the F1 parameters.

The doubly homozygous genotypes (U
) use four of these, while genotypes which are heterozygous at one locus and homozygous at the other use nine (U A A V AB , U B B V AB , U AB V A A , and U AB V B B ).

Parameters for Transgressive Hybrid Phenotypes
Recall that a Turing instability for a two-variable activator-inhibitor system (e.g., Eqs. 2 and 3) requires linear stability of the reaction terms, f (u, v) and g(u, v).That is, the Jacobian defined by must be stable (both eigenvalues of A(u ss , v ss ) must have negative real part) when evaluated at the steady state satisfying 0 = f (u ss , v ss ) and 0 = g(u ss , v ss ).We constrained our simulations to unique steady-state solutions, though the model can produce systems with three steady states.Additionally, the matrix that arises from linearizing the full reaction-diffusion system, must be unstable (i.e., at least one eigenvalue has positive real part) for some spatial frequency k.It can be shown that this requires d = D v /D u > 1 (see Murray 2001, chap.2).
For our model of patterning in diploid Mimulus (Eq.4), the 4 × 4 Jacobian matrix is obtained by linearizing the reaction terms (Eq.5).For compactness, these can be written as where The Jacobian takes the form with diagonal blocks and off-diagonal blocks In these expressions, U u , U v , . . ., V v are evaluated at the steady state (u ss , v ss , ûss , vss ) solving 0 = f (u ss , v ss , ûss , vss ), and similarly for g, f , and ĝ (Eq.11).A valid multigenerational parameter set (i.e., one resulting in a transgressive F1 phenotype) must yield three stable Jacobian matrices: J A , J B , and J F1 .Additionally, the matrix that results from linearizing the diploid reaction-diffusion system, must be unstable for the F1 hybrid, yet stable for both parents.When Eq. 15 is rewritten in terms of the relative diffusion coefficient, d = D v /D u , there exists a critical value of d, denoted d , above which the matrix J will be unstable for some spatial frequency k (a Turing bifurcation).Assume without loss of generality that d A ≤ d B .In that case, a valid multi-generational parameter set has critical diffusion ratios that satisfy , we obtain a multi-generational parameter set that results in Turing-stable parents and a Turing-unstable F1 hybrid (see Appendix B for further discussion).

Results
An intriguing observation of the experimental Mimulus system is that inbred (i.e., homozygous) parent flowers with unpatterned petal phenotypes produce reliable patterned phenotypes in F1 hybrids (recall Sect.1).In particular, crosses between solid colored yellow morph M. cupreus and M. l. variegatus lead to spotted F1 hybrids (Fig. 1).To see if this phenomenon could be reproduced by our multi-generational model (Sect.2), we sought parameter sets that yield this transgressive phenotype.We require parameters for which Parents A and B are Turing stable, while the F1 hybrid is Turing unstable, as discussed above.

Unpatterned Parents Can Produce Patterned Hybrid Offspring
Figure 4A shows a representative simulation in which parents A and B are unpatterned while the F1 hybrid is patterned.These simulations are performed using an alternating-direction implicit Crank-Nicholson-like numerical scheme with no flux boundary conditions.As a proxy for the concentration of anthocyanin pigment, the steady-state value of total activator concentration (u + û) is shown in yellow-to-red pseudo color.The total inhibitor concentration (v + v) is not shown because it does not influence the visual appearance of Mimulus petals, but in these simulations it is in phase with the activator.
Figure 4B shows the simulated pheonotypes of the F2 generation.The 3 × 3 grid organizes the results by genotype.The parent and hybrid F1 cases are recapitulated on the diagonal.As in the Mimulus experimental system, the simulated F2 hybrid population exhibits a wide variety of phenotypes.For this parameter set, 5 of the 6 phenotypes unique to the F2 generation are patterned (Turing unstable) while 1 is solid (Turing stable).The patterned phenotypes consist of spots with different intensity, size, and wavelength (e.g., spots in U B B V A A and U B B V AB are similar in size, but those in U B B V AB are lighter and closer together).
Figure 5 shows four more simulated F2 hybrid populations derived from unpatterned parents whose cross yields a patterned F1 hybrid (each displayed in a 3 × 3 grid akin to that in Fig. 4B).In three of the four examples (panels B, C, and D), the unique F2 phenotypes include one or more individuals with a pattern distinct from the F1 hybrid.The number of patterned individuals and variety of patterns among the F2 hybrids can be low (panel A) or high (panel D) depending on model parameters.3) produces unpatterned phenotypes when using homozygous parent parameters.Using heterozygous parameters corresponding to the cross between unpatterned parents (P A × P B → F1), the model produces patterned phenotypes in F1 offspring.Yellow-to-red pseudo color represents the total activator concentration (i.e., low to high, u + û).B Simulated breeding of F1 hybrids (F1 × F1 → F2) yields various phenotypes, each corresponding to one of nine F 2 genotypes (U A A V A A , U AB V AB , etc.).These diploid genotypes result from all possible combinations of parent haplotypes (denoted by subscripts A and B) that are composed of two alleles each for activator (U ) and inhibitor (V ) loci.Note: the F2 individuals with genotypes U A A V A A , U AB V AB , and U B B V B B are identical to parent A, the F 1 hybrid, and parent B, respectively.Parameters are as in

Parents with Identical Unpatterned Phenotypes Can Produce Patterned Offspring
The multi-generational Turing model produces a transgressive patterned F1 phenotype and a broad distribution of F2 phenotypes (Sect.3.1) consistent with the experimental Mimulus system.Note that the cross between M. l. variegatus and M. cupreus (Fig. 1) involves unpatterned parents that are phenotypically distinct (M.l. variegatus is pink, M. cupreus is yellow).This is recapitulated in the simulated cross of Fig. 4, in which Parent A has a high concentration of activator (as does M. l. variegatus) while Parent B has a low activator concentration (like yellow M. cupreus).
Figure 5A shows that our model formulation can produce transgressive F1 phenotypes even when the unpatterned parents are phenotypically similar (one yellow and one light orange).This observation highlights the transgressive nature of the patterned F1 and F2 phenotypes produced by the model.That is, the maxima of spatially inhomogeneous activator concentration in patterned F1 and F2 individuals can exceed the spatially homogenous activator concentrations of both unpatterned parents (e.g., genotype U AB V AB in Fig. 6A).In some cases, unpatterned F2 individuals exhibit activator concentration that is more extreme than either parent (e.g., genotype U B B V A A in Fig. 6C).

A Possible Mechanism for the Emergence of Transgressive F1 Hybrid Phenotypes
Figure 5A shows a transgressive F1 hybrid (red spots on yellow background) from a cross between phenotypically similar unpatterned parents (one yellow and one light orange).Interestingly, Fig. 6 shows that phenotypically identical unpatterned parents (U A A V A A and U B B V B B ) can also yield patterned F1 hybrid offspring (U AB V AB ) and novel patterned F2 phenotypes.We will discuss this limit case in detail because it suggests one possible mechanism for the emergence of transgressive phenotypes in hybrid Mimulus.
The parents in Fig. 6 are phenotypically identical (unpatterned yellow) because, in the parameter set used, the rate constants do not depend on the allele type (α , and similarly for v).The eight binding constants used in the parent simulations (both doubly homozygous genotypes) also do not depend on allele (κ For parameters with this symmetry, the reaction-diffusion systems representing parent A and B are the same; consequently, the parent simulations yield identical phenotypes (unpatterned yellow, see U A A V A A and U B B V B B in Fig. 6).The simulation of the doubly heterozygous F1 hybrid involves these parental parameters as well as eight additional binding constants chosen to satisfy κ AB uu = κ B A uu < κ A A uu = κ B B uu (and similarly for κ uv , κ vu , and κ vv ).This choice is consistent with small structural differences in the transcription factors from alleles A and B leading to less effective binding to the regulatory side associated with the alternate allele type.
The central panel of Fig. 6 (U AB V AB ) shows that a patterned F1 hybrid is possible under the assumptions made in the previous paragraph.In this simulation, the F2 hybrids include five distinct phenotypes, four of which are patterned.Each of the three patterned phenotypes that are unique to the F2 population is produced by two equivalent genotypes (U Note that these genotype pairs result from swapping one or both of the homozygous loci (activator Figure 6 demonstrates that transgressive phenotypes in doubly heterozygous F1 hybrids can emerge when the TFs derived from one parent are less effective at regulating TFs derived from the other parent (with different allele type).For example, the inequality κ B A vu < κ A A vu means that the activator of allele type A has lower affinity for the binding site of inhibitor B than A. We will refer to regulation between identical and distinct allele types as cis and trans, respectively.In cis regulation, both the tran-Fig.6 Transgressive patterns arising from phenotypically identical parents.A The diploid reaction-diffusion model (Eq.4) can produce transgressive patterns in the F1 generation even when the homozygous parents have identical unpatterned phenotypes.In the simulated cross (P A × P B → F 1 ), the heterozygous parameter set is identical to the parents except for trans-interaction binding affinities (e.g., uu , and similarly for κ uv , κ vu , and κ vv ).B The F 2 hybrids include five distinct phenotypes, four of which are patterned (Color figure online) scription factor and its regulatory binding site are derived from the same allele type.In trans regulation, the transcription factor regulates the production of the alternate allele's gene product.To explore this mechanism further, observe that the parameter symmetries of this limit case allow the reaction-diffusion equations for the F1 hybrid to be contracted.Writing the total concentration of activator and inhibitor as ū = u + û and v = v + v, Eq. 6 is equivalent to The assumption that trans binding is less effective than cis binding implies that 0 ≤ m * * < 1. Setting m * * = 1 gives equations for the phenotypically identical parents.
Figure 7 summarizes a parameter study for which the only distinction between parents and F1 hybrid parameters is a decreased trans binding affinity for one of the four regulatory pathways (κ uu , κ uv , κ vu , or κ vv ). Figure 7A shows that reduction of the trans efficacy of the activator regulating inhibitor (decreasing m vu ) decreases the critical diffusion ratio d .This suggests a hypothesis for the emergence of transgressive phenotypes in an F1 hybrid derived from phenotypically identical parents: a Turing bifurcation occurs as a result of a decrease in the inhibitor production rate, κ vu (1+m vu ).In the parents, κ vu (1 + m vu ) = 2κ vu , but in the hybrid κ vu ≤ κ vu (1 + m vu ) ≤ 2κ vu .Figure 7C shows that a similar result is obtained with a reduction of the trans efficacy of the inhibitor binding the activatory regulatory site (decreasing m uv ).
Figure 7 demonstrates that reduction in trans efficacy of inhibitor regulating activator (m uv ) or activator regulating inhibitor (m vu ) can lead to instability in the F1 hybrid, but this is not the case for m uu or m vv .That is, decreasing m uv (or m vu ) decreases d , but decreasing m uu (or m vv ) increases d .These results depend on the parameters chosen for the phenotypically identical parents (Eq.16 with m * * = 1).Nevertheless, Fig. 7 is representative of the most common outcomes we have observed for admissible parameter sets (i.e., Turing stable parents, Turing unstable hybrid).

Discussion
This paper presents a multi-generational Turing model of pattern formation in hybrid Mimulus.Model development was motivated by experimental observations of petal phenotypes across three generations: parents, F1, and F2 (Fig. 1).Our model formulation is explicitly diploid, i.e., there is a representation of multiple copies of each genetic loci-two for activator, two for inhibitor (Fig. 3).The model reproduces transgressive phenotypes in F1 hybrids between two unpatterned parents, i.e., red × yellow → red spots on yellow (Fig. 4).Consistent with experiments, simulated selfing of this F1 hybrid often yields a distribution of phenotypes in the F2 hybrid population (Fig. 5).
It is instructive to compare our model with the inheritance of an activator-inhibitor reaction-diffusion system motivated by pattern formation in zebrafish (Miyazawa et al. 2010).This prior work, which is the only other example of a multi-generational Turing model to be found in the literature, presumes that zebrafish phenotypes are controlled by a large number of genetic loci.Assuming additive inheritance, model F1 hybrids use parameter values intermediate to parental values.This approach generated novel phenotypes in F1 hybrids.However, the F1 hybrid retains no information regarding the more extreme parental parameter values; it is therefore unable to reproduce a distribution of phenotypes in an F2 population.In contrast, the multi-generational model presented here allows for loci to be homozygous or heterozygous for parental allele type.Simulated selfing of F1 hybrids leads to nine distinct F2 genotypes and a wide range of phenotypes, consistent with monkeyflower experiments (Fig. 5).This is the preferred approach for a model of the anthocyanin pathway in Mimulus, which is known to be controlled by a small number of genes.
Under the model put forth by Miyazawa et al. (2010), offspring phenotype will necessarily be intermediate to parental phenotype.However, it is well documented that hybridization can lead to offspring phenotype that is more extreme than either parent (see Hochholdinger and Baldauf 2018;Birchler et al. 2010, for reviews of hybrid vigor and heterosis).We find that our explicitly diploid model can yield offspring with phenotypes more extreme than either parent (Fig. 5A).As a limiting case, we performed simulations in which parents are phenotypically identical.These simulations suggest a hypothetical mechanism for the emergence of hybrid phenotypes involving transcription factors whose regulatory binding is weaker between parental alleles versus within parental alleles.In principle, this hypothesis could be empirically tested using ChIP-Seq to identify binding sites and relative strengths of binding (see Liang and Kele 2012, for details).
This paper is a proof of concept that emphasizes a specific biological question, model formulation, and numerical simulations.Each multi-generational simulation presented here is, in fact, nine distinct reaction-diffusion systems and solutions and their numerically calculated steady states.These systems are related to one another in a complex manner through combinations of parameter assignments that account for Mendelian inheritance of multiple alleles and gene products.Because there are no continuously changing (bifurcation) parameters in our representation of inheritance, it is not clear to us what form a deeper mathematical analysis would take.
Possible mechanisms of pattern formation through inheritance will necessarily depend on the model equations and parameters used.Although there is experimental evidence for an activator-inhibitor reaction-diffusion system in Mimulus Ding et al. (2020), the anthocyanin pathway is known to include two accessory proteins, bHLH and WD40, and at least two copies of the inhibitor RTO that may bind and sequester bHLH (Yuan et al. 2014;Ding et al. 2020).In future work, this paper's approach to diploidy and inheritance could be used to study patterning mediated by a gene regulatory network model with different, and more realistic, features.

Appendix A: Equations and Contractions for F2 Genotypes
In the equations for the diploid F2 population, the four-variable system can be contracted by one equation per instance of homozygosity.For the parent with homozygous genotype U A A V A A , the four-variable system (Eq.7) can be contracted to a two-variable system (Eq.8).The equations used for parent B (U B B V B B ) are identical with the replacement of B for A. In contrast, simulations of the doubly heterozygous F1 hybrid with genotype U AB V AB require all four equations (Eq.6).Six additional genotypes are unique to the F2 population.The equations for these genotypes (contracted when possible) are enumerated below.

Fig. 1
Fig. 1 Transgressive petal phenotypes in crosses between Mimulus species.A Inbred (i.e., homozygous) parent flowers with unpatterned petal phenotypes produce reliable patterned phenotypes in F1 hybrids.B Self-pollination of F1 M. l. variegatus × M. cupreus hybrids yield anthocyanin patterns of varying size, density, and intensity of spots (Color figure online)

Fig. 2
Fig.2Schematic diagram of the activator-inhibitor gene regulatory network (haploid).The horizontal lines represent two unlinked genetic loci that code for transcription factors: an activator u and an inhibitor v. Filled circles represent transcription factor binding sites with affinities κ uu , κ uv , κ vu , and κ vv .The parameter α u is the baseline production rate of activator u.The parameter γ u is the increase in production rate that is regulated by the binding of the activator and the inhibitor (magenta curves).Similarly, α v and γ v are the baseline and regulatable production rate of inhibitor v (see Eq. 3) (Color figure online)

Fig. 4
Fig.4Multi-generational simulation of transgressive petal pattern phenotypes.A The diploid reactiondiffusion model (Eq. 4 and Fig.3) produces unpatterned phenotypes when using homozygous parent parameters.Using heterozygous parameters corresponding to the cross between unpatterned parents (P A × P B → F1), the model produces patterned phenotypes in F1 offspring.Yellow-to-red pseudo color represents the total activator concentration (i.e., low to high, u + û).B Simulated breeding of F1 hybrids (F1 × F1 → F2) yields various phenotypes, each corresponding to one of nine F 2 genotypes (U A A V A A , U AB V AB , etc.).These diploid genotypes result from all possible combinations of parent haplotypes (denoted by subscripts A and B) that are composed of two alleles each for activator (U ) and inhibitor (V ) loci.Note: the F2 individuals with genotypes U A A V A A , U AB V AB , and U B B V B B are identical to parent A, the F 1 hybrid, and parent B, respectively.Parameters are as in Table 4 (Color figure online)

Fig. 5
Fig. 5 Four more examples of simulated selfing of transgressive patterned F1 hybrids.A-D are ordered to show increasing pattern variety of the F2 hybrid populations.A No new patterns arise in the F2 generation.B One new F2 pattern is similar to the F1 hybrid.C Three new F2 patterns with spots larger and more intense than the F1 hybrid.D Four new F2 patterns, both labyrinthine and spotted phenotypes with various wavelengths, spot size, and spot intensity (Color figure online) ) where the superscripted A and B are dropped because α * = α A * = α B * (similarly for β and γ ) and κ * * = κ A A * * = κ B B * * .In this equation, m * * is the ratio of trans to cis binding constants, m * * = κ trans * * /κ cis * * , where κ cis * * = κ A A * * = κ B B * * and κ trans * * = κ AB * * = κ B A * * represent affinities within and between parental haplotypes, respectively.

Fig. 7
Fig. 7 Reducing the trans-interaction efficacy affects the critical diffusion ratio d .The trans efficacy of transcription factor binding is the ratio of two affinities (m uu , m uv , m vu , and m vv ).In A the parameter m vu (horizontal axis) is defined as m vu = κ trans vu /κ cis vu where κ cis uv = κ A A uv = κ B B uv and κ trans uv = κ AB uv = κ B A uv represent affinities within and between parental haplotypes, respectively.The critical diffusion coefficient ratio (d = D v /D u , vertical axis) locates a Turing instability (i.e., spotted patterns occur above the blue curve).Reducing the trans efficacy m vu decreases d .Reducing the trans efficacies m uu (B), m uv (C), and m vv (D) impact d with different directions and intensities.Parameters as in Fig. 6 (Color figure online)

Table 1
Assignments of rate constants for P A × P B → F1 hybrid

Table 2
Assignment of binding constants for parents and F1 hybrid superscript pair denotes the allele type of the TF that is binding

Table 3
Assignment of rate and binding constants for the F2 hybrid genotype U AB V B B

Table 4 (
Color figure online)