Modeling repertoire variation

To describe the stochastic dynamics of an individual clone, we define a probabilistic rule relating its frequency f′ at time t′ to its frequency f at an earlier time t: G(f′, t′|f, t). In this paper, t and t′ will be pre- and post-vaccination time points, but more general cases may be considered. It is also useful to define the probability distribution for the clone frequency at time t, ρ(f ) (Fig 1A).

The true frequencies of clones are not directly accessible experimentally. Instead, sequencing experiments give us number of reads for each clonotypes, n, which is a noisy function of the true frequency f, described by the conditional probability P(n|f ) (Fig 1B). Correcting for this noise to uncover the dynamics of clones is essential and is a central focus of this paper.

Our method proceeds in two inference steps, followed by a prediction step. First, using same-day replicates at time t, we jointly learn the characteristics of the frequency distribution ρ(f ) (Fig 1A) and the noise model P(n|f ) (Fig 1B). Second, by comparing repertoires between two time points t and t′, we infer the parameters of the evolution operator G(f′, t′|f, t ), using the noise model and frequency distribution learned in the first step (Fig 1C). Once these two inferences have been performed, the dynamics of individual clones can be estimated by Bayesian posterior inference. These steps are described in the remaining Results sections. In the rest of this section, we define and motivate the classes of model that we chose to parametrize the three building blocks of the model, schematized in Fig 1: the clone size distribution ρ(f), the noise model P(n|f), and the dynamical model G(f′, t′|f, t).

This method differs from existing approaches of differential expression detection [4, 5] in at least three ways. First, it can explicitly account for the finite count of cells with a given clonotype. That level of description does not exist in differential expression. Second, it follows a Bayesian approach, which gives the posterior probability of expansion of particular clones, rather than a p -value (although see [7] for a recent Bayesian approach to differential expression). Third, it includes information about the clone size distribution as a prior to assess the likelihood of expansion, and can thus extract information about clonal structure and diversity. A detailed description of classical differential expression analysis is given in the Methods section.

Inferring the noise profile from replicate experiments

To study variations arising from experimental noise, we analysed replicates of repertoire sequencing experiments. The tasks of learning the noise model and the distribution of clone sizes are impossible to dissociate. To infer P(n|f), one needs to get a handle on f, which is unobserved, and for which the prior distribution ρ(f) is essential. Conversely, to learn ρ(f) from the read counts n, we need to deconvolve the experimental noise, for which P(n|f) is needed. Both can be learned simultaneously from replicate experiments (i.e. f′ = f), using maximum likelihood estimation. For each clone, the probability of observing n read counts in the first replicate and n′ read counts in the second replicate reads:

where θ_null is a vector collecting all the parameters of both the noise model and the clone size distribution, namely θ_null = {f_min, ν} for the Poisson noise model, θ_null = {f_min, ν, a, γ} for the negative binomial noise model, and θ_null = {f_min,ν, a, γ, M} for the two-step noise model.

While Eq 5 gives the likelihood of a given read count pair (n, n′), we need to correct for the fact that we only observe pairs for which n + n′ > 0. In general, many clones in the repertoire are small and missed in the acquisition process. In any realization, we expect n + n′ > 0 for only a relatively small number of clones, N_obs ⪡ N. Typically, N_obs is of order 10⁵, while N is unknown but probably ranges from 10⁷ for mouse to 10⁸ − 10¹⁰ for humans [14, 15]. Since we have no experimental access to the unobserved clones (n = n′ = 0), we maximize the likelihood of the read count pairs $$ $(n_i,n_i^{\prime })$ , i = 1, …, N_obs, conditioned on the clones appearing in the sample:

While the condition N〈f〉 = 1 ensures normalization on average, we may instead require that normalization be satisfied for the particular realization of the data, by imposing:

where N is estimated as N = N_obs/(1 − P (0, 0)). The first term corresponds to the total frequency of the unseen clones, while the second term corresponds to a sum of the average posterior frequencies of the observed clones. Imposing either Eq 7 or N〈f〉 = 1 yielded similar values of the parameter estimates, $$ ${\hat{\theta }}_{\text{null}}$ .

To test the validity of the maximum likelihood estimator, Eq 6, we created synthetic data for two replicate sequencing experiments with known parameters θ_null under the two-step noise model, and approximately the same number of reads as in the real data. To do so efficiently, we developed a sampling protocol that deals with the large number of unobserved clones implicitly (see Methods). Applying the maximum likelihood estimator to these synthetic data, we correctly inferred the ground truth over a wide range of parameter choices (S1 Fig).

Next, we applied the method to replicate sequencing experiments of unfractioned repertoires of 6 donors over 5 time points spanning a 1.5 month period (30 donor-day replicate pairs in total). For a typical pair of replicates, a visual comparison of the (n, n ′) pairs generated by the Poisson and two-step noise models with the data shows that the Poisson distribution fails to explain the large observed variability between the two replicates, while the two-step model can (Fig 2A–2C). The normalized log-likelihood of the two-step model was slightly but significantly higher than that of the negative binomial model, and much larger than that of the Poisson model (Fig 2D). The two-step model was able to reproduce accurately the distribution of read counts P(n ) (Fig 3A), as well as the conditional distribution P(n′|n ) (Fig 3B), even though those observables were not explicitly constrained by the fitting procedure. In particular, P(n) inherits the power law of the clone frequency distribution ρ(f ), but with deviations at low count numbers due to experimental noise, which agree with the data. Also, the two-step model outperformed the negative binomial noise model at describing the long tail of the read count distribution for clones that were not seen in one of the two replicates (see S2 Fig).

Fig 4 shows the learned values of the parameters for all 30 pairs of replicates across donors and timepoints. While there is variability across donors and days, probably due to unknown sources of biological and methodological variability, there is a surprising degree of consistency. Despite being inferred indirectly from the characteristics of the noise model, estimates for the number of cells in the samples, M, are within one order of magnitude of their expected value based on the known concentration of lymphocytes in blood (about one million cells per sample). Likewise, f_min is very close to the smallest possible clonal frequency, N_cell, where N_cell ≈ 4 · 10₁₁ is the total number of T cells in the organism [16].

The inferred models can also be used to estimate the diversity of the entire repertoire (observed or unobserved). The clone frequency distribution, ρ(f), together with the estimate of N can be used to estimate Hill diversities (see Methods):

In Fig 5 estimates, we show the values, across donor and days, of three different diversities: species richness, i.e. the total number of clones N (β = 0); Shannon diversity, equal to the exponential of the Shannon entropy (β = 1); and Simpson diversity, equal to the inverse probability that two cells belong to the same clone (β = 2). In particular, estimates of N ≈ 10⁹ fall between the lower bound of 10⁸ unique TCRs reported in humans using extrapolation techniques [14] and theoretical considerations giving upper-bound estimates of 10¹⁰ [15] or more [17].

Learning the repertoire dynamics from pairs of time points

Now that the baseline for repertoire variation has been learned from replicates, we can learn something about its dynamics following immunization. The parameters of the expansion model (Eq 4) can be set based on prior knowledge about the typical fraction of responding clones and effect size. Alternatively, they can be inferred from the data using maximum likelihood estimation (Empirical Bayes approach). We define the likelihood of the read count pairs (n, n′) between time points t and t′ as:

where $$ ${\theta }_{\text{exp}}=\{\alpha ,s_0,\bar{s}\}$ characterizes ρ_s(s ) (Eq (4)) with $$ $\bar{s}$ parametrizing ρ_exp(s), and where $$ ${\theta }_{\text{null}}={\hat{\theta }}_{\text{null}}$ is set to the value learned from replicates taken at the first time point t. The maximum likelihood estimator is given by

This maximization was performed via gradient-based methods. In Methods we give an example of an alternative semi-analytic approach to finding the optimum using the expectation maximization algorithm.

In addition to normalization at t, we also need to impose normalization at t′:

with ρ(f′|n, n′) ∝ ∫dfρ(f)G(f′|f)P(n|f)P(n′|f′) is the posterior distribution of the f′ given the read count pair. In practice, we impose Z = Z′, where Z is the normalization of the first time point given by Eq 7. Intuitively, this normalization constraint sets s₀ so that the expansion of a few clones is compensated by the slight contraction of all clones.

We first tested the method on synthetic data generated with the expansion model of Eq 9, with an exponentially distributed effect size for the expansion with scale parameter, $$ $\bar{s}$ :

where Θ(s′) = 1 if s′ > 0 and 0 otherwise. We simulated small, mouse-like and large, human-like repertoires (number of clones, N = 10⁶ and N = 10⁹; number of reads/sample N_reads = 10⁴ and N_reads = 2 ⋅ 10⁶, respectively), using ν = 2 and f_min satisfying N〈f〉_ρ(f) = 1. The procedure consisted of sampling frequencies and log fold factors N times, normalizing by the empirical sum, and then sampling reads from the corresponding measurement distributions, P(n|f ) (see Methods for details). Inference on these data produced a pair of estimates $$ $({\bar{s}}^{*},{\alpha }^{*})$ . For the parameter-free Poisson measurement model, we analyzed the differential expression model, Eq (9), over a range of biologically plausible parameter values. In Fig 6A, we show the parameter space of the inference from two time points of a single mouse repertoire generated with $$ $({\bar{s}}^{*},{\alpha }^{*})=(1.0,{10}^{-2})$ and $$ $s_0=s_0(\alpha ,\bar{s})$ fixed by the normalization constraint Z′ = Z. The sampling procedure was repeated and the set of inferred pair estimates were plotted. The errors are distributed according to a diagonally elongated ellipse (or ‘ridge’), with a covariance following the inverse of the Hessian of the log-likelihood. The imprecision of the parameter estimates is due to the small number of sampled responding clones. With α* = 0.01 and N_obs ≈ 10₄ sampled clones, only a few dozens responding clones are detected. For human-sized repertoires, millions of clones are sampled, which makes the inference much more precise (see Fig 6A, inset).

Once learned, the model can be used to compute the posterior probability of a given expansion factor by marginalizing f, and using Bayes’ rule,

We illustrate different posterior shapes from synthetic data as a function of the observed count pairs in S3 Fig. We see for instance that the width of the posterior narrows when counts are both large, and that the model ascribes a fold-change of s₀ to clones with n′ ⪅ n.

Note that the value of the true responding fraction α is correctly learned from our procedure, regardless of our ability to tell with perfect certainty which particular clones responded. By contrast, a direct estimate of the responding fraction from the number of significantly responding clones, as determined by differential expression software such as EdgeR [13], is likely to misestimate that fraction. We applied EdgeR (see Methods) to a synthetic repertoire of N = 10⁹ clones, a fraction α = 0.01 of which responded with mean effect $$ $\bar{s}=1$ , and sampled with N_read = 10⁶. EdgeR found 6, 880 significantly responding clones (corrected p-value 0.05) out of N_obs = 1,995,139, i.e. a responding fraction 6, 880/1, 995, 139 ≈ 3 ⋅ 10⁻³ of the observed repertoire, and a responding fraction 6, 880/10⁹ ≈ 7 ⋅ 10⁻⁶ of the total repertoire, underestimating the true fraction α = 10⁻².

Inference of the immune response following immunization

Next, we ran the inference procedure on sequences obtained from human blood samples across time points following yellow fever vaccination. To guide the choice of prior for s, we plotted the histograms of the naive log fold-change ln n′/n (Fig 6B). These distributions show symmetric exponential tails, although we should recall that these are likely dominated by measurement and sampling noise. Yet, the difference between the pair of replicates (black) and the pre- and post-vaccination timepoints (red) motivates us to model the statistics of expansion factors as:

with typical effect size $$ $\bar{s}$ . We also tested other forms of the prior (asymmetric exponential, centered and off-centered Gaussian), but they all yielded lower likelihoods of the data (Table 1).

We applied the inference procedure (Eq 10) between the repertoires taken the day of vaccination (day 0), and at one of the other time points (day -7, day 7, day 15, and day 45) after vaccination. Since there are two replicates at each time point, we can make 4 comparisons between any pair of time points. The results are shown in Fig 6C in log-scale.

Same-day comparisons (day 0 vs day 0) gave effectively zero mean effect sizes ( $$ , below the discretization step of the integration procedure), or equivalently α ≈ 0, as expected. Comparisons with other days yielded inferred values of α and $$ $\bar{s}$ mostly distributed along the same ‘ridge’, as observed on synthetic data (Fig 6A), with variations across replicates and donors. The mean effect size $$ $\bar{s}$ is highest at day 15, where the peak of the response occurs, but is also substantially different from 0 at all time points except day 0 (including before vaccination at day −7), with often high values of α. We speculate that these fluctuations reflect natural variations of the repertoire across time, experimental batch effects, as well as biological variability due to differences in the affinities and precursor frequencies of responding clonotypes. As a consequence of the natural diversity, values of the responding fraction α are not learned with great precision, as can be seen from the variability across the 4 choices of replicate pair, and are probably gross overestimations of the true probability that a naive T cell responds to an infection, which is believed to be of order 10⁻⁵ − 10⁻³ [18].

Identifying responding clones

The posterior probability on expansion factors ρ(s|n, n ′) (Eq 13, S4 Fig) can be used to study the fate and dynamics of particular clones. For instance, we can identify responding clones as having a low posterior probability of being not expanding P_null = P(s ≤ 0|n, n′) < 0.025. P_null is the Bayesian counterpart of a p-value but differs from it in a fundamental way: it gives the probability that expansion happened given the observations, when a p-value would give the probability of the observations in absence of expansion. We can define a similar criterion for contracting clones.

To get the expansion or contraction factor of each clone, we can compute the posterior average and median, 〈s〉_n,n′ = ∫ds sρ(s|n, n′) and s_median (F(s_median|n, n′) = 0.5, for the cumulative density function, $$ $F\left(s\right|n,n^{\prime })={\int }_{-\infty }^s\rho \left(\tilde{s}\right|n,n^{\prime })\mathrm{d}\tilde{s}$ ), corresponding to our best estimate for the log fold-change. In Fig 7A, we show how the median Bayesian estimator differs from the naive estimator s_naive = ln n′/n. While the two agree for large clones for which relative noise is smaller, the naive estimator over-estimates the magnitude of log fold-changes for small clones because of the noise. The Bayesian estimator accounts for that noise and gives a more conservative and more realistic estimate.

Fig 7B shows all count pairs (n, n′) between day 0 and day 15 following yellow fever vaccination, with red clones above the significance threshold line P_null = 0.025 being identified as responding. Expanded clones can also be read off a plot showing how both P_null and 〈s〉_n,n′ vary as one scans values of the count pairs (n, n ′) (S5 Fig).

Given the uncertainty in the expansion model parameters $$ ${\theta }_{\text{exp}}=(\bar{s},\alpha )$ , we wondered how robust our list of responding clonotypes was to those variations. In Fig 7C, we show the overlap of lists of strictly expanding clones (P(s ≤ 0|n, n′)<0.025) as a function of θ_exp, relative to the optimal value $$ ${\hat{\theta }}_{\text{exp}}$ (black circle). The ridge of high overlap values exactly mirrors the ridge of high likelihood values onto which the learned parameters fall (Fig 6D). Values of $$ ${\hat{\theta }}_{\text{exp}}$ obtained for other replicate pairs (square symbols) fall onto the same ridge, meaning that these parameters lead to virtually identical lists of candidates for response.

The list of identified responding clones can be used to test hypotheses about the structure of the response. For example, recent work has highlighted a power law relationship between the initial clone size and clones subsequent fold change response in a particular experimental setting [19]. We can plot the relationship in our data as the posterior mean log fold change versus the posterior initial frequency, f (Fig 7D). While the relationship is very noisy, emphasizing the diversity of the response, it is consistent with a decreasing dependency of the fold change with the clone size prior to the immune response.

The robustness of our candidate lists rests on their insensitivity to the details of how the model explains typical expansion. In S3 Fig, we show how the posterior belief varies significantly for count pairs (0, n′), n′ > 0, across a range of values of $$ $\bar{s}$ and α passing along the ridge of plausible models (Fig 7C). A transition from a low to high value of the most probable estimate for s characterizes their shapes and arises as $$ $\bar{s}$ becomes large enough that expansion from frequencies near f_min is plausible, and the dominant mass of clones there makes this the dominate posterior belief. Thus, these posteriors are shaped by ρ_s(s) at low $$ $\bar{s}$ , and ρ(f) at high $$ $\bar{s}$ . Our lists vary negligibly over this transition, and thus are robust to it.

Code

All code used to produce the results in this work was custom written in Python 3 and is publicly available online at https://github.com/mptouzel/bayes_diffexpr.

Normalization of the clonal frequencies

Here we derive the condition for which the normalization in the joint density is implicitly satisfied. The normalization constant of the joint density is

with δ(Z − 1) being the only factor preventing factorization and explicit normalization. Writing the delta function in its Fourier representation factorizes the single constraint on $$ $\stackrel{\to }{f}$ into N Lagrange multipliers, one for each f_i,

Crucially, the multi-clone integral in Eq (15) over $$ $\stackrel{\to }{f}$ then factorizes. Exchanging the order of the integrations we obtain

with $$ $\langle e^{\mu f}\rangle ={\int }_{f_{\text{min}}}^1\rho \left(f\right)e^{\mu f}\mathrm{d}f$ . Now define the large deviation function, $$ $I\left(\mu \right)≔-\frac{\mu }{N}+log\langle e^{\mu f}\rangle$ , so that

Note that I(0) = 0. With N large, this integral is well-approximated by the integrand’s value at its saddle point, located at μ* satisfying I′(μ*) = 0. Evaluating the latter gives

If the left-hand side is equal to 〈f〉, the equality holds only for μ * = 0 since expectations of products of correlated random variables are not generally products of their expectations. In this case, we see from Eq (19) that $$ $\mathcal{Z}=1$ , and so the constraint N〈f〉 = 1 imposes normalization.

Null model sampling

The procedure for null model sampling is summarized as (1) fix main model parameters, (2) solve for remaining parameters using the normalization constraint, N〈f〉 = 1, and (3) starting with frequencies, sample and use to specify the distribution of the next random variable in the chain.

In detail, we first fix: (a) the model parameters (e.g. {α, a, γ, M}), excluding f_min; (b) the desired size of the full repertoire, N; (c) the sequencing efficiency (average number of UMI per cell), ϵ, for each replicate. From the latter we get the mean number of reads per sample, $$ $N_{\text{reads}}^{\text{eff}}=ϵM$ . Note that the actual sampled number of reads is stochastic and so will differ from this fixed value.

We then solve for remaining parameters. Specifically, f_min is fixed by the constraint that the average sum of all frequencies, under the assumption that their distribution factorizes, is unity:

This completes the parameter specification.

We then sample from the corresponding chain of random variables. Sampling the chain of random variables of the null model can be performed efficiently by only sampling the N_obs = N(1 − P(0, 0)) observed clones. This is done separately for each replicate, once conditioned on whether or not the other count is zero. Samples with 0 molecule counts can in principle be produced with any number of cells, so cell counts must be marginalized when implementing this constraint. We thus used the conditional probability distributions P(n|f) = ∑_m P(n|m)P(m|f) with m, n = 0, 1, …. P(n′|f) is defined similarly. Note that these two conditional distributions differ only in their sampling efficiency, ϵ. Together with ρ(f), these distributions form the full joint distribution, which is conditioned on the clone appearing in the sample, i.e. n+ n′ > 0 (denoted $$ $\mathcal{O}$ ),

with the renormalization accounting for the fact that (n, n′) = (0, 0) is excluded. The 3 quadrants having a finite count for at least one replicate are denoted q_x0, q_0x, and q_xx, respectively. Their respective weights are

Conditioning on $$ $\mathcal{O}$ ensures normalization, $$ $P\left(q_{x0}\right|\mathcal{O})+P\left(q_{0x}\right|\mathcal{O})+P\left(q_{xx}\right|\mathcal{O})=1$ . Each sampled clone falls in one the three regions according to these probabilities. Their clone frequencies are then drawn conditioned on the respective region,

Using the sampled frequency, a pair of molecule counts for the three quadrants are then sampled as (n, 0), (0, n′), and (n, n′), respectively, with n and n′ drawn from the renormalized, finite-count domain of the conditional distributions, P(n|f, n > 0).

Using this sampling procedure we demonstrate the validity of the null model and its inference by sampling across the observed range of parameters and re-inferring their values (see S1 Fig).

Computing Fisher information for constrained maximum likelihood problem

The replicate model parameters are θ = (ν, a, γ, log₁₀ M, log₁₀ f_min). Let C(θ) = Z(θ) − 1 be the constraint equation such that we wish to satisfy C(θ) = 0. Let θ* denote the parameters maximizing the likelihood subject to C(θ) = 0. Then the hyperplane orthogonal to the gradient ∇_θ C(θ*) and passing through θ* is the local subspace in which the constraint is satisfied. The projection of Hessian of the log likelihood, H, into this subspace is given by,

where the matrix $$ $P=\stackrel{\to }{n}{\stackrel{\to }{n}}^{\top }$ projects onto $$ $\stackrel{\to }{n}$ , the unit vector co-linear with ∇_θ C(θ*). The inverse of H has one zero eigenvalue; the remaining eigenvalues characterize the Fisher information at the constrained optimum. Error bars for Fig 2 are the projections of the corresponding ellipsoid onto the respective parameter axes.

When computing error bars for the diversities, we use the standard deviation of the statistics of a Monte Carlo estimate of the log diversities obtained via parameter value samples from the multivariate Gaussian approximation of the likelihood using the projected Hessian, $$ $\hat{H}$ .

Comparison to differential expression analysis

Differential expression deals with RNA-seq data, which reports the bulk expression of a large number of genes in a population of cells, and aims to detect significant differences in expression across different populations, either at different times, or under different conditions.

Repertoire sequencing (RepSeq) and expression analysis aim at inferring fundamentally different quantities, although both do it through the number of reads per gene. In differential expression analysis, one is interested in reconstructing the level of expression of particular genes, which are the same in all cells, while in RepSeq one is interested in the number of cells expressing a given clonotype. Thus, in RepSeq the number of transcripts will depend on the number of cells carrying that clonotypes, but also on their expression level, which is assumed to be clonotype-independent but noisy. There are thus three levels of noise in RepSeq: cell sampling noise, expression noise, and mRNA capture noise. By constrast in differential expression there is expression noise, cell-to-cell variability, and capture noise. These sources of noises combine in a different manner than in RepSeq.

edgeR [5], a classical differential expression analysis software, proceeds by learning a noise model using a negative binomial model for expression noise from two identical conditions. Then, comparing RNA-seq data from two datasets, it evaluates a p-value corresponding to the probability that the observed difference in expression between the two datasets has occured just because of noise. We applied edgeR treating each clonotype as a separate gene.

Obtaining diversity estimates from the clone frequency density

For a set of clone frequencies, $$ ${\left\{f_i\right\}}_{i=1}^N$ , the Hill family of diversities are obtained from the Rényi entropies, as D_β = exp H_β, with $$ $H_{\beta }=\frac{1}{1-\beta }ln\left[{\sum }_{i=1}^N{f}_i^{\beta }\right]$ . We use ρ(f) to compute their ensemble averages over f, again under the assumption that the joint distribution of frequencies factorizes. We obtain an estimate for D₀ = N using the model-derived expression, N_obs + P(n = 0)N = N, where N_obs is the number of clones observed in one sample, and $$ $P(n=0)={\int }_{f_{\text{min}}}^{1}P(n=0|f)\rho \left(f\right)\mathrm{d}f$ . For β = 1, we compute exp(N〈 − f log f〉_ρ(f)) and for β = 2, we use 1/(N〈f²〉_ρ(f)).

Differential model sampling

Since the differential expression model involves expansion and contraction in the test condition, some normalization in this condition is needed such that it produces roughly the same total number of cells as those in the reference condition, consistent with the observed data. One approach (the one taken below) is to normalize at the level of clone frequencies. Here, we instead perform the inefficient but more straightforward procedure of sampling all N clones and discarding those clones for which (n, n′) = (0, 0). A slight difference in the two procedures is that N_obs is fixed in the former, while is stochastic in the latter.

The frequencies of the first condition, f_i, are sampled from ρ(f) until they sum to 1 (i.e. until before they surpass 1, with a final frequency added that takes the sum exactly to 1). An equal number of log-frequency fold-changes, s_i, are sampled from ρ(s). The normalized frequencies of the second condition are then $$ $f_i^{\prime }=f_i{e}^{s_i}/{\sum }_j{f}_j{e}^{s_j}$ . Counts from the two conditions are then sampled from P(n|f) and P(n′|f′), respectively. Unobserved clones, i.e. those with (n, n′) = (0, 0), are then discarded.

Inferring the differential expression prior

To learn the parameters of ρ(s ), we performed a grid search, refined by an iterative, gradient-based search to obtain the maximum likelihood. We tested different forms of prior shown in Table 1.

For a more formal approach, expectation maximization (EM) can be employed when tractable. Here in a simple setting, we demonstrate this approach of obtaining the optimal parameter estimates from the data by calculating the expected log likelihood over the posterior and then maximizing with respect to the parameters. In practice, we first perform the latter analytically and then evaluate the former numerically. We choose a symmetric exponential as a tractable prior for this purpose:

with $$ $\bar{s}>0$ , and no shift, s₀ = 0. The expected value of the log likelihood function, often called the Q-function in EM literature, is where $$ ${\bar{s}}^{\prime }$ is the current estimate. Maximizing Q with respect to $$ $\bar{s}$ is relatively simple since $$ $\bar{s}$ appears only in $$ ${\rho }_{\text{exp}}\left(s\right|\bar{s})$ which is a factor in $$ $P(n,n^{\prime },s|\bar{s})$ . For each s, so that $$ $\frac{\partial Q\left(\bar{s}\right|{\bar{s}}^{\prime })}{\partial \bar{s}}={\sum }_{i=1}^{N_{\text{obs}}}{\int }_{-\infty }^{\infty }\mathrm{d}s\rho \left(s\right|n_i,n_i^{\prime },{\bar{s}}^{\prime })\frac{\partial log\left[{\rho }_{\text{exp}}\left(s\right|\bar{s}))\right]}{\partial \bar{s}}=0$ implies so that $$ ${\bar{s}}^{*}=\frac{1}{N_{\text{obs}}}{\sum }_{i=1}^{N_{\text{obs}}}{\bar{s}}_{(n_i,n_i^{\prime })}$ , where

The latter integral is computed numerically from the model using $$ $\rho \left(s\right|n,n^{\prime },{\bar{s}}^{\prime })=P(n,n^{\prime },s|{\bar{s}}^{\prime })/{\int }_{-\infty }^{\infty }P(n,n^{\prime },s|{\bar{s}}^{\prime })\mathrm{d}s$ . Q is maximized at $$ $\bar{s}={\bar{s}}^{*}$ since $$ . Thus, we update $$ ${\rho }_{\text{exp}}\left(s\right|\bar{s})$ with $$ $\bar{s}\leftarrow {\bar{s}}^{*}$ . The number of updates typically required for convergence was small.

The constraint of equal repertoire size, Z′ = Z can be satisfied with a suitable choice of the shift parameter, s₀, in the prior for differential expression, ρ_s(s), namely s₀ = −ln Z′/Z. The latter arises from the coordinate transformation s ← Δs + s₀, and adds a factor of $$ $e^{s_0}$ to all terms of Z′.