Sporulation-competent plasmodial cells of wild type strain LU897 × LU898 (Starostzik and Marwan, 1998) were obtained as previously described (Starostzik and Marwan, 1998; Rätzel et al., 2020). A total of 2.8 gram of plasmodial mass was applied to a 14 cm Ø Petri dish that contained 90 ml of semi-rich Golderer agar (Golderer et al., 2001), based on a salt solution of 0.01% (w/v) niacin, 0.01% (w/v) niacinamide, 0.1% (w/v) CaCO₃, and 0.14 mM CuCl₂, supplemented with 5 g peptone from meat (Sigma Aldrich), 0.75 g yeast extract (Becton, Dickinson & Co.), and 3.9 mM glucose per liter, adjusted to pH 4.6 with concentrated HCl. After starvation for 7 days at 22°C in complete darkness, sporulation was induced with a 15 min pulse of far-red light (λ ≥ 700 nm, 13 W/m²) (Starostzik and Marwan, 1998). Before and at 1-h time intervals after the start of the far-red pulse, samples were taken in duplicate at arbitrarily chosen but distant positions on the plate. Each sample was obtained by picking an agar plug of 1.13 cm² with the cut bulb of a disposable Pasteur pipette (EA62.1; Carl Roth, Karlsruhe, Germany). The plasmodial mass on the agar plug was scraped off with a pipet tip and, by cutting the tip, transferred into a vial of glass beads immersed in liquid nitrogen (Figure 1). After extraction of RNA and removal of contaminating DNA (Marquardt et al., 2017), the relative abundance of the mRNAs of 35 genes, differentiation marker and reference genes (Hoffmann et al., 2012; Supplementary Table 2) was analyzed by gene expression profiling (GeXP), a multiplex RT-PCR method (Hayashi et al., 2007) as previously described (Rätzel and Marwan, 2015; Marquardt et al., 2017).

To correct for differences in the concentration of total RNA and in the efficiency of the RT-PCR reaction, the gene expression values were normalized to the median of the estimated relative concentrations of mRNAs of the 35 genes in each RNA sample. Each normalized expression value was subsequently normalized to the geometric mean of all values obtained for a given gene, and this was performed separately for each gene.

Data were analyzed and processed with a revised and extended pipeline written in R (R Core Team., 2016), based on the previously described script (Rätzel et al., 2020). The normalized gene expression data were clustered and significant clusters were determined with the help of the Simprof algorithm (Clarke et al., 2008) as provided by the clustsig package (Whitaker and Christman, 2014). Expression patterns were visualized in the form of a heatmap generated by the heatmap.2 function, provided as part of the gplots package (Warnes et al., 2016). Changes in gene expression over time were visualized by multidimensional scaling based on Euclidean distance (Gower, 1966) with the help of the cmdscale function provided as part of the stats package v3.5.1 (R Core Team., 2016). Petri nets were constructed from single cell trajectories of gene expression as previously described (Rätzel et al., 2020). Each trajectory is a temporal sequence of gene expression states, where each state corresponds to a Simprof significant cluster. Petri nets specified in ANDL format (Abstract Net Description Language) (Heiner et al., 2013) were imported into Snoopy (Rohr et al., 2010) and graphically displayed by running the Sugiyama layout algorithm (Sugiyama et al., 1981). Petri net places, each representing a gene expression state, were colored according to the relative temporal stability of the expression state or according to the relative frequency with which each gene expression state occurred. Petri net transitions, corresponding to transits between states were colored according to the frequency with which each transit occurred or according to the data subset in which the transit occurred. These parameters were computed and coloring was performed by automatic editing of Snoopy files encoded in xml format, again with the help of a R script. The raw data and the complete computational pipeline used in this study is provided as part of the Supplementary Material.

In previous studies we have shown that the gene expression pattern in samples taken at the same time from different sites of a plasmodium covering a standard Petri dish (9 cm Ø) did not change within the limits of accuracy of the measurements. Accordingly, repeated sampling of the same plasmodial cell yields true time series (Rätzel, 2015; Werthmann and Marwan, 2017). To allow more samples to be taken without consuming too much of the plasmodial mass, we now prepared plasmodia on 14 cm Ø Petri dishes, increasing the surface area covered by the plasmodial mass by 2.4-fold, and took smaller samples by punching agar plugs of 1.13 cm² per sample, to harvest a small portion of the initial total plasmodial mass. To test whether the homogeneity in gene expression is impaired or even lost in the larger plasmodia, we took 9 or 16 samples at the same time from approximately evenly spread sites of a plasmodium, and estimated the gene expression pattern twice in the RNA of each sample, to obtain one technical replicate of each measurement (Aselmeyer, 2019; Driesch, 2019). This allowed to estimate the biological variation in gene expression within a plasmodium as compared to the technical accuracy of the measurements. In order to correct for potential differences in the efficiency of the RT-PCR between reactions, the expression value for each gene was normalized to the median of the expression values of all genes measured in the sample (for details see Materials and Methods). To estimate the technical accuracy of the measurements, the relative deviation of first and second measurement from the mean of the two measurements was estimated for each assayed gene in each of the retrieved plasmodial samples. The frequency distribution of all values was almost symmetrical with a tail consisting of a small number of low values, obviously as a result of inefficient RT-PCR reactions. To estimate the degree of homogeneity in gene expression within a plasmodial cell, we asked to which extent the expression values for the 9 or 16 samples taken from the same plasmodium deviated from their median. To restrict the influence of technical artifacts on the result, we considered the subset of the data where first and the second measurement of the same plasmodial sample deviated not more than two-fold from the mean of the two values. In three of the total of 46 analyzed plasmodia (30 far-red stimulated; 16 dark controls), individual samples deviated from the rest of the samples of the same plasmodium by more than a factor of two. As errors in sample preparation could not be ruled out, these three plasmodia were excluded and the remaining data set of 43 plasmodia (28 far-red stimulated; 15 dark controls) was analyzed taking the mean of 1st and 2nd measurement for each gene in each sample. Among the total of 7160 values, 98% of the symmetric frequency distribution (Supplementary Figure 1) were between 0.48- and 2.10-fold deviation of the median of all values of the respective plasmodium (Supplementary Table 1). There was no obvious difference between dark controls and far-red stimulated plasmodia which were measured at 6 h after the light pulse when genes were already differentially regulated (Supplementary Figure 1 and Supplementary Table 1), indicating that even during the period where the mRNA abundance changed in time, the homogeneity in gene expression levels is maintained. Visual inspection of outliers within the distribution did not reveal any candidates for specific genes that might be inhomogeneously expressed.

In summary, the gene expression values throughout a plasmodium deviated not more than approximately two-fold from the median of all samples from the same plasmodium and were thus within the limits of the technical accuracy of the measurements, even under conditions were genes were in the process of being up- or down-regulated. These differences measured between samples were minor as compared to the differential regulation where the expression level of genes changed in the order of ten to more than hundred-fold (Supplementary Figure 2). These results are consistent with the results of the time-series experiment, where for each time point, two samples were retrieved and analyzed from the same plasmodial cell (see below).

As the assayed genes were evenly expressed and changed evenly in time throughout the large plasmodia, at least within the limits of accuracy of the measurements, we took time series at 1 h time intervals. In order to assay, in each experiment, for the homogeneity and synchrony in gene expression throughout the plasmodium, we took two samples at each time point from different, arbitrarily chosen but distant sites of the plasmodium (Figure 1). In far-red stimulated plasmodia, the first samples (referred to as the 0 h samples) were taken at the start of the experiment, i.e., immediately before application of the 15 min pulse of far-red light. All subsequent samples were taken at 1 h time intervals until 10 h after the start of the experiment [At 5 to 6 h after the far-red pulse cells have passed the commitment point, while visible morphogenesis starts several hours later by entering the transient nodulation stage at about 11 h after the pulse (Hoffmann et al., 2012)]. In the dark controls, the far-red stimulus was omitted. Gene expression in each plasmodial sample taken at a given time point was analyzed twice by GeXP-RT-PCR, where the measurement and the corresponding technical replicate are referred to as 1st and 2nd measurement for sample #1, and 3rd and 4th measurement for sample #2, respectively. Data were normalized as described above.

The technical quality of the measurements was estimated separately for the two data sets, each comprising the data of the samples collected at the 11 time points of the time series. For each plasmodial sample, the relative deviation of the two measurements (1st and 2nd, or 3rd and 4th) from the mean of the two measurements was estimated. The frequency distributions of the deviations and corresponding quantil values indicated that the technical qualities of 1st and 2nd, as well as 3rd and 4th measurement were virtually identical with 95% of the values differing less than a factor of two from each other (Figure 2 and Table 1).

The degree of spatial variability of gene expression within a plasmodium was estimated by combining the data sets for the first and the second sample of a plasmodium taken at each time point of the time series. The frequency distribution of the deviation of each measurement from the mean of 1st, 2nd, 3rd, and 4th measurement of the two samples taken from each plasmodium at any time point was virtually identical to the frequency distributions obtained for the technical replicates, indicating that gene expression within the analyzed plasmodia varied at maximum within the limits of accuracy of the measurements (within a factor of 2 in 95% of the samples). This conclusion is based on the comparison of the quantile distributions of the data sets (Figure 2 and Table 1) considering a total of 36,540 data points.

With this data set, we investigated how expression changes as a function of time in the individual plasmodial cells. The gene expression pattern of a plasmodial cell at a given time point was obtained as the mean of the four expression values of each gene measured in the two plasmodial samples picked at that time point.

For visual representation of the data set and of single-cell trajectories of gene expression, we performed multidimensional scaling (MDS) to obtain a data point for the expression pattern of each cell at each time point. Single-cell trajectories of gene expression are shown in Figure 3. Notably, the gene expression patterns of un-stimulated cells (dark controls) changed as a function of time with the highest variability along coordinate 2 of the MDS plot Figures 3A,B. Trajectories of far-red stimulated cells (Figures 3C,D) moved from the left side to the right side of the plot, while the shape of individual trajectories varied to a certain extent, indicating that the response of the cells was similar though not identical. Obviously, the trajectories of six of the eight far-red-stimulated plasmodia of experiment #1 traversed a considerably larger area of the MDS plot (Figure 3C) as compared to the other stimulated cells, indicating a larger variation in gene expression during the response to the stimulus. The extent of variation is accordingly obvious when the bulk of data points is placed in the same plot (Figure 4A). To search for genes that may account for the scattering along coordinate 2, we visually inspected the individual time series displayed in the form of a heat map (Supplementary Figure 2). In addition to the genes that were clearly up- or down-regulated in response to the stimulus, the messages of four genes, hstA, nhpA, pcnA, and uchA, in the following called pcnA-group genes, changed over time in some of the plasmodia, but there was no obvious consistent relationship to the time point of stimulus application. When only genes were included in the analysis that were clearly up- or down-regulated in response to the stimulus (Supplementary Figure 2, see also Figure 8), the data points of the MDS plot were indeed less scattered (Figure 4B). A qualitatively similar result was obtained plotting the up-regulated and the down-regulated genes separately (Figures 4C,D), with some more variation in the expression of the down-regulated genes. However, expression of the pcnA-group over time (Figure 4E) was clearly different from the up- or down-regulated genes. Expression of the pcnA-group genes was different between cells from experiments #1 and #2 as seen from the trajectories of the cells (Supplementary Figure 3), suggesting that ongoing internal processes in cells of experiment #1 might even influence their response to the light stimulus. Indeed, according to the corresponding MDS plots of the bulk data points (Supplementary Figure 4) the response of the up- and down-regulated genes was less uniformly in cells of experiment #1 (Supplementary Figure 4C) as compared to those of experiment #2 (Supplementary Figure 4D).

For a further analysis, we performed hierarchical clustering of the expression data for all assayed 35 genes, differentiation marker and reference genes (Supplementary Table 2), and identified significantly different clusters of expression patterns with the help of the Simprof algorithm (Clarke et al., 2008; Supplementary Figure 5). The temporal sequences of gene expression patterns classified as Simprof significant clusters defined a trajectory for each individual cell and revealed significant differences between cell trajectories (Table 2). To relate gene expression states and trajectories we constructed a Petri net (bipartite graph) as previously described (Werthmann and Marwan, 2017; Rätzel et al., 2020), by representing each gene expression state by a place and the temporal transit between two states by a transition (Figure 5). A single token marking one place of the Petri net indicates the current gene expression state of a cell. The token moves from its place to a downstream place when the transition, connecting the two places through directed arcs, fires. As each transition is connected to exactly two places (one pre-place and one post-place), tokens are neither formed nor destroyed when moving through the net, so the gene expression state of the cell remains unequivocally defined at any time. The coherent Petri net obtained this way represents a state machine predicting possible developmental trajectories in terms of Markov chains of gene expression states (Rätzel et al., 2020).

The basic modeling principles are summarized in Table 3. We observe the following structural properties of the model which we call ‘Waddington landscape Petri net’:

For technical reasons we add immediate transitions starting alternative trajectories, in order to get a statistical distribution of states in which the experiments have started or will start with a given probability.

In contrast to most state-of-the-art pseudo-time series approaches found in the literature (Saelens et al., 2019), the structure of the Waddington landscape Petri net is not restricted to a partial order, meaning it is neither restricted to a directed acyclic graph nor to a tree. Instead we obtain what is known in Petri net theory as ‘state machine,’ also called in other communities ‘finite state machine’ or ‘finite automata,’ which may involve cycles.

A state machine with one token and its reachability graph, or Markov chain for stochastic Petri nets, are isomorph (i.e., have the same structure, there is a 1-to-1 correspondence); to put it differently: our (stochastic) Petri net represents the Markov chain of states the cells assume in the course of their developmental trajectory and accordingly on their walk through the Waddington landscape. We assume that the Petri net represents the corresponding region of the Waddington landscape predicting possible developmental paths a single cell can follow, which of course yields a state machine.

Representing Markov chains as Petri nets comes with a couple of advantages.

First, Petri nets are equipped with the concept of T-invariants, which belong to the standard body of Petri net theory from very early on (Lautenbach, 1973). We consider T-invariants as crucial in terms of biological interpretation of the generated net structures (Sackmann et al., 2006; Heiner, 2009). The computation of T-invariants is rather straightforward for state machines; due to their simple structure it holds:

Second, modeling the differentiation-inducing stimuli, what we have not done so far, would turn some of the free choice conflicts into non-free choice conflicts, which involves, technically speaking, leaving the state machine net class. To unequivocally identify transits that are stimulus-dependent, we need a higher data density which we will hopefully achieve in one of our next experiments. With stimulus-dependent transitions, the constructed Petri nets and their Markov chains do not coincide anymore, instead the Markov chains as well as the reachability graph are directly derived from the Petri nets and may be analyzed by standard algorithms. Finally, our Petri net approach paves the way for the actual ultimate goal of our future work - reconstructing the underlying gene regulatory networks based on the reachability graphs encoded by the Waddington landscape Petri nets.

Figure 6 displays a Petri net assembled from cell trajectories considering the set of 35 genes. The graphical representation laid out using the Sugiyama algorithm emphasizes the directionality of concurrent processes (Sugiyama et al., 1981). The net indicates that cells started in different states (connected to C0; see Legend to Figure 5 for details) and, after stimulation by far-red light, proceeded to a small set of terminal states (places C71, C76, C77, C78, C79, C95) via multiple, more or less highly connected intermediate states. To facilitate the interpretation of the Petri net, we have colored the places and transitions according to different criteria. In Figure 6A, the transitions being specific to cells of experiments #1 or #2 and for dark controls or light-stimulated cells in the respective experiments, are colored differently. Each place is colored according to the relative frequency of its corresponding gene expression state, indicating that some states occurred more frequently than others. From this representation it is obvious that cells of experiments #1 and #2 form different branches, in part projecting onto different terminal states, reflecting accordingly different developmental trajectories to commitment and sporulation. Figure 6B displays the same Petri net, but with a different color coding. Here, transitions are colored according to how frequent the corresponding transits occurred, indicating that some paths were more frequently taken than others. Places are colored according to their relative stability, defined as the average residence time of a cell in the respective state (see Methods for details). Coloring indicates that cells reached a terminal state through states of different stability, e.g., meta-stable intermediates.

Un-stimulated cells (Table 2, dark controls) spontaneously switched between significantly different states of gene expression. Their trajectories gave three disconnected Petri nets (Figure 7A; the three nets were connected to C0 for technical reasons, see legend to Figure 5).

Petri nets of Figures 6A,B, 7A indicate that the expression pattern in both, stimulated and un-stimulated cells developed predominantly in forward directions while there were some transits back to previous states creating so-called transition-invariants (T-Invariants; Figure 5). We asked whether stimulus-independent temporal expression differences like those observed in the subset of the pcnA-group of genes might have added to this directedness. Therefore, we constructed a Petri net from trajectories based on significant clusters, this time exclusively clustering the subset of up- and down-regulated genes. Basic features found in the Petri nets of up- and down-regulated genes were similar to the ones found in the Petri nets for the full set of 35 genes: Trajectories formed parallel main branches, there were intermediate nodes of different stability, of different connectedness, and hence states that occurred with different frequency (Supplementary Figure 6). In contrast, there was a high number of minimal T-invariants that heavily involved places representing gene expression states that occurred in un-stimulated cells. A Petri net built by considering only transits that occurred in un-stimulated cells (Figure 7B) was nearly covered with T-Invariants, indicating spontaneous, reversible alterations in the expression of up- and/or down-regulated genes. Again, states of gene expression displayed different stability. Considerable variation in the expression level of the up-and down-regulated genes is even most obvious from the heat map of initial states from which trajectories emerged and of terminal states that were observed during the experiment (Figure 8). The high density of T-invariants in the un-stimulated cells suggests that similarly, the T-invariants involving places corresponding to light-stimulated cells are due to gene expression changes that do spontaneously occur before cells are caught by a new attractor formed in response to the far-red stimulus (see Discussion).

Figure 7 also shows that selecting sets or subsets of genes for hierarchical clustering and subsequent Petri net construction may yield Petri nets of different structure delivering accordingly non-redundant information on corresponding subsets. This is also shown in Table 4 for the set of 35 genes and for subsets, the down-regulated, up-regulated, down- and up-regulated, and the pcnA-group of genes. Except for the pcnA-group, the number of places per gene was approximately the same. The number of transitions per gene however became less with more genes considered. This suggests that up-regulation, down-regulation and even expression of the pcnA-group of genes are at least partly coordinated or co-regulated processes. The average number of minimal T-invariants per gene compared for the different groups of genes (Table 4) suggests that reversibility observed for the subsets vanishes when more genes are considered, obviously due to combinatorial effects.

We have argued that Petri net modeling disentangles the complex response of cell reprogramming (Rätzel et al., 2020) and predicts feasible developmental pathways through the Waddington landscape, resulting in significantly distinct single cell trajectories. To reveal similarities and differences of the expression kinetics of individual genes, we plot the geometric mean of the concentration values of the mRNA in a cluster logarithmically, normalized to its concentration at the start of the experiment (t = 0 h) as a function of time for any single cell trajectory. In this kind of plot, the time course of the mRNA of each gene starts at the same point, while the slope of the curve indicates the x-fold change in mRNA abundance over time. Plotting subsets of genes suggests that trajectories through different regions of the Petri net of Figure 6 indeed emerge from qualitatively different expression kinetics, and that genes are also differently regulated relative to each other when different trajectories are compared. In the example shown in Figure 9, pldA is early up-regulated in quite a number of trajectories, followed by pwiA and finally by ligA and rgsA that appear strongly correlated at least in some of the plots. Qualitatively different patterns of regulation relative to each other are also evident for the three phospholipase D-encoding genes (Supplementary Figure 7). The pldA gene is up-regulated while pldB and pldC are down-regulated. In some of the trajectories, the initial change in the concentration of the pldA and pldC mRNAs is inverse as compared to the overall time course. A more comprehensive representation with more genes displayed makes similarities and differences between trajectories even more obvious (Supplementary Figure 8). Here, we observe a phenomenon, which is also seen in Table 2, namely that cells remain in a certain state for some time. This occurs predominantly in the unstimulated cells but is also seen in some of the light stimulated cells, e.g., in those that proceed to state C95 (Figures 6A,B). Cells seemingly are trapped in a meta-stable state (e.g., C96, C69, C100, etc.; Figure 6B) for some time until the developmental program proceeds. We presumably will need more data to see whether this is an artifact which occurs by the discretization of gene expression through clustering. Conversely, discretization might help to identify tipping points for the differential regulation of gene expression as the plots in Supplementary Figure 8 suggest.