SymSim: simulating multi-faceted variability in single cell RNA sequencing

The abundance of new computational methods for processing and interpreting transcriptomes at a single cell level raises the need for in-silico platforms for evaluation and validation. Simulated datasets which resemble the properties of real datasets can aid in method development and prioritization as well as in questions in experimental design by providing an objective ground truth. Here, we present SymSim, a simulator software that explicitly models the processes that give rise to data observed in single cell RNA-Seq experiments. The components of the SymSim pipeline pertain to the three primary sources of variation in single cell RNA-Seq data: noise intrinsic to the process of transcription, extrinsic variation that is indicative of different cell states (both discrete and continuous), and technical variation due to low sensitivity and measurement noise and bias. Unlike other simulators, the parameters that govern the simulation process directly represent meaningful properties such as mRNA capture rate, the number of PCR cycles, sequencing depth, or the use of unique molecular identifiers. We demonstrate how SymSim can be used for benchmarking methods for clustering and differential expression and for examining the effects of various parameters on their performance. We also show how SymSim can be used to evaluate the number of cells required to detect a rare population and how this number deviates from the theoretical lower bound as the quality of the data decreases. SymSim is publicly available as an R package and allows users to simulate datasets with desired properties or matched with experimental data.

distribution with a fixed mean of 1 and a standard deviation .
is the within-population σ σ variability parameter and can be set by the user (for the results in this section is set to 0.5) . σ The coordinates of the gene effect vectors can be interpreted as the dependence of its kinetics on the levels of EVFs. For instance, a positive value means that higher concentration of the corresponding EVF can give rise to a higher on rate of a certain promoter (if the EVF and gene effect vectors are both for parameter k on ). The gene effect values are first drawn independently from a standard normal distribution. We then replace each gene effect with a value of zero with probability , thus ensuring that every gene is only affected by a small subset η of EVFs. The sparseness parameter can be set by the user; in this paper we set it to a fixed η value of 0.7.
To ensure that the parameters used for simulation fall into realistic ranges, we estimate the distribution of kinetic parameters of genes from real data (Methods). The estimation is done by fitting a Beta-Poisson distribution to imputed experimental data. As our reference, we used a UMI based dataset of 3,005 cortex cells by Zeisel et al 15  the experimental data). For the purpose of this analysis, the UMI based data was imputed using scVI 4 and the non-UMI data was imputed using ZINB-WaVE 17 (as scVI is only applicable to large data sets). We performed the parameter estimation separately in each of the three largest clusters in the cortex dataset (each cluster is assumed to represent a relatively homogenous subpopulation), and on the entire T cell data (a single condition, which did not contain obvious clusters) and obtained similar distribution ranges (Figures 2b and S1b). These ranges are also in line with observations from other experiments, using smFISH [18][19][20][21][22][23][24][25] or transcription inhibition same cell over time, thus potentially misleading methods for cell state annotation and differential expression. To increase the overall extent of bimodality in the data, we divide (decrease) all k on and k off values by (Figure 2c, yellow arrow). The parameter Bimod can take value from 0 10 Bimod to 1. This way, other properties such as burst frequency ( k on /(k on +k off ) ) and synthesis rate ( s ) remain the same. In Figure 2d we show the effect of varying the Bimod parameter on gene-expression distribution in a simulated homogenous population. Expectedly, as Bimod increases, so does the number of bimodal genes, as well as the average Fano factor ( Figure   S1a).

The second knob: extrinsic variation via extrinsic variability factors (EVFs)
While the first knob focuses on variation within a homogeneous set of cells, the second knob allows the user to simulate multiple, different cell states. This added complexity is achieved by setting different EVF values for different cells, in a way that allows users to control cellular heterogeneity and generate discrete sub-populations or continuous trajectories. To this end, SymSim represents the desired structure of cell states using a tree (which can be specified by the user), where every subpopulation (in the discrete mode) or every cell (in the continuous mode) is assigned with a position along the tree. Different positions in the tree correspond to different expected EVF values, and the expected absolute difference between the value of an EVF of any two cells is linearly proportional to the square root of their distance in the tree (Supplementary Material Section 2).
When SymSim is applied in a discrete mode, the cells are sampled from the leaves of the tree.
The set of cells that are assigned to the same leaf in the tree form a subpopulation, and their EVF values are drawn from the same distribution. As above, we draw these EVF from a normal distribution, where the mean is determined by the position in the tree and the standard deviation is defined by the parameter . When SymSim is applied in a continuous mode, the cells are σ positioned along the edges of the tree with a small step size (which is determined by branch lengths and number of cells; Methods). The EVF values are then drawn from a normal distribution where the mean is determined by the position in the tree, and the standard deviation is defined by (Figure 3a). σ To facilitate the correspondence between EVF values and distances in the tree we use a Brownian motion procedure 27 (Methods; Figure 3a). Specifically, for each EVF we set the mean value at the root of the tree to a fixed number (default set root node to 1) and then perform Brownian motion along the branch. Figure 3a illustrates this process using populations 2 and 3 in the tree as an example. Notably, in the continuous mode, this formulation can give rise to a rich set of patterns of changes in gene expression from root ('progenitor cells') to leaves ('target cells'), including the commonly observed impulse profile 28,29 (Figure S1c-d). As an alternative, we also implemented a mode for simulating continuous data by which gene expression from root to leaves is determined explicitly by an impulse function. This might be preferable if the user would like to generate smoother changes in gene expression, or specific temporal patterns.
In the following analyzes we use the Brownian motion model.
Notably, SymSim only generates a subset of EVFs from the tree, while the remaining ones are drawn from the same distribution for all subpopulations (Figure 3a). The tree-sampled subset, which we term Diff-EVFs (Differential EVFs) represents the conditions or factors which are different between sub-populations, and they usually account for a small proportion of all the EVFs. The number of Diff-EVFs can be set by the user. The results in this section were produced with 60 EVFs, 20% of them are Diff-EVFs.
With this formulation, users can control the extent of between-population variation by setting the branch lengths of the input tree, and combine it with a desired level of within-population variation by setting the parameter . Notably, both and the square root of σ σ branch lengths in the tree are in units of EVF values. It is therefore the case that for any two positions in the tree, the ratio of square root tree distance to determines the separability σ between the respective distributions of the values assigned to any given Diff-EVF (Supplementary Material Section 2). As illustration, Figure 3 depicts the tSNE plots of cells from the same input tree with different in either a discrete ( Figure 3b) or continuous (Figure 3c) σ mode. Notably, both panels show that the tSNE plots reflect the structure of the input tree well.

The third knob: technical variation
A large part of the variation observed in scRNA-seq data sets stems from technical sources [30][31][32] . The technical confounders reflect noise, reduced sensitivity and bias that are introduced during sample processing steps such as mRNA capture, reverse transcription, PCR amplification, RNA fragmentation, and sequencing. In order to introduce realistic technical variation into our model, we explicitly simulate the major steps in the experimental procedures.
We implemented two library preparation protocols: (1) full length mRNAs profiling without the use of UMIs (e.g., with a standard SmartSeq2 33 ); and (2) Profiling only the end of the mRNA molecule with addition of UMIs (e.g., 10x Chromium 34 ). The former protocol is usually applied for a small number of cells and with a large number of reads per cell, providing full information on transcript structure 35 . The latter is normally applied for many cells with shallower sequencing, and it is affected less by amplification and gene length biases 30 .
The workflow of these steps are shown in Figure 4a (Methods). Starting from the simulated true mRNA content of a given cell (namely, number of transcripts per gene, sampled from the stationary distribution of the promoter kinetic model), the first step is mRNA capture, where every molecule is retained with probability . ( Figure S2a-b).
It has been previously shown that estimation of gene-expression levels from full length mRNA sequencing protocols has amplification biases related to sequence-specific properties like gene length and GC-content 30,37 , whereas the use of UMIs can correct these biases 37,38 . In data is from 16 ). In cases where UMIs are used, gene length effects are also modeled during amplification, but these effects are mitigated since each molecule is counted at most once. We therefore do not observe gene length bias in the UMI-based simulated data, similarly to the experimental data ( Figure S2c, real data is from 15 ). Finally, we model batch effects with multiplicative factors that are gene-and batch-specific. In Figure S2d, we show the same population of cells are separated by batches. To simplify the discussion at the remainder of this paper, we assume that the data comes from a single batch.
In Figure 4c, we show the comparison between the simulated true mRNA content of one cell and the observed counts obtained with or without UMI. We consider two scenarios: the first scenario represents a study with a low technical confounding and the second one represents a highly confounded dataset. Parameters which differ between these "good" and "bad" cases in this example include capture efficiency ( α ), extent of amplification bias ( MaxAmpBias ) and sequencing depth ( Depth ). Using "bad" technical parameters introduce more noise to true counts, and compared to the non-UMI simulation the UMIs reduce technical noise. The histograms of true counts and four versions of simulated counts are shown in Figure S2e. Using quantile-quantile plots (Q-Q plots; Figure S2f) further demonstrates that UMIs help in maintaining a better representation of the true counts in the observed data. of genes between simulated and experimental data. Notably, we observe a certain level of inaccuracy in matching the SD at the lower ends, which can be due to lowly expressed genes.

Fitting parameters to real data
Indeed, when we filter out lowly expressed genes, the matching of SD improve substantially (Q-Q plots shown in Figure S3a). Furthemore, we conducted a similar analysis by training Splatter 11 with the same experimental datasets as input, and found that SymSim matches this data significantly better ( Figure S3b). Finally, we inspected the relationship between mean (across all cells) and detection rate (fraction of cells in which the gene is detected) from the SymSim simulations, and observed a similar trend as in the experimental data ( Figure S3c).

Using SymSim to evaluate methods for clustering and differential expression
SymSim can be used to benchmark methods for single cell RNA-Seq data analysis as it provides both observed counts and a reference ground truth. In this section we demonstrate the utility of SymSim as tool for benchmarking methods for clustering and differential expression in a sample consisting of multiple subpopulations, using the structure depicted in Figure 3a σ=0.8 and σ=1, we can tell that when σ is high enough to make the clustering challenging, further increasing σ does not yield obvious changes ( Figure 5b). We observe a similar trend of saturation, inspecting increasing levels of capture efficiency (α), especially with scVI. Comparing the methods to each other, we see that scVI has the highest ARI in most cases and that PCA and SIMLR are comparable with SIMLR being slightly better when the rare population accounts for 5% of all cells and the other way around when the size of rare population increases to 10% of all cells.
Our mechanism for simulating multiple populations automatically generates differentially expressed (DE) genes between populations (in the discrete setting; Figure 3b) or along pseudotime (in the continuous setting; Figure 3c). In the following, we use SymSim to benchmark methods for detecting DE genes, focusing on the discrete setting. We use two criteria to define the ground truth set of DE genes. The first criterion is that the number of Diff-EVFs that are associated with a non-zero gene effect value (which we denote as nDiff-EVFgene; Figure 6a) should be larger than zero . This criterion is motivated by our model of transcription regulation: the kinetic parameters of a gene are affected by extrinsic factors, and changes to extrinsic factors might therefore lead to changes in the number of transcripts.
Indeed, when we compare the true simulated gene expression values between subpopulations (i.e., before introducing technical confounders with the third knob), we get a uniform (random) distribution of p-values for genes with no Diff-EVFs, and an increasing skew as nDiff-EVFgene increases ( Figure 6b, using Wilcoxon test); Figure S4a shows that the log fold change of gene-expression between subpopulations increases with nDiff-EVFgene . Nevertheless, in some cases, the actual expression values of genes with nDiff-EVFgene>0 might not differ since the effects of different Diff-EVFs or the effects of modifying different kinetic parameters may cancel out. Differential expression might also be blurred by a high within-population variability. We therefore added an additional constraint and require that all DE genes must have a sufficiently large log fold difference in their simulated true simulated expression levels (Methods).
An important distinguishing feature of SymSim is that it provides an intuitive way for generating case studies for DE analysis that consist of multiple subpopulations with a predefined structure of similarity. To illustrate this, consider populations 1, 2, and 4 ( Figure 6c), which form a hierarchy (2 and 4 are closer to each other and similarly distant from 1). This user-defined structure is reflected in the sizes of the sets of DE genes, obtained respectively from populations 1 vs 2 (1212 genes), 1 vs 4 (1204 genes) and 2 vs 4 (680 genes). Consistent with the hierarchy, the first two gene sets are overlapping and larger than the third one.
As an example for a benchmark study, we used four methods to detect DE genes: edgeR 40  and t-test improve in their relative performance, compared to edgeR and DESeq2. When increasing capture efficiency, all methods gain performance except for the case of AUROC with 300 cells. In that case, the drop in AUROC for some methods is caused by inflation in p-values as α increases, which results in lower specificity ( Figure S4b). Notably, we noticed that the adjusted p-values from DESeq2 can have many missing entries (NAs), especially when α is low (and thus counts are low), and therefore we used its unadjusted p-values in Figure 6d-e.
However, this assignment of NAs in practice filters out genes which do not pass a certain threshold of absolute magnitude (explained in DESeq2 vignette 42 ). To make use of this filtering, we conducted an additional analysis where we used the adjusted p-values for DESeq2 and compare it to all other methods using only the non-filtered (non NA) genes ( Figure S4c). As expected, the performance of all methods (and specifically DESeq2) improves when considering only this set of genes, and converges to high values already at lower capture efficiency rates.
To summarize, we find that edgeR has the best overall performance, with the t-test rank second followed by Wilcoxon test. We note that the aim of this section is to demonstrate the use of SymSim for methods benchmarking instead of performing a comprehensive comparisons of methods. Nevertheless, our ranking is consistent with results from a recent paper which evaluated 36 methods for DE analysis with single cell RNA-Seq data 43 .

Experimental Design
Deciding how many cells to sequence is a decision many researchers face when designing an experiment, and the optimal number of cells to sequence highly depends on the nature of the biological system under investigation and the respective technical hurdles. A previous approach to this problem 44 assumes that the goal of the experiment is to identify subpopulations of cells and provides a theoretical lower bound for the problem. This bound considers the aspect of counting cells (namely, sequencing enough representative cells from each subpopulation), but it does not account for the identifiability of each subpopulation, which may be hampered by both technical and biological factors as well as the performance of clustering algorithms.
In the following we demonstrate how SymSim can be used to shed more light on this important problem. Importantly, in its current form SymSim does not use real data to model between-population variability. We therefore interpret the results in a relative manner ---how do different variability factors shift the required number of cells, compared to each other and to the theoretical lower bound. Our example focuses on a case of one rare subset, represented by cells from population 2 (using the same tree in Figure 6c; note that one can easily generalize this procedure to multiple rare subpopulations). We simulate observed counts with numbers of cells ( N ) ranging from 600 to 7000. These simulations were based on the parameters fit to the cortex dataset 15 with varying levels of σ and α (100 simulations per parameter configuration).
We applied the same three clustering methods as described in the previous section ( k-means with scVI or PCA and SIMLR). We say that a given algorithm was successful in "detecting the rare population" if at least 50 cells from this set are assigned to the same cluster, and form at least 70% of the cells in that cluster. We use these labels to compute an empirical success probability P for each algorithm and each parameter configuration. To get an upper bound on performance that better reflects the data (rather than the algorithm), we take the maximum P at each configuration, and apply cubic spline smoothing (gray curves, Figure 7a-d).
In each plot we also include the theoretical limit which only requires the presence of at least 50 cells from the rare subpopulation (Methods). The theoretical curve (which is independent of all parameters except N) reaches almost 1 at N=1400. Conversely, the empirical curves vary dramatically, based on parameter values. For an easy case of low within-population variability (σ=0.6) and high capture efficiency (α=0.1) the empirical upper bound curve is close to the theoretical one (Figure 7a). This curve clearly decreases when increasing the effect of either nuisance factor (Figure 7b,c). The reduction is substantially more dramatic when both nuisance factors increase (Figure 7d).
To understand the implications on the number of cells required in a given setting, we calculated how many cells are required, in each configuration, to achieve a success rate of 0.75 ( P =0.75, Figure 7e). As expected, the resulting numbers can be much higher than the theoretical lower bound.

Discussion
SymSim has the following features which are advantageous over existing simulators: (i) We simulate true transcript counts from a kinetic model that can be interpreted in terms of transcript synthesis rate, promoter activation and deactivation. (ii) When generating multiple discrete or continuous populations, instead of generating biological differences through directly altering the true transcript count distribution, we set Diff-EVFs, which can be interpreted as biological conditions that cause the differences between subpopulations of cells. This is a more natural and realistic way to simulate biological transcriptional differences. (iii) The EVF formulation provides an intuitive way to specify and simulate complex structures of cell-cell similarity, without the need for manual specifications of the numbers of DE genes 11 . (iv) When generating observed counts, we simulate key steps in real experimental protocols, which automatically gives us dropout events, length bias, and distribution of library sizes. We also provide choices to use UMI based protocols or non-UMI full length mRNA protocols, as the properties of data output from these two categories can be very different.
The main input parameters to SymSim are self-explanatory with their own biological or technical meanings, which users can adjust to match an experimental dataset of interest. While the procedure of parameter fitting was developed in order to generate simulated datasets with similar properties, it may also provide additional insight, as the parameters are biologically or technically interpretable. For instance, comparing the parameters fit to the the UMI and non-UMI datasets in this study we note that the capture efficiency inferred to the latter is much higher

Estimating kinetic parameters from real data
We estimated kinetic parameters from experimental data using an MCMC approach. For each gene, its expression X depends on p , the proportion of time it is on , and the mRNA synthesis rate s .The parameter p itself is a random variable determined by the kinetic parameters k on and k off . We model p as a Beta distributed variable with shape parameters k on and k off . We model X as a Poisson distributed variable with parameter p*s . The distribution of X is then identical to the distribution calculated using the Master Equation 45 . The downsampling effect is modeled as a Binomial sampling with X being the number of trials, and f being the probability that a transcript is sampled for sequencing.
We fit this model to the experimental data using the Gibbs sampler implemented in RJAGS.
At every iteration, we sample each parameter from its marginal posterior conditional on the value of all other parameters. To meet the assumption that all cells share the same kinetic parameters we divide cells by clustering that is performed in the original study and fit the model to counts in a single cluster of cells at a time. We also use imputed read counts, rather than the raw read counts. We use scVI 4 and ZINB-WaVE 17 for the imputation. Since MCMC is dependent on initial conditions, we fit the model independently three times, for each cell cluster and each imputation method. We thinned the MCMC chain to reduce the effect of autocorrelation, and combined all results to obtain the final distribution of kinetic parameters.

Ranges of kinetic parameters from literature
We look into literature for the ranges of kinetic parameters k on , k off , and s which are experimentally measured [18][19][20][21][22][23][24][25] . The range of burst size, or s , from these studies ranges from 2-4000. And the k on and k off values ranges from 0.0001 to 1 per minute, and the half-life of mRNA varies from 1-10 hours, which correspond to 0.001 to 0.01 per minute. This means that k on / d and k off / d could take values from 10 -2 to 10 3 . The specific parameter values reported by these studies are in Supplementary Material Section 1.

Simulation of discrete and continuous populations
The structure of populations can be represented by a tree and the user can input the tree in length and the number steps is the depth of the tips. For a given tip and a given EVF, the value we sample at the tip is used as the mean of a Gaussian distribution to sample the values for that EVF for all cells in that population with standard deviation σ (Figure 3a).
For the continuous mode we also provide an alternative option to the Brownian motion model: using impulse functions for modeling the path-specific variation. When impulse function is used, for cells sampled from branches that are not on the root-tip path for the specific EVF, they are sampled from a univariate normal with mean equal to the EVF value at their most recent common ancestor with the varying path, and standard deviation σ .

Simulation of technical steps from mRNA capturing to sequencing
We simulate two categories of library preparation protocols, one does not use UMIs (unique molecular identifiers) 34 and sequences full length mRNAs ( using procedures in Smart-seq2 33 as template), and the other uses UMIs and sequences only the 3' end of the mRNA (using the Chromium chemistry by 10x Genomics as template). In the pre-amplification step, we provide option of using linear amplification to mimic the CEL-seq protocol. As shown in Figure 4a, we take one transcript with 16 molecules as an example. To implement the UMIs, each original molecule has a variable to its count at each step. The technical steps include the following: 1) Capturing step: molecules are captured from the cell with probability α.
2) Pre-amplification step: if using non-CEL-Seq protocols, this step involves N rounds of PCR amplifications. We introduce sequence-specific biases during amplification, which includes transcript length bias and other bias assigned randomly. Parameter lenslope can be used to control the amount of length bias, and MaxAmpBias is used to tune the total amount of amplification bias. If using CEL-Seq protocol this step is the in vitro transcription (IVT) linear amplification. Note that for simplicity, this pipeline omits several steps, including reverse transcription, and library cleaning up.

Simulation of amplification biases
During PCR amplification of the full length cDNAs, the PCR amplification rate (namely, probability to be amplified) can vary for different transcripts. As a result, some transcripts are over-and some are under-amplified. This causes the unwanted amplification bias. To simulate this, for each gene, its PCR amplification rate is set to a sum of a basal amplification rate (the input parameter rate_2PCR , which equals to the average amplification rate across all genes) plus a bias term B . The bias term B ranges from -MaxAmpBias to MaxAmpBias , where MaxAmpBias is a user-specified parameter to represent the total amount of amplification bias in our system.
B is composed of two categories of biases: biases related to transcript length (referred to as gene length bias, denoted by ) and biases caused by other factors (denoted by ).

B length B rand
We use a linear function to model the gene length bias: we first bin all gene lengths into nbins bins, and get the average length in each bin: . The length bias term l , l , ..., l ) L = ( bin(1) bin (2) bin (nbins) associated with a gene in bin i is set to: The parameter lenslope controls the extent of gene length bias. To ensure that does B length not exceed MaxAmpBias , the parameter lenslope should be smaller than . We then set the second term to a random value in the axAmpBias/(nbins Therefore, for a given gene with length , its PCR amplification rate is: l This rate is used in all rounds of PCRs in the pre-amplification step. The biases then get amplified as more PCR cycles are performed, where transcripts with higher amplification rate will likely get more molecules. Assigning a UMI to each molecule before amplification allows us to collapse all molecules with the same UMI after amplification, so different amplification rates will not affect the final molecule counts. For Figure 4b, is set to 0.023, is enslope l axAmpBias M set to 0.3, is set to 20, and rate_2PCR is set to 0.7. bins n

Fitting simulation parameters to real data
To find the best matching parameters to a real dataset, we simulate a database of datasets with a grid of parameters over a wide range. For each simulated dataset, we calculate the following statistics: mean, percent-non-zero, standard deviation of genes over all cells. Then given a real dataset, we find the simulated dataset which have the most similar distributions of the statistics to the real data, and return the corresponding parameter configurations.
The parameters and their ranges for simulating the two databases are as following:

Applying dimensionality reduction and clustering methods
We apply three different dimensionality reduction methods to cluster cells simulated from multiple discrete populations: PCA, scVI and SIMLR. PCA is the naive baseline method that is also the most commonly seen in single-cell RNA-seq analysis. scVI is a more recent method that uses a zero-inflated negative binomial variational auto-encoder model to infer latent space for each single cell. For both the first two methods, cluster identities are then assigned using k-means clustering. The third method, SIMLR, performs dimensionality reduction and cluster identity iteratively to maximize cluster separation.

Simulation of differentially expressed (DE) genes
Diff-EVFs give rise to differences between populations as well as DE genes between

Detection of differentially expressed (DE) genes
DE genes in observed counts are detected respectively with edgeR, DESeq2, Wilcoxon test and Student t-test. For edgeR, we used the quasi-likelihood approach (QLF) with cellular detection rate as covariate. For DESeq2, we use "local" for the fittype parameter, and we evaluate its performance respectively based on the p-values and adjusted p-values which serve as filtering of genes.

Binomial model to calculate the probability of detecting a population
Assuming all sequenced cells are correctly assigned to its original population, the probability that at least x cells are detected from a population only depends on the binomial sampling.
Denote the total number of cells by N and the proportion of the cells in the given population by r , the probability that at least x cells are detected for the population is: This formula is used to generate the black curves in Figure 7a-d.

Sources of experimental data
We use two experimental datasets throughout this paper, one is from 16

which does not use
UMIs and the other is from 15 which uses UMIs. The former datasets profiles Th17 cells under various conditions and in our paper we use a subpopulation of 130 TGF-β1+IL-6 cells. We refer to this dataset as the Th17 dataset. The second dataset profiles 3006 cerebral cortex cells. The authors found nine classes in these cells. In this paper, to get distributions of kinetic parameters to map our simulated parameters to the same distribution, we perform parameter estimation respectively on 1) 628 cells sampled from the oligodendrocyte class; 2) 715 cells sampled from the CA1 pyramidal neurons; 3) 296 cells sampled from the S1 pyramidal neurons. To verify that SymSim can simulate data with similar statistics with given experimental dataset, we use all 948 oligodendrocyte cells.         The data is the same as that plotted in Figure 4c.