Introduction

Leaves, petals, and plants are full of vibrant colours and beautiful patterns. Gardeners try their best to defend their plants from botanical diseases. Few would think that viruses can create beauty in our world. However, tulips infected with the tulip breaking virus (TBV) can generate petals with streaks, stripes, or flames (see Fig. 1a–b).

Fig. 1: Broken Tulips.
figure 1

a Tulipa Absalon variety, one of the few maintained “truly” broken bulbs that is striped from the tulip breaking virus rather than genetics. Yellow denotes a loss of pigments, and suggests those regions are infected. Reproduced with permission from The Mount Vernon Ladies’ Association40. b Tulip of the variety Carnaval de Rio. Image taken by Kwang Mathurosemontri, available on unsplash. c A schematic of the components of the TBV model ((1)–(3)). Along the base of the petal are layers of cells, each with varying concentrations of anthocyanin pigment due to the virus consuming the substrate and reducing pigment synthesis.

Full size image

Tulips were brought from Turkey to the Netherlands in the 1590’s by the ambassador Ogier Ghislain de Busbecq1. The Dutch botanist Carolus Clusius is credited with having planted the first tulips in Holland. They were so popular that his garden was frequently robbed until he eventually gave away his collection, spreading tulips throughout the Netherlands1. During his studies, Clusius noticed that some tulips had striped patterns, but these were also smaller, weaker, and less likely to reproduce. It seemed that their beauty had come with frailty, and he suggested they were diseased1.

The beauty of these broken bulbs enchanted the Dutch population during the 17th century, causing a soar in tulip prices and a subsequent plummet. This phenomenon is referred to as Tulipomania, and some economists allude to it as the first recorded financial bubble2,3,4. During the Tulipomania, it was recognized that colour breaking was transmissible from broken tulips to healthy, solid tulips by aphids. The first identified virus, the Tobacco Mosaic Virus (TMV), was discovered in 18925. Only later, in 1928, could Cayley and McKay6 identify the virus that is responsible for the broken tulip patterns, the tulip breaking virus (TBV). This fascinating history of broken tulips is gracefully portrayed by Deborah Leipziger in her poem The Broken Tulip7:

$$\begin{array}{ll}For\,decades,no\,one\,knew \hfill & What\,causes\,tulips\,to\,\,{\hbox{``}}{break} {\hbox{''}}\, \hfill \\ what\,caused\,the\,flaring \hfill & The\,mosaic\,virus \hfill \\ the\,feathering\,of\,tulips, \hfill & carried\,by\,aphids \hfill \\ Parrot\,like, \hfill & infects\,bulbs \hfill \\ Red\,on\,orange \hfill & and\,the\,flower\,breaks \hfill \\ Peppermint\,red\,on\,white \hfill & its\,hold\,on\,one\,color, \hfill \\ Black\,on\,tangerine\,- \hfill & the\,primary\,color\,suppressed \hfill \\ The\,eruption\,into\,flame \hfill & and\,lighter\,colors\,bleeding\,through \hfill \\ for\,broken\,tulips\,like \hfill & \\ Absalom\,and\,Mabel \hfill & the\,beauty\,of\,a\,curse\hfill \end{array}$$

Striping on tulips is not necessarily caused by TBV. Other virus species, like Nepovirus and Potexvirus, and different members of the Potyvirus family, have also been identified on striped tulips8,9. However, the main source for tulip breaking remains potyviruses like TBV, Lily mottle virus, and Rembrandt tulip-breaking virus8,10. Despite the knowledge of these viral strains, the process of how a virus infection of the petal can lead to spatial pattern formation remains to be understood. Karin Moelling11 formulated this as an open question in her book from 2017. She writes on pages 223-224: “Alan Turing ... described the mathematics leading to stripes. ... an activator and a long-range inhibitor are interacting—but I do not know which of the two is the virus”.

The answer to Dr. Moelling’s question is a bit more complicated than expected. We shall see that the pattern forming mechanism in broken tulips is that of an activator-substrate model (also called activator-depletion model)12. Unlike an activator-inhibitor system, we have a long-range substrate instead of having a long range inhibitor. In that context, the virus is the activator, which inhibits anthocyanin expression.

Although patterns are prevalent everywhere in biology, from the stripes on a zebra to the structure of fingers, their root factors are only understood in very limited cases. In 1952, the famous computer scientist and mathematician Alan Turing asserted that some patterns could be caused by chemicals that react and diffuse across a spatial ___domain13,14; specifically, by two chemicals undergoing some type of feedback in their rate of formation and diffusing in such different rates (one faster than the other) as to produce diffusive instability. Proposed reaction-diffusion equations have been successful in mimicking seashells13, fish stripes15, and mammalian skin patterns16,17. In 1977, Wolpert et al. proposed a simple gene activation model that could produce distinct thresholds for patterns in a diffusive environment18. In their model, cells retain positional information by sensing the concentration level along a gradient. From then on, these two theories have been competing with each other. However, given the great variety of biological functions, Turing and Wolpert’s theories have been shown to act in concert in many cases, to produce biological patterns19,20.

Considering animal skin patterns, an additional challenge arose. While an animal grows, its skin surface also grows and changes, and this often occurs on the same time scale as the pattern forming process. Edmund Crampin and collaborators developed methods to apply pattern forming models on domains that change over time21. Their method transforms a system of reaction-diffusion equations on a one-dimensional growing ___domain to a modified system on a fixed ___domain. Our model for the TBV will be based on both the Wolpert-type gradient sensing and the Turing instability mechanisms on a growing ___domain.

Turing instability has already been applied to model the formation of patterns in plants; in particular, to describe the spot formation in the flowers of Monkeyflowers (Mimulus)22,23. Monkeyflowers contain spots on their nectar guides, which consist of high levels of anthocyanin. This example gives great credence to using reaction-diffusion equations to model patterns in plants.

Here, we present a mathematical model that provides a non-linear dynamics explanation for the formation of petal patterns in broken tulips. The model is based on the interplay of virus infection with cell pigmentation gene expressions for anthocyanin and comprises a system of partial differential equations (PDEs), balancing the effects of viral infection, viral use of cell resources, and inhibition of gene expressions for cell pigments. The PDE system includes two famous pattern forming mechanisms, a Wolpert mechanism of pattern formation in spatial gradients18, combined with a Turing instability mechanism14. These dynamics are embedded in a growing ___domain of a tulip petal. We use the methods developed by Crampin et al.21 to introduce a growing ___domain into the PDE formulation, and develop a parametrization of a growing tulip petal to predict the formation of stripe patterns in broken tulips.

Building the TBV model

To develop the model, we describe the spatio-temporal dynamics and biological reactions between three components of an infected tulip petal, namely the pigment, the virus, and the viral resources. The natural pigments in the tulip are represented with the concentration of anthocyanin, T(x, t), the primary pigment in red, black, and purple tulips24. Anthocyanin does not diffuse freely within nor between petal cells. Rather, it is transported upon its synthesis into vacuoles for storage, from which we can see a flower’s bright colours25. The infection is induced by viral components, which consists either of TBV or other viruses, measured in viral load V(x, t), that diffuse across plant cells filled with substrate, whose concentration is described by our last variable, B(x, t). As “building blocks", it represents resources inside cells such as proteins, amino acids, hormones, and nutrients that are used to build new virions. A similar substrate approach is found in Moreno’s work26.

Assuming that only the viral components and the substrate, but not anthocyanin, diffuse across a spatial ___domain [0, L(t)], which grows in time as the plant grows, and that reactions are not dependent on ___location, temperature, or soil acidity, the model describing the dynamics of anthocyanin T(x, t), the virus load V(x, t), and the substrate B(x, t), is given as

$$\frac{\partial T}{\partial t}=\underbrace{\left(\frac{1}{1+{r}_{0}V}\right){\rho }_{0}}_{{{{\rm{synthesis}}}}\,{{{\rm{input}}}}}-\underbrace{{\delta }_{T}T}_{{{{\rm{decay}}}}}+\underbrace{\frac{{\rho }_{T}{T}^{2}}{{k}_{T}^{2}+{T}^{2}}}_{{{{\rm{pigment}}}}\,{{{\rm{autoregulation}}}}},$$
(1)
$$\frac{\partial B}{\partial t}={D}_{B}\frac{{\partial }^{2}B}{\partial {x}^{2}}+\underbrace{\sigma \left(1-\frac{B}{{K}_{B}}\right)B}_{{{{\rm{substrate}}}}\,{{{\rm{production}}}}}-\underbrace{{\rho }_{B}{V}^{2}B}_{{{{\rm{substrate}}}}\,{{{\rm{consumption}}}}},$$
(2)
$$\frac{\partial V}{\partial t}={D}_{V}\frac{{\partial }^{2}V}{\partial {x}^{2}}+\underbrace{{\rho }_{V}\left(\frac{1}{1+{r}_{T}T}\right){V}^{2}B}_{{{{\rm{virus}}}}\,{{{\rm{production}}}}}-\underbrace{\gamma V}_{{{{\rm{decay}}}}},$$
(3)

on a ___domain [0, L(t)] with homogeneous Neumann boundary conditions (see Supplementary Methods—Section 2). The parameters of the model are defined in Table 1 within Supplementary Methods—Section 1.

Let us briefly describe each of the terms in model ((1)–(3)). The production and decay of the anthocyanin T in equation (1) is modelled using the well-known Wolpert model18. It is a simple yet accurate model describing the production of a chemical via DNA synthesis in cells. The parameter ρ0 denotes a constant proportional to the signal level that activates the genes responsible for the biosynthesis of T. As the presence of viruses is suspected to reduce the biosynthesis of anthocyanin27,28, this constant is reduced by a factor \({(1+{r}_{0}V)}^{-1}\). The exact mechanism causing the reduction or inhibition of anthocyanin biosynthesis is not known; however, Lesnaw et al.1 propose that viral gene products may interfere (as a consequence of protein-protein interactions) with activators or suppressors that regulate the biosynthesis of anthocyanin; or that the virus infection may induce co-suppression (posttranscriptional gene silencing; PTGS) of the anthocyanin pathway regulatory genes, a phenomenon first discovered when attempts were made to overexpress the chalcone synthase gene (an enzyme in the anthocyanin biosynthesis pathway) in transgenic petunia29. The parameter δT denotes the natural degradation of anthocyanin. As there is evidence of an autoregulatory feedback loop of transcription factors responsible for the regulation of anthocyanin biosynthesis17,30,31,32,33, we describe a simple autoregulation mechanism for anthocyanin in the last term of equation (1), where ρT denotes the rate of growth that is halved when T  =  kT. In the absence of the virus (V = 0), the tulip naturally regulates the anthocyanin concentration onto a homogeneous and constant homoeostatic value.

For the substrate equation (2), we use a simple logistic growth model with growth rate σ and carrying capacity KB. Since viral replication is proportional to the substrate consumption and cooperative in nature26,34, the viral cooperative interaction appears in both the substrate consumption term of equation (2) and the virus production term of equation (3). This latter term is reduced by a factor \({(1+{r}_{T}T)}^{-1}\) to reflect the anti-viral properties of anthocyanin27,28,35. The virus is naturally degraded at a rate γ. Finally, for the equations of B and V we add diffusion terms. A schematic of the components of the TBV model is depicted in Fig. 1c.

A comprenhensive choice of parameter values for the TBV model is described in the Supplementary Methods—Section 1 and a sensitivity analysis in Supplementary Methods—Section 7. Also, we nondimensionalize the model in Supplementary Methods—Section 2; this allows us to reduce the number of parameters and use the nondimensional model to carry out a multiscale analysis (see Supplemental Methods—Section 3) and numerical simulations (see Supplemental Methods—Section 6).

Incorporating the tulip petal growth into the model

The base of the flower petal can be identified as a one-dimensional interval, [0, L(t)], that lengthens over time and creates layers of a petal that continues to grow during a TBV infection. In order to examine how the growth and curvature of tulip petals evolve, we take a closer look of a tulip petal in Fig. 2a, which shows its fan-like veinal structure. This structure leads us to assume that cells aligned with the same perpendicular arc to the veins have the same age.

Fig. 2: The growth and curvature of tulip petals.
figure 2

a A visible light image of a petal from "Tulipa Seadov" printed on paper. Lines were drawn from evenly spaced points on the 6.35 cm proximodistal (base to tip) axis, and then perpendicular to each vein. b Estimated circular arc lengths of the perpendicular arcs to the petal veins with respect to their position on the proximodistal axis. The fitted logistic curve L(t), given by Equation (4). c The transformation of the intervals [0, L(t)] from straight blue lines into green arcs according to the circular parametrization. d Radii of the circular arcs in (c). A Wolpert function (blue) and a cubic polynomial (orange) were fitted. e Dynamics of the TBV model on the ___domain [0, L(t)] with logistic growth given by Equation (4) over 35 days, with parameters L0 = 0.5 cm, κ = 0.727, ξ = 9.936 and τ = (35 ⋅ 24 h)/6.35 = 132.28 h. Time is given in hours and space is in cm. f Dynamics of the TBV model on the corresponding parametrized arcs (curved green circular arcs in (c) with radii R(t) = − 0.037205(t/τ)3 + 0.454(t/τ)2 − 1.565(t/τ) + 2.682, with the same value for τ). Time applies to the intersection of each circular arc with the center line.

Full size image

Even though tulip petals come in a wide array of shapes and sizes, we choose to illustrate, in Fig. 2a, a tulip petal of the variety "Tulipa Seadov”, and use this triumph tulip as a representative sample. Perpendicular arcs to the petal veins are drawn equidistant along the proximodistal (base to tip) axis of the petal. These arcs are parametrized, for the purpose of generating the simplest growing tulip petal shaped ___domain, as circular arcs (see Supplementary Methods—Section 5 for details). This parametrization allows us to estimate the circular arc lengths, illustrated in Fig. 2b, where the x-axis, denoting the distance along the proximodistal axis, is associated with the growth time. The estimated arc length grows logistically and thus fitted with the function

$$L(t)=\frac{{L}_{0}{e}^{\frac{\kappa t}{\tau }}}{1+\frac{1}{\xi }\left({e}^{\frac{\kappa t}{\tau }}-1\right)},$$
(4)

with parameters κ and ξ, and timescale τ, determined by the time it takes for the tulip flower to start growing and mature. The typical time for this is 30 to 40 days36,37. Assuming that the growth is homogeneous and knowing the length of the proximodistal axis of a sample petal is 6.35 cm, we can create functions of petal growth with respect to time by relating the proximodistal length to the flowering time. These functions return arc lengths that correspond to the lengths of the growing ___domain [0, L(t)] for our TBV model. This correspondence, shown in Fig. 2c, is obtained by estimating the arc radii. Figure 2d displays these radii values and their non-monotonic behaviour fitted with a Wolpert and a cubic polynomial function (see Supplementary Methods—Section 5 for more details).

Emergent stripe patterns on a growing tulip petal

Armed with the analytic function (4) describing the ___domain growth, we can now solve the TBV model (1)–(3) on a growing tulip shaped ___domain. For this purpose, we implement in the TBV model (1)–(3), the procedure developed by Crampin et al.21 to transform a reaction-diffusion model on a growing ___domain [0, L(t)] into a modified model on a fixed ___domain (see Supplementary Methods—Section 4 for details).

To illustrate the results, we start with an initial state where concentration of viral components is randomly distributed with a standard deviation of 0.1 around its kinetic steady state. We let the dynamics of the TBV model evolve first on the growing ___domain of length L(t), indicated as straight blue lines in Fig. 2c, which results in the emergence of spatial patterns of its components shown in Fig. 2e. Thereupon, we identify these spatially patterned dynamics onto the corresponding parametrized arcs (curved green lines in Fig. 2c), and obtain the stripe patterns on a tulip shaped petal ___domain, shown in Fig. 2f.

The simulations carried out so far, which start with an initial condition at the petal base, presuppose that all the new layers of the tulip petal emerging from its base share the same initial state. We call this case, the Forward Model. We can also carry out simulations starting with an initial condition at the outer rim, and consider the case in which a new layer of the TBV model components emerging from the petal base is assumed to be solely influenced by the previous layer. This case is named the Backward Model. Note that even though the simulations for the two models have opposite initial domains (one at the base and the other at the tip of the petal), the oldest part of the tulip petal is the outer rim in both cases. In the forward model, the time lapse of the simulation for every layer of the petal coincides with the chronological growth of the petal; in other words, the initial simulated layer is the base of the petal and the last simulated layer is the outer rim, which then determines the age of the whole petal. In the backward model, the time lapse of the simulation coincides with the appearance timeline of every layer; in other words, the most outer layer, being the first to appear, is the initial simulated layer, whereas the last simulated layer is the base of the petal, which then determines the age of the whole petal.

In Fig. 3, we arranged six simulated tulip petals infected by TBV on different ___domain length sizes for the two above mentioned models, and notice fewer stripes in smaller domains and more stripes in larger domains. This is expected for an activator-substrate mechanism. Another striking result comes from the pattern formation differences between the two models. Observing the region close to the outer rim of the petals, we notice that for the forward case the effect of the viral pattern formation dynamics kicks in as soon as the petal starts growing (petals in Fig. 3a), whereas the effect of the virus on the formation of patterns is felt by the petal some time after for the backward case, resulting in a transient period without any patterns (petals in Fig. 3c). Interestingly, this difference in tulip petal patterns was captured during the Dutch Golden Age by artists who wished to preserve the beauty of broken tulips in their paintings (Fig. 3b, d).

Fig. 3: Simulations of tulip petals infected by TBV.
figure 3

Initial conditions are the model’s kinetic steady state concentrations except for the virus concentration, which is randomly distributed with a standard deviation of 0.1 around its kinetic steady state. a Dynamics of the TBV model on tulip petals with varying initial base lengths that grow forwards and towards the broad rim. This forward dynamic assumes that all the new layers of the tulip petal emerging from its base share the same initial state. The initial length L0 is listed above each petal in centimetres. L0 = 0.526 is realistic for forwards growth, while L0 = 4.80 is realistic for backwards growth. b The tulip Kamelot van Wena, painted by an unknown Dutch Artist in the 17th century. Reproduced with permission from Norton Simon Art Foundation41. c Dynamics of the TBV model on tulip petals with varying initial outer rim lengths that grow backwards and towards the petal base. This backward dynamics assumes that new layers of the TBV model components emerging from the petal base is solely influenced by the previous layer. d The tulip Witte Merveljie, painted by an unknown Dutch Artist in the 17th century. Reproduced with permission from Norton Simon Art Foundation42. There are intriguing similarities between the petals of the paintings in (b and d) and the simulated petals of (a and c), respectively.

Full size image

Turing and Wolpert mechanisms working in unison explain stripe patterns

A thorough multiscale and Turing stability analysis of the TBV model ((1)–(3)) is carried out in Supplementary Methods—Section 3. Based on our comprehensive choice of parameter values, it comes to light that the substrate and the virus act on a faster time scale than the pigmentation. The substrate-virus sub-system is responsible for the formation of patterns, and we can say they form a pre-pattern. The slower T equation then amplifies this pre-pattern by expressing high pigment levels for low virus concentrations and low pigment levels for locations of high virus load. In other words, the Turing model reactions are created upstream of the positional information19.

Discussion

We set out to solve a flower mystery that has served as a source of inspiration for artists and writers since the 17th century, the delightful stripe patterns of broken tulips. This mystery started to reveal itself in 1926, when it was discovered that the stripes are caused by the tulip breaking virus infection of the tulip petals. However, how this infection leads to stripe patterns has been an open question ever since.

In this article, we develop a mathematical model to explain the underlying instability that leads to stripe patterns, thereby answering the above open question and contributing to the full unfolding of this 350 year old mystery.

The TBV model (1)–(3) identifies the dynamic interaction between the tulip breaking virus and its resources as a substrate-virus system responsible for the formation of Turing-like patterns and taking place at a faster time scale than the pigment dynamics, which is described by a Wolpert positional information process. This interplay between these two pattern formation mechanisms (Turing’s and Wolpert’s) acting together on different time scales, and embedded on a growing petal shape ___domain parametrized with respect to the veinal structure growth of tulip petals, results in the emergence of stunning stripe patterns.

In addition to portraying the beauty of broken tulips within a pattern formation frame built with the robust theories of Turing and Wolpert working in unison, we have expanded upon and applied Crampin’s method of modelling growing domains to a real, biological model, namely the TBV model. The model assumes that a petal grows sequentially, layer on layer of cells, and therefore taken as a one dimensional ___domain that grows uniformly along its axis. More applications of growing domains have been long overdue because many biological surfaces are not static. Hopefully this paper can reveal the potential of this powerful method, that may extend to more biological systems, such as the stripes generated by genetically bred tulips.

Also, we hope some modifications to the proposed mathematical model, such as having diffusion terms for all dependent variables or having changes in the expressions of the reaction terms, can lead to models that may apply to more diverse problems of plant patterning.

We assume that the presence of the TBV reduces the biosynthesis of anthocyanin. Some possible mechanisms causing this reduction, such as viral gene products interfering with anthocyanin biosynthesis regulators or the virus infection inducing co-suppression, have been proposed1, but they haven’t been proven. Nonetheless, the proposed mathematical model presented in this article did not require to incorporate a specific mechanism causing the reduction of anthocyanin biosynthesis in the presence of the virus.

Nowadays, “true” broken bulbs infected with viruses are not available for popular distribution on the public market; however, they can be found in nature and also acquired in biological resource centers for potential life science research purposes such as trying to find out the precise biological mechanism leading to the breaking of tulips from an experimental virological perspective. Fortunately, for fans of their beauty, natural striped tulip varieties exist. It is uncertain if the TBV model can be applied to engineered tulip varieties. One can make a case that transgenes38, MYB proteins24, other proteins, or a combination of them can serve as a surrogate breaking virus within natural varieties since they can suppress genes and interfere with pigment regulation. However, natural regulators should not consume the substrate at a rate we expect from viruses.

Finally, the beauty of the patterns generated by the TBV model is astonishing. We could generate a wide variety of artworks by simply changing few parameters and colours. Most Turing patterns are spots or hatches, but branching stripes can be another addition to the repertoire of mathematical art.

Statistics and reproducibility

The values for the lengths and radii of the tulip petal measurements are available on github39. The striped tulip patterns can be reproduced on any commercially available software such as matlab, maple, mathematica etc., using the parameter values outlined in the Supplemental Methods Section 1 and the computer code that is explained in the Supplemental Methods Section 6 and provided on github39.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.