Model equations

We construct a mathematical model where a stable distribution of the spacers in CRISPR arrays emerges as a result of interaction between bacteria and phages. Our model takes into account the key mechanisms of the molecular processes of spacer acquisition and loss as well as host death and replication [2, 3, 13–17, 20, 22–26]. Note that for some bacteria possessing spacers, their CRISPR machinery inside the cell can be non-active (e.g. orphan arrays, loss of cas genes, decay in CRISPR repeats), so the cell will not able to acquire and actively uptake novel spacers. [27] In our model, we always assume that all bacteria with CRISPR can actively use this machinery.

For simplicity, the model assumes that each microbial genome contains a single CRISPR array. Since our main focus is on the statistical distribution of the number of spacers present in the CRISPR array, we do not distinguish here between different types of spacers or their origin. We use a discrete time modelling approach based on Markov chains, following the model flowchart provided in Fig 1. A detailed description of the model is given in the supplementary material. We assume that the maximal number of spacers that the bacteria may have cannot exceed some fixed large number N. We group all individual cells into classes corresponding to the total number of spacers present in their CRISPR array, i ∈ {1, …, N}. The abundance of class i at time t is described by F_i(t). Numbers F_i(t) compose the state vector F(t) = 〈F₁(t), …, F_N(t)〉. We assume the total population size is settled on some constant environmental carrying capacity, i.e. Σ_i F_i(t) = F_const.

Fig 1

Simplified flowchart of the mathematical model: Adjustment of class size i at time iteration t.

The first half-iteration describes the infection in the case the cell meets a phage (with probability p). Possible outcomes of this stage are: gaining j − i spacers (with probability Q_i,j), death of the cell, and survival without gaining any spacers (with probability 1 − p). The second half-iteration mimics the replication stage. Replication can result in loss of i − j spacers in the daughter cell. The parent cell is assumed to always keep the same number of spacers. The loss or conservation of spacers number at this stage are described, respectively, by Δ_j,i and S_i. Details on the parameterization of Q_i,j, Δ_j,i, and S_i are provided in the supplementary material.

Each discrete time step t is split up into two stages (Fig 1). The first stage is related to infection by the phage, which can result in either failure or success of the immune defence; signifying either death of the cell or its acquisition of spacers, respectively. The given stage is modelled by:

The second stage consists of the replication of the surviving cells, continuing until the carrying capacity is reached and the death of those cells whose immune systems failed to resist infection is compensated for. We assume that during replication the parent cell—that inheriting the leading DNA strand—does not lose any spacers. On the contrary, the daughter cell, inheriting the lagging strand of the parent genome, may lose spacers. This is described by the following equation

The elements of the transition matrices A⁽¹⁾ and A⁽²⁾ depend on the number of spacers i as well as model parameters. In the supplementary material, we show that under the model assumptions, Q_i,j, Δ_j,i, and S_i can be parameterised as follows

Our model does not differentiate spacers in regard to their position in the CRISPR array. The assumption is supported by the work by McGinn and Marraffini [30] which shows that spacers located downstream in the CRISPR array can still provide strong immunity when they constitute some critical proportion of the culture infected (in this study we assume that this is always the case).

We made the following assumptions. The number of spacers added during one iteration event is assumed to be exponentially distributed, described by P(μ_k, i − j), P(μ_m, i − j), or P(μ_n, i − j) with mean values given by μ_k, μ_m, and μ_n, respectively. The use of exponential law to describe adding of spacers has some empirical support in the literature, for example, in the study by Li and co-authors [31] the band intensity in the agarose gels is a monotonic function of the number of DNA copies, it seems to behave in a similar way as the exponential distribution, at least qualitatively (see Fig 2B in the cited paper).

Microbial viruses have very high mutation rates and a single mutation in the protospacer or the protospacer adjacent motif (PAM) in the phage genome may allow the virus to escape detection [13, 26]. However, inefficient spacers containing several mismatches might prime new spacer acquisition more efficiently, compared to naive adaptation [32]. In the case of the perfect match of spacer and the phage protospacer we expect successful priming to occur (in this case often referred to as Interference-driven spacer acquisition); this has empirical evidence [25, 33]. In this light, we assume that μ_n ≈ μ_m ≫ μ_k. New spacers are recruited from different regions of the phage genome [16], so that a phage overcoming resistance of a host with an initiated CRISPR array is very unlikely. During replication events we assume that the cell inheriting the leading DNA strand does not lose any spacers. Daughter cells, inheriting the lagging strand of the parent genome, lose spacers according to an exponential distribution with mean value μ_δ(i) = S_Li, where S_L is the fraction of existing spacers that will be lost. The loss of spacers is modelled via an exponential distribution which has partial support in previous studies. For example, the study by Kupczok and Bollback [34] combining a theoretical model and empirical data shows that the maximum of the probability density function (known as the mode of the distribution) of the loss of spacers is observed for either 0 and 1 spacers. Larger numbers of lost spacers are also possible, but at a smaller probability (see also some empirical demonstration of the possibility of loss of multiple spacers [35]). The exponential distribution qualitatively describes this behaviour and we use it here for mathematical convenience.

The parameter values will generally vary across biomes and ecosystems, depending on the abundance and diversity of viruses, mutation and replication rates, etc. We investigated the dependence of the equilibrium probability distribution of the number of spacers i present in CRISPR on the values of these parameters (for mathematical details see the supplementary information). To compare our model with the empirical distributions more appropriately, we address the case of bacterial and archaeal genomes where only a single CRISPR array is present.

Finally, we should highlight that the above mathematical model can also describe a more complicated scenario of the abortive infection mechanism of CRISPR. According to this mechanism, even for the perfect match of spacer, the phage would kill the infected bacteria but will get a reduced burst size [27, 36]. In this case, instead of a single cell, our model will consider a group of neighbouring cells as a combined entity. For such cells the probability to survive the infection will be higher and will be determined by an averaged over space (volume) parameters. In this case, the spatial structure of the microbial community and spatial infection patterns will be taken implicitly. The scale-invariance of the distribution of spacers can justify this interpretations: the entire system behaves similarly to the subsystems even taken at the microlevel.

Metagenomic samples used and collection of spacers

We used all publicly available (3,858) metagenomic datasets from the IMG/M system [37]. We manually classified all samples into 13 distinct habitat types (Table 1) based on the associated metadata, following the criteria described in IMG/VR [38]. 93% of the datasets contained a geographic location (longitude and latitude) which is visualized in Fig 2A. Visualization was done using the Processing programming language (https://processing.org) and a freely available equirectangular projection of the world map (http://eoimages.gsfc.nasa.gov) was used as a background image. Sample points are positioned by latitude and longitude coordinates of Biosamples obtained from GOLD [39].

Fig 2

(a) World map displaying 3,585 out of 3,858 public metagenomic datasets, for which geographic location is reported, and classified into 13 distinct habitat types (see the color code in the caption) from which we extracted the CRISPR spacer information provided by IMG/M. The total number of samples from a single location is shown with different circle sizes, as indicated in the inset box. For more detail see Table 1 The World map taken from https://visibleearth.nasa.gov/images/57752. (b) The presence of power law distributions of spacers within CRISPR arrays of individual metagenomes. Each point represents all the metagenomes containing a given number of CRISPR arrays. The ordinate value shows the proportion of metagenomes for which either a truncated power law or a non-truncated power law was the best fit.

Table 1

Metagenomic datasets classified by habitat types.

Metagenome Environment	Number of Samples
Air	21
Aquaculture	2
Aquatic Sediment	82
Engineered	344
Fresh Water	599
Host (human)	503
Host (other)	289
Host (plant)	275
Marine	747
Non-Marine Saline and Alkaline	152
Terrestrial (soil)	638
Terrestrial (other)	56
Thermal Springs	152
Total	3,858

We used the total number of spacers per contig (CRISPR-Cas array) and per sample generated by the IMG/M system, that uses a modified version of the CRISPR Recognition Tool (CRT) [40] algorithm described in [41].

Complete prokaryotic genomes

All available bacterial and archaeal (6,214) genomes were downloaded from GenBank [42]. We used the CRISPR Recognition Tool (CRT) [40] with default settings to detect CRISPR arrays. As in the case of metagenomes Fig 3B shows the number of CRISPR arrays (y axis) possessing given number of spacers (x axis) regardless of the number of CRISPR-Cas systems in genomes. Fig 3C shows the same data but only for genomes equipped with single CRISPR array.

Fig 3

Spacer distributions in empirical genome sequences generally exhibit power law behaviour.

Statistical distribution of the number of spacers across CRISPR arrays in the data combined from (a) 3,858 metagenomes, (b) 6,214 complete archaeal and bacterial genomes, (c) 237 complete genomes possessing only a single CRISPR array. The curves are fitted using the power law function with an exponential decay p(i) ∼ i^−α e^−λi and are shown by solid lines. The dashed line shows the slope of the power law function. The shaded area shows 95% confidence bands.

A better approach to the distribution of spacers would take into account the distribution of the occurrence of CRISPR systems per genome. This can be done for complete genomes. We should say, however, that the available database of completely sequenced prokaryotic genomes in GenBank seems to be heavily biased towards a limited number of medically important and model organisms [43]. Adjusting the effects of this bias from the empirical data is a non-trivial task which itself can be subjective. Our metagenome dataset seems less biased (e.g. see Table 1). In this case, however, extracting complete genomes would largely reduce the number of points. Nevertheless, distributions in the three cases belong to the same family and are qualitatively identical (see the Results section).

Statistical analysis

Data on the number of spacers in CRISPR arrays in metagenomes, as well as in complete genomes, were handled with Python. We used the Python package powerlaw [44] for analysis of fat tailed distributions in order to adequately fit the patterns observed in natural CRISPR arrays and our numerical simulations to appropriate probability density functions. To generate the results in Figs 3 and 2B, we used the discrete powerlaw.fit object. The goodness of the truncated_powerlaw fit was compared with other candidate heavy tailed distributions using distribution_compare, which calculates Loglikelihood ratio for two given distributions p₁ and p₂ as $$ \mathcal{R}={\sum }_{i=1}^n(\text{ln}[p_1\left(x_i\right)-p_2\left(x_i\right)])$$. P-value is calculated as $$ p=\frac{1}{\sqrt{2\pi n{\sigma }^2}}({\int }_{-\infty }^{-\left|\mathcal{R}\right|}{e}^{\frac{-t^2}{2n{\sigma }^2}}dt+{\int }_{\left|\mathcal{R}\right|}^{+\infty }{e}^{\frac{-t^2}{2n{\sigma }^2}}dt)=\frac{2}{\sqrt{\pi }}{\int }_{\frac{\mathcal{R}}{\sigma \sqrt{2n}}}^{+\infty }{e}^{-t^2}dt$$, where $$ {\sigma }^2=\frac{1}{n}{\sum }_{i=1}^n([p_1\left(x_i\right)-p_2\left(x_i\right)]-[\frac{1}{n}{\sum }_{i=1}^n\text{ln}{p}_1\left(x_i\right)-\frac{1}{n}{\sum }_{i=1}^n\text{ln}{p}_2\left(x_i\right)\left]\right)$$ [21, 44]. Loglikelihood ratios (positive if the data is more likely in the truncated_powerlaw distribution) and significance p-values are listed in Table 2. If p is small (common practice is to consider p < 0.05) then the value for $$ \mathcal{R}$$ is unlikely to be a chance result and its sign can be trusted as an indicator of which distribution fits the data better. Pure power law distribution is a subcase of truncated power law, whereas exponential distribution is a subcase of stretched exponential. Thus the method has to be correct for comparisons of both nested and non-nested distributions. Loglikelihood ratio test is shown to be correct for both cases for the distributions we considered [45].

Table 2

Positive Loglikelihood ratios indicate that the Truncated Power Law fits the empirical data better than the candidate distribution function.

Comparison of the truncated power law p(x) ∼x^−α e^−λx to other heavy-tailed candidate distributions in fitting empirical distribution of spacers.

	Metagenomes		Complete Genomes
	Loglikelihood ratio	P-value	Loglikelihood ratio	P-value
Exponential $$	62981.65	0.0	3877.26	4.90x10−159
Stretched Exponential $$	1306.33	8.23x10−3	35.15	1.17x10−16
Longnormal, μ > 0 $$	2909.08	1.47x10−56	150.64	7.82x10−28
Power Law p(x)∼x−α	536.83	0.0	58.05	0.0

Mathematical model and simulations

Detailed derivation of the expressions for the components of the transition matrices in A^(1,2) is provided in the supplementary material. The iterative numerical simulation was implemented in MATLAB. The maximal number of spacers was considered to be N = 1, 000, the total population size was set as F_const = 10⁶. At each step, the replication rate was multiplied by the scaling factor N_ν to preserve the constant population density at F_const; the value of N_ν is computed using the corresponding formula in the supplementary material. To find the stationary distribution F* we both used direct iterations (after roughly t = 4000 the distribution becomes approximately stationary) and solved the matrix equation F* = F*A = F*A⁽¹⁾ A⁽²⁾. Note that the actual value of F_const does not affect the modelling results: its re-scaling by the factor of 10, 100, etc, does not modify the final distribution of spacers if the simulation time is large enough. Moreover, the statistical distribution (measured as percentage of spacers) would remain the same even in the absence of the upper bound for F_const. Fitting of curves to stationary distributions obtained in the model was done using the Python package powerlaw [44]. The folder S1 Matlab Code contains a MATLAB code which produces the statistical distributions of spacers corresponding to Fig 4 in the Results section.

Fig 4

Statistical distribution of spacers in CRISPR array predicted by the mathematical model with linear parameterisations of q_i and ν_i.

In all graphs, colored lines with dots represent the final distribution F* in the model. The thin, dashed, dotted, and dash-dotted lines show the fitting of the power law distribution with exponential cutoff given by p(i) ∼ i^−α e^−λi. In the labels, likelyhood stands for Loglikelyhood ratio of truncated power law distribution to the second best fit. (a,b) For the fixed parameter values p = 0.7, h = 1, s = 0.4 and μ_n = 10, increasing the fraction of spacers lost during replication from S_L = 1/8 to S_L = 1/7.1 increases the power law exponent from α from 1.03 to 3.1. (c) For the fixed parameter values p = 0.7, h = 0.8, s = 0.2, g = 0.6 and μ_n = 10, decreasing the average number of spacers gained due to viral infection μ_n from 20 to 11.4 increases the power law exponent from 1.91 to 2.57. (d) For the fixed parameter values p = 0.7, h = 0.8, s = 0.2, S_L = 1/6 and μ_n = 14, decreasing the probability of survival in the case of a mutated protospacer from 0.799 to 0.7795 increases α from 1.85 to 2.57. In all of the above cases we assume that 100μ_k = μ_m = μ_n. Decreasing the value of μ_m has a similar effect as decreasing μ_n (see Fig F in S1 Text).

Revealing statistical distribution of spacers from data

We conducted extensive empirical data analysis of the statistical distribution of the number of spacers across CRISPR arrays in 3,858 metagenomes [37] from various geographical locations and ecosystem types (2,189,103 CRISPRs and 11,724,296 spacers in total). The results are presented in (Fig 3A). Fitting the data points with standard curves shows that the distribution follows a power law distribution with exponential decay given by p(i) ∼ i^−α e^−λi, with α = 2.57 and λ = 0.004 ±0.0013. Our estimate for uncertainty of the power law exponent provide the 95% confidence interval to be 2.47 < α < 2.67.

In order to exclude the influence of possible artefacts of sequencing or assembling technique for metagenomic reads, which may corrupt the resulting pattern, we also analyzed the distribution of spacers across CRISPR arrays in all complete archaeal and bacterial genomes available in GenBank (6,214 genomes, 101,471 CRISPRs, 439,829 spacers) [42]. Our results are shown in Fig 3B. We found that the distribution can be approximated by a power law with exponential cutoff for the parameters α = 2.27, and λ = 0.004 ± 0.0014 (the 95% confidence interval for α is 2.20 < α < 2.34). We also plot the distribution of spacers for genomes where only a single CRISPR array is available (237 genomes, 1,094 spacers) [42]. We found that for bacteria and archaea, the power exponent is α = 1.95 ± 0.31 and the exponential decay factor is λ = 0.005 ± 0.0027 (Fig 3C). The apparent disparity with the results obtained for complete genomes with various numbers of CRISPR arrays and metagenomic data as well as a higher confidence interval for α in Fig 3C is likely to be explained by the small sample size.

We examined the goodness-of-fit for the distributions of the exponential, lognormal, stretched exponential and power law type (all frequently used in the literature [44]) and found that the power law distribution with an exponential cutoff fits our data better than the other distribution functions considered, with p-values well below 0.05 (Table 2). The exponential cutoff in the observed power law distribution occurs due to the natural limitations to the number of spacers a CRISPR array may accumulate. Details on the statistical analysis and sample collection are provided in Materials and Methods.

The power law behaviour and scale-invariance in the distribution of spacers in CRISPR can also be seen directly when CRISPR arrays in individual metagenomes are compared to data for the entire planet (Fig 2B). Individual metagenomes typically contain some number of CRISPRs and distribution of spacers per array can be revealed even within single metagenome. We find that in more than two-thirds of the individual metagenomes, the power law gives a better fit than other heavy-tailed distributions. However, our analysis shows that the corresponding p-value was < 0.05 in each comparison for only 7.75% of the datasets. This can be explained as a result of the generally low signal-to-noise ratio when there are too few CRISPRs and spacers within a metagenome.

Finally, we show (Fig 2A) the geographic location of the majority of the available metagenomic datasets (93%). The samples are classified according to the type of the habitats (13 types overall), with the distribution of the metagenomic datasets across the habitat types shown in Table 1. Thus our results regarding distributions of spacers hold for microbial communities regardless of their geographic locations and/or habitat types (see also Table A in S1 Text). As soon as the number of CRISPR arrays in the sample is sufficient for statistical inference, truncated power law distribution is clearly observed, even within a single metagenomic sample. Interestingly, combining data for several metagenomes into a single data set results in the same pattern which can be explained by the fractal nature of the considered distribution type.

Spacer distribution predicted by the model

Using the introduced model, we run iterations for large times t until the distribution of spacers becomes stationary. The eventual stationarity of F in the model is guaranteed since the combined transition matrix A = A⁽¹⁾ A⁽²⁾ is aperiodic and irreducible by construction. The stationary distribution F* is given by the principle eigenvector of A corresponding to the maximal eigenvalue of 1. Here we investigate the possible distributions F* which emerge in various scenarios for the dependence of the replication rate ν_i and the spacer matching probability q_i on the array size i.

We first consider constant values q_i = q and ν_i = 1, i.e. the replication rates and the effectiveness of the immune system do not depend of the size of the CRISPR library. We find that the resultant distribution F* is normal or skew normal in agreement with the central limit theorem. We also observe that the average number of spacers in the array increases with an increase of the average spacer gain μ_n and with a decrease of the fraction of existing spacers lost during replication S_L. The corresponding graphs are shown in Fig D in S1 Text). Including the dependence of q_i on i (we consider a monotonically increasing function q_i) results in a right-centered distribution similar to a normal distribution (supplementary material Fig E(i) in S1 Text). On the other hand, by keeping q_i constant and assuming the replication rate to be a decreasing function of i (e.g. a linear function), we obtain a left-skewed distribution F*; however, the distribution function still rapidly converges and does not exhibit fat-tailed behavior (Figs E(ii), E(iii) and E(iv) in S1 Text).

We find that in order to obtain a power law-like distribution of spacers in the CRISPR array one needs to assume dependence on i in both q_i and ν_i. For example, by considering the simplest (linear) parameterisations q_i = (i − 1)/N and ν_i = (N + 1 − i)/N we observe power law distributions with an exponential cutoff, as shown in Fig 4A and 4B. This observation can be quantitatively verified by fitting the curves to the data points predicted by the model (the estimated exponents α, and λ are shown in the figure). We also use statistical analysis to compare some other candidates for heavy-tailed distributions using the maximum likelihood principle [44]. We find that a power law with an exponential cutoff is the best fit in the case where a left-centered heavy-tailed distribution emerges and persists under small variation of the parameters. The fat tailed behaviour is not strongly dependent upon the concrete parameterisations of q_i and ν_i: the pattern persists for some monotonically increasing functions q_i (e.g. logistic) and monotonically decreasing functions ν_i (e.g., Gaussian centered at zero). The assumptions of monotonicity in the dependence of q_i and ν_i are natural since the possession of a large CRISPR array may incur costs for the cell, and a longer CRISPR library would normally signify a larger biodiversity of spacers, signifying a higher probability of finding an appropriate spacer [16, 20, 46].

Finally, we investigated the dependence of the spacer distribution exponents, α and λ on key model parameters. Overall, we found that a truncated power law distribution can be observed in the model within a large range of model parameters. An increase in the fraction of spacers lost during replication S_L, was seen to generally lead to a gradual increase in the power law exponent α and a decrease in the maximum number of spacers in the population i_max (i.e. the number of non-empty classes F_i, see Fig 4A and 4B). An increase in the average number of spacers gained μ_n generally leads to smaller power law exponents α (Fig 4C), and the same dependence is found for variation of the parameters μ_m and μ_k (Fig K in S1 Text). A reduction in g, the probability that an inefficient spacer does not cause the death of the cell, results in an increase in α (Fig 4D). We should note that variation in the model parameters (S_L, μ_n, μ_m, μ_k, h, g) can result in either gradual or abrupt, bifurcation-like changes in the slope of the spacer distribution. Corresponding examples are provided in the supplementary material (Figs F-K in S1 Text). Finally, by varying key parameters, the model is able to reproduce the CRISPR spacer distributions shown in Fig 3 for metagenomes (α = 2.57), all bacterial genomes (α = 2.27), and bacterial genomes with a single CRISPR array only (α = 1.95). The corresponding graphs are provided in the supplementary material.