1 Introduction

Fig. 1
figure 1

Trajectories of HSV virions for Donor F17 at exogenous antibody concentrations \(0\,\upmu \text {g}/ \text {mL}\) (left) and \(1.0\,\upmu \text {g}/ \text {mL}\) (right). Top row: The displacement of HSV virions in the x-direction. The time indicated in the horizontal axis is shifted for each path so that \(t=0\) corresponds to the moment the path is first observed. Bottom row: All two-dimensional HSV virion trajectories overlaid and plotted in a single frame. For all sub-figures, the trajectory frame rates are 15 observations per second

Full size image

There are several mechanisms by which antibodies (Ab) produced by the immune system can interfere with and even prevent viral infection after an invasion. Antibodies have long been known to bind to surface epitopes on invading virions, rendering the pathogen ineffective either by blocking the epitope from binding to receptors on target cells, or signaling to other immune cells/molecules to inactivate the virus or destroy virus-infected cells. Recent experiments have revealed a previously under-appreciated mechanism: physical hindrance of virion motion and potentially the complete immobilization of virions in mucus secretions that lie on the epithelium (Wang et al. 2014; Newby et al. 2017). Specifically, the presence of virion binding, immunoglobulin G (IgG) antibody, was shown to directly decrease the mobility of the herpes simplex virus (HSV) virions in human cervicovaginal mucus (CVM) (Wang et al. 2014; Schroeder et al. 2018), as well as influenza and Ebola virus-like particles in human airway mucus (Yang et al. 2018). An example of the effect can be seen in Fig. 1, where we display virion trajectories for two populations of HSV virions, originally studied in Wang et al. (2014). The left and right columns show virion movement in the presence of low and high Ab concentrations, respectively. The degree of activity in the low Ab concentration is notably higher.

The possibility of using IgG to hinder the motion of different viruses in mucus provides a novel strategy for immunologists to develop methods to prevent and/or treat viral infection (Newby et al. 2017; Witten and Ribbeck 2017). Population-scale experimental methods have shown that Ab are slightly less mobile in mucus than in phosphate-buffered saline (Olmsted et al. 2001). The reduced diffusivity of Ab in mucus has been attributed to weak transient bonds between individual Ab and the polymeric microstructure of mucus, or “mucin mesh” (Olmsted et al. 2001). Meanwhile, many virions have been shown to diffuse unimpeded in mucus in the absence of a detectable Ab concentration (Olmsted et al. 2001; Wang et al. 2014). For this reason, the observation that virion mobility in CVM is impeded in the presence of Ab (even across the menstrual cycle) implies there must be some physicochemical mechanism at work (Wang et al. 2014; Schroeder et al. 2018).

Fig. 2
figure 2

A schematic depiction of the proposed immobilization process of virions, green circles, by antibodies, blue ‘Y’s, in a mucosal medium. Virions become immobilized when “enough” antibodies are bound to the virions and the mucosal fibers, gray lines. Arrows indicate Ab interacting solely with the mucin fibers. Figure originally presented in Wang et al. (2014)

Full size image

Recently, the authors and collaborators have explored the possibility that Ab can work in tandem with the mucin mesh to hinder diffusing virions. (See Fig. 2 for an idealized schematic of the interactions.) In theory, as a virion diffuses through mucus, an array of Ab can accumulate on its surface. When a sufficient number of virion-bound Ab form low-affinity bonds to the mucin mesh, the virion can become tethered and essentially trapped. This hypothesis was introduced by Olmsted et al. (2001) and confirmed by Wang et al. (2014), by Newby et al. (2017), and by Schroeder et al. (2018). In 2014, Chen et al. (2014) introduced a stochastic/deterministic hybrid model for the immobilization of human immunodeficiency virus (HIV) by IgG in CVM and demonstrated the potential impact of the tandem effect of Ab–virion binding and Ab–mucus transient binding on the ability of viral populations to cross, enter, and pass through a thin mucosal layer. Later, Wessler et al. (2015) used numerical simulations to explore combinations of Ab–virion and Ab–mucus reaction kinetics that produce an optimal effect. Newby et al. (2017) further demonstrated that very-low-affinity Ab–mucus bonds optimize trapping of diffusing nanoparticles using experimental and simulated data along with providing theoretical arguments.

Underlying these mathematical models is a Switching Diffusion Hypothesis: that the chemical reactions responsible for virion (or nanoparticle) immobilization are reversible and, as a consequence, virions should switch between diffusive and immobilized states. When compared to the experimentally observable timescale of 10–20 s, the Ab–mucin kinetic rates are expected to be fast, while the Ab–virion kinetic rates are expected to be slow (see Table 1). It is not clear, however, whether the state switching between diffusion and immobilization should be on a faster or slower timescale than the observable 10–20 s.

In recent modeling efforts, (Chen et al. 2014; Wessler et al. 2015), the effect of surface-bound Ab on the diffusivity of a virion was assumed to be incrementally multiplicative. That is to say, there is a constant \(\alpha \in (0,1)\) such that the diffusivity has the following state-dependent form.

$$\begin{aligned} {\hbox {Incremental Knockdown Hypothesis}: } \,\, D(S(t),N(t)) = \alpha ^{N(t)} D. \end{aligned}$$
(1)

Here D is the diffusivity of the virion in mucus in the absence of Ab, N(t) is the number of virion-bound Ab, and S(t) is the subset of Ab simultaneously bound to the mucin mesh. This reduction in diffusivity is independent of S(t) because the number of simultaneously bound Ab changes so rapidly (relative to the number of bound Ab), the virion only feels the average effect of these changes, which is captured by the number of bound Ab, N(t). The parameter \(\alpha \) can be expressed in terms of the Ab–mucin binding and unbinding rates (\(m_{\text {on}}\) and \(m_{\text {off}}\), respectively) and the effective concentration [M] of binding sites on the surfaces of mucin fibers. If \(m_{\text {on}}[M]\) and \(m_{\text {off}}\) are very large, so that there are many on-and-off switches per second, then an effective diffusivity arises with a so-called knockdown factor \(\alpha = m_{\text {off}}/ (m_{\text {on}}[M] + m_{\text {off}})\) (Chen et al. 2014). In this way, we say that the Incremental Knockdown Hypothesis follows from assuming that the dynamics is in a Fast Switching Regime. That is to say, in this modeling regime, one assumes that [diffusion \(\rightleftarrows \) immobilization] switching is faster than the times between experimental observations and faster than what we plan to use as a simulation time step. We depict a typical trajectory of a virion under this hypothesis in Fig. 3(a). A virion rapidly changes between the immobilized (red) state and freely diffusing states (green). The resulting path has a reduced effective diffusivity that is well approximated by Eq.  1, and the virion exhibits qualitatively less movement than a virion predominately in the freely diffusing state (seen in blue).

Table 1 Parameters and known values incorporated in the model
Full size table
Fig. 3
figure 3

Top row a–c The path of a virion assuming it takes one (green trajectory) or ten (blue trajectory) simultaneously bound Ab for immobilization. Red intervals correspond to periods of immobilization. Bottom row d–f The virion–Ab–mucin dynamics that govern the movement of the simulated virion directly above it. Within each frame, the number of bound Ab N(t) is shown by the purple trajectory and the subset of these Ab that are simultaneously bound to the mucin fibers S(t) assuming a low threshold, \(T= 1\), and higher threshold, \(T = 10\), shown by the green trajectory and blue trajectory, respectively. The binding rate cascade factor c increases from left to right: \(c = 1\), \(c = 20\) and \(c = 200\), respectively. Other model parameters used in the simulation are \(([A]_0, [A]_{\text {exo}}, N_*) = (0.2\upmu \text {g}/ \text {mL}, 0.1\upmu \text {g}/ \text {mL},120)\). The mathematical model is fully described in Sect. 2.4

Full size image

Recent particle tracking experiments now make it possible to analyze virion behavior as it is modulated by various concentrations of Ab (Wang et al. 2014). In Fig. 1, we display two populations of HSV virions diffusing in CVM with 0 \(\upmu \text {g}/ \text {mL}\) and 1 \(\upmu \text {g}/ \text {mL}\) concentrations of exogenous HSV-binding IgG. There is qualitatively less virion movement in CVM with higher concentrations of Ab, but, as we argue below using path-by-path analysis, the trajectories of individual virions appear to resemble either that of a strictly immobilized virion or a strictly freely diffusing virion. This absence of observable switches between immobilized and freely diffusing states might seem to ratify the Fast Switching Hypothesis. However, closer analysis of the freely diffusing particles shows that the diffusivity of freely diffusing virions is essentially the same across all exogenous Ab concentrations. This contradicts the Incremental Knockdown Hypothesis, which predicts the diffusivity should decrease with increasing Ab concentration. While there are essentially no observable switches, and the diffusivity of the free population is not incrementally affected by Ab concentration, we find that the proportion of completely immobilized virions is unmistakably increasing with Ab concentration. (See also Wang et al. 2014.) This suggests an alternative hypothesis: we are in a Slow Switching Regime where switching takes place fast enough (less than the incubation period of thirty minutes) so that the experiments display different movement patterns, but slow enough (more than 20 s) so that we do not see switches in the observational time window.

In this work, we develop and implement the tools necessary for making the preceding claims. To be specific, we use clustering analysis to partition virion paths into a few distinct behavioral patterns. We implement a Bayesian switch-point detection algorithm to assess the prevalence of switches in mobile virions. We develop a Markov chain model for virion–Ab–mucin interactions for use in our characterization of the dependence of virion motility on Ab concentration. A critical feature of this model is the possibility that virion immobilization requires multiple simultaneously surface-bound Ab and that a single virion–Ab–mucin binding event might lead to a cascade of such binding events, which would serve to enhance trapping. Using uncertainty quantification techniques, we explore the limitations of the available data, but argue there is a reasonable parameter regime that is fully consistent with experimental observations.

2 Data Collection, Statistical Methods, and Mathematical Model

2.1 Data Collection

Single particle tracking data of HSV virions were collected from seven different CVM samples at five added doses of exogenously anti-HSV-1 IgG (0, 0.033, 0.1, 0.333, 1.0) \(\upmu \text {g}/ \text {mL}\) with an incubation period of half an hour to one hour. For each sample, the virions were tracked for a duration of 20 s. The x-position and the y-position of all traces were observed at a time interval of \(\delta =1/15s\) . For all the experiments, the fluorescent viruses tracked were the HSV-1 mutant 166v, containing a VP22-GFP tegument protein that serves as an internal fluorophore. The diffusional motion of the fluorescent virions mixed into fresh human cervicovaginal mucus was visualized using an EMCCD camera (Evolve 512; Photometrics, Tucson, AZ) mounted on an inverted epifluorescence microscope (AxioObserver D1; Zeiss, Thornwood , NY), equipped with an Alpha Plan-Apo 100 x /1.46 NA objective, environmental (temperature and CO\(\_2\)) control chamber and an LED light source (Lumencor Light Engine DAPI/GFP/543/623/690). The position time series of virions were recorded using the MetaMorph imaging software resulting in videos (512 \(\times \) 512 pixels, 16-bit image depth) with temporal resolution of 66.7 ms and spatial resolution of 10 nm for 20 s. The x, y position time series were extracted from recorded videos using MetaMorph software. For a more detailed description of the collection process, see the Methods section in Wang et al. (2014).

2.2 Statistical Tools for Virion Trajectory Analysis

We used standard statistical techniques to assess whether the behavior of each virion is consistent with the defining properties of Brownian motion (stationarity with Gaussian-independent increments) and to infer physical parameters.

2.2.1 Test for Gaussianity and Independence of Increments

We used normal quantile–quantile (qqnorm) plots to qualitatively verify that the path statistics are approximately Gaussian. The qqnorm plots for the increment processes had approximately linear relationships for all particles indicating that the x and y increment processes for all particles could be described as Gaussian. To construct such plots, we used the \(\texttt {qqnorm()}\) function in the R programming language found in the stats library.

Noting that if a Gaussian process has uncorrelated increments then the increments are independent, we tested for independence of increments by quantifying the statistical significance of their correlation. Let \(\{U_i(k) := X_i(k\delta ) - X_i((k-1)\delta )\}_{k=1}^n\) and \(\{V_i(k) := Y_i(k\delta ) - Y_i((k-1)\delta )\}_{k=1}^n\) denote the ith particle’s x and y increment processes, respectively. For the ith particle, we estimated the correlation between the x and y increment processes separated h time steps apart using the sample autocorrelation function, \({\mathcal {A}}_i(h; U)\) and \({\mathcal {A}}_i(h; V)\) used in the R programming language. If there are n increments of uniform duration \(\delta \), then for a time lag of \(h\delta \)

$$\begin{aligned} {\mathcal {A}}_i(h; X) := \frac{\tfrac{1}{n} \sum _{k = 1}^{n-h}\big ( U_i ((k+h)) - \overline{U_i}\big )\big ( U_i(k) - \overline{U_i} \big ) }{\tfrac{1}{n}\sum _{j= 1}^{n}(U_i(k) - \overline{ U_i})^2}, \end{aligned}$$
(2)

where \(\overline{U_i} := \frac{1}{n} \sum _{k=1}^n U_i(k)\) (Venables and Ripley 2013). We say the ith particle’s increment processes are anti-persistent (persistent) if both \({\mathcal {A}}_i(h = 1; X)\) and \({\mathcal {A}}_i(h = 1; Y)\) are below (above) the critical value for a \(95\%\) significance level and independent otherwise.

2.2.2 Mean-Squared Displacement

The primary statistical tool for describing a population of microparticle paths is the so-called ensemble mean-squared displacement (MSD), which we denote \(\langle {\mathcal {M}}(t) \rangle \). To calculate it, we first compute a pathwise MSD for each trajectory (denoted \({\mathcal {M}}_i(t)\) for the ith path) and then take an average over these functions. If there are n steps that are uniform of duration \(\delta \), then as defined in Qian et al. (1991),

$$\begin{aligned} {\mathcal {M}}_i(k\delta ) := \frac{1}{n - k + 1} \sum _{j = 0}^{n- k}\big |{X_i((j+k)\delta ) - X_i(j \delta )\big |}^2. \end{aligned}$$

For t between the time points \(\{k \delta \}\), we define \({\mathcal {M}}_i(t)\) by linear interpolation. The slope of the MSD displayed on a log–log scale provides an estimate for each particle’s diffusive exponent, \(\nu \), in the large time regime \(({\mathcal {M}}_i(t) \sim C t^\nu )\). Following standard particle tracking nomenclature, an individual path is said to be Brownian if \(\nu = 1\), subdiffusive if \(\nu \in (0,1)\), and stationary if \(\nu = 0\).

2.2.3 Effective Diffusivity

A fundamental quantity to measure for a Brownian path is its diffusivityD. If (X(t), Y(t)) is the 2d position of the particle at time t, then its diffusivity is defined to be \(D := \lim _{t \rightarrow \infty } {\mathbb {E}}(X^2(t)+Y^2(t)) / 4t\). For a Brownian path with n steps of uniform duration \(\delta \), the maximum likelihood estimator (MLE) for its diffusivity has the form

$$\begin{aligned} D_{\mathrm {eff}}:= \frac{1}{4\delta n} \sum _{j = 1}^{n} \big ( U(j \delta )^2 + V(j \delta ) ^2\big ), \end{aligned}$$
(3)

shown in Appendix A. We refer to \(D_{\mathrm {eff}}\) as the path’s effective diffusivity. We note that this effective diffusivity is only a consistent estimator for D if the path has all the characteristics of Brownian motion, namely stationary, independent, Gaussian increments.

However, as seen in Fig. 4(a)–(c) there are many paths with anti-correlated increments. For such a process, “diffusivity” is not well defined. For those paths that can be described by Brownian motion, \(D_{\mathrm {eff}}\) does not account for observational error. Using the method proposed by Vestergaard et al. (2014), we found a typical variance of the localization error to be minor (of order \(0.01\,\upmu \hbox {m}^2/\hbox {s}\)) as compared to the effective diffusivity (of order \( 1 \upmu \hbox {m}^2/\hbox {s}\)) for such virions, see the SI Section 2 for further details and supporting figure. Nevertheless, we use \(D_{\mathrm {eff}}\) as a descriptor for these paths because this serves the purpose to distinguish between the particles in two different states by the clustering methods described below.

For a given collection of N particles, the ensemble effective diffusivity is the weighted average effective diffusivity of the tracked particles in the sample, denoted \(\langle D_{\text {eff}} \rangle \). When evaluating population statistics in particle tracking experiments, if particle paths are weighted independent of path length, then it has been shown that there is a bias toward highly mobile particles, further discussed in Sect. 2.3.1 (Wang et al. 2015). Based on that analysis, we report the effective diffusivity of an ensemble by taking an average weighted by path lengths. Let \(D_{\mathrm {eff}}^i\) denote the effective diffusivity of ith freely diffusing virion, which has path length \(n_i\). Then

$$\begin{aligned} \langle D_{\mathrm {eff}} \rangle := \sum _{i=1}^N\omega _i D_{\mathrm {eff}}^i \quad \text { where } \omega _i = \frac{n_i}{\sum _{j = 1}^N n_j }. \end{aligned}$$
(4)

2.2.4 Bias-Corrected and Accelerated Percentile (\(BC_a\)) Confidence Interval Method

We constructed confidence intervals for ensemble statistics based on the bootstrapping \(BC_a\) method due to its second-order accuracy and invariance under transformations. See Efron and Tibshirani (1994) for the formulation of confidence intervals using this method. We used the boot.ci() function in the R programming language found in the in the boot library to obtain the \(BC_a\) confidence intervals for the ensemble statistics as follows. First, we simulated 10,000 booted samples (with replacement) from an ensemble of N tracked particles weighted by the particle path lengths. The \(BC_a\) confidence interval is then the usual confidence interval constructed using this population of (weighted) bootstrap samples.

2.3 Classification Scheme for Virion Paths

For each donor and concentration, we employed a hierarchical clustering algorithm to separate the HSV virions into distinct clusters based on a set of pathwise statistics: x-increment ACF, y-increment ACF, and \(\log 10\) transform of the effective diffusivity. We defined the dissimilarity measure between pairs of virions i and j by a weighted Euclidean distance d(i, j) with weights of 1 / 4, 1 / 4, and 1 / 2 for the differences in \({\mathcal {A}}_i(1; X)\), \({\mathcal {A}}_i(1; Y)\), and \(\log 10(D_{\mathrm {eff}})\), respectively. The dissimilarity measure between clusters was set to be the average linkage. That is to say, the dissimilarity between clusters R and Q is defined to be

$$\begin{aligned} d(R,Q) = \frac{1}{|R||Q|}\sum _{i \in R, j \in Q} d(i,j). \end{aligned}$$
(5)

Hierarchical clustering is an agglomerative clustering method (Kaufman and Rousseeuw 2009). The algorithm is initialized by setting each data point as a distinct cluster. During each iteration, clusters are merged together to minimize the dissimilarity between all clusters. The algorithm stops when all data points are in a single cluster. This process is depicted graphically through the dendrogram where clusters merge at a height equal to dissimilarity between them. We obtained the k cluster by cutting the resulting dendrogram at the uniform height yielding k clusters.

Because hierarchical clustering is an unsupervised clustering method, in which the number of clusters is not known a priori, the number of clusters has to be chosen by the practitioner. One such method of obtaining the number of clusters is called the elbow method. In this approach, the within-sum-of-squares values of the clusters (WSS) is computed and plotted for a range of cluster numbers. The WSS decreases with the number of clusters, and typically there is a bend or “elbow” in the graph that guides the selection of an appropriate number of clusters.

In all cases, there was a major drop in WSS from one to two clusters, but it was rarely clear how to specify the elbow among \(k = 2, 3,4\) or 5 clusters. We chose to use \(k = 4\) in almost all cases because the results were consistent with our biophysical intuition that there might be Freely Diffusing, Immobilized, Subdiffusive, and Outlier states. A few examples of each class are displayed in Fig. 1 of supplemental materials. Although we allowed for four clusters when labeling each cluster with a biological classification, clusters were typically merged together as the Subdiffusive and Outlier class were few in number.

2.3.1 “Frame-by-frame” Method to Compute Empirical Distribution of Each Cluster

It has been shown in Wang et al. (2015) that fast-moving particles are overestimated on shorter timescales in 2d particle tracking. This bias toward the fast-moving particles arises due to individual fast particles leaving and reappearing in the focal plane as distinct traces and to new particles entering and leaving the focal plane throughout the duration of the experiment. To minimize overestimating the freely diffusing population, we employed the “frame-by-frame” method developed in Wang et al. (2015) to compute the fraction of each population present in the data. The “frame-by-frame” method assigns each tracked particle a weight based on the number of frames the particle appears in the field of view, whereas in the conventional method each particle has the uniform weight of one. Under this weighting system, for a sample of size N, the weighted sample proportion of the ith state is given by

$$\begin{aligned} \hat{p_i} = \sum _{k = 1}^N\omega _k \delta _{ik} \qquad \text {for }\quad \omega _k = \frac{n_k}{\sum _{k = 1}^N n_i} \end{aligned}$$
(6)

where \(\delta _{ik}\) is the Kronecker delta function.

2.4 Mathematical Model for Asymptotic Probability of Immobilization

We mathematically model the dynamics of a virion under the Switching Diffusion Hypothesis by the following SDE:

$$\begin{aligned} \mathrm {d}\varvec{X}(t) = \sqrt{2 D\big (N(t), S(t)\big )} \mathrm {d}\varvec{W}(t) \end{aligned}$$
(7)

where \(\varvec{W}(t)\) is standard 2d Brownian motion and the state-dependent diffusivity, D(N(t), S(t)), depends on two time-dependent processes: N(t), the number of antibodies bound to the surface of a focal virion at time t, and S(t), the subset of these antibodies simultaneously bound to mucin binding sites at time t. We establish a threshold parameter T. A virion is defined to be immobilized if there are at least T simultaneously bound antibodies, \(S(t) \ge T,\) and defined to be freely diffusing if there are fewer than T simultaneously bound antibodies, \(S(t)< T. \) Under this convention, the time-dependent diffusivity is given by

$$\begin{aligned} D\big (N(t), S(t)\big ) = {\left\{ \begin{array}{ll} D &{} 0 \le S(t) < T\\ 0 &{} T \le S(t) \le N(t) \end{array}\right. } \end{aligned}$$

where the constant D is the diffusivity of the virion in mucus in the absence of Ab. In the following sections, we present a mathematical model that describes the asymptotic probability of the immobilized state when exposed to varying exogenous antibody concentrations.

2.4.1 Model Assumptions

Based on the initial population clustering analysis, there appears to a subpopulation of virions that do not interact with the antibodies. We define q to be the probability that a given virion will interact with the Ab population. Second, for the sake of simplicity, we assume that Ab–virion binding sites operate independently from each other. However, we allow for cooperativity among the Ab in binding to the mucosal environment. Once the virion has T simultaneously bound Ab–mucin–virion interactions \((S \ge T)\), the surface-bound antibodies might bind to the mucin fibers differently than if the virion was freely diffusing. We parameterize this by a multiplicative change in Ab–mucin binding rate through the introduction of the dimensionless parameter c. If \(c >1\), the parameter has a cascade effect, aiding in the immobilization process (D’Orsogna and Chou 2009; Goychuk et al. 2014; Holcman and Schuss 2014; Matsuda et al. 2014, Newby et al. 2017).

2.4.2 A Markov Chain Model for Virion–Ab–Mucin Dynamics

Let \(N_*\) denote the number of independent Ab binding sites on the surface of an HSV virion. Antibodies bind and unbind from these sites at rates \(k_{\text {on}}\) and \(k_{\text {off}},\) respectively, with dissociation constant \(k_{\text {d}}:= k_{\text {off}}/k_{\text {on}}\).

Virion-surface-bound antibodies interact with the surrounding mucosal medium, binding to and unbinding from mucin binding sites, at rates \(m_{\text {on}}\) and \(m_{\text {off}}\), with dissociation constant \(m_{\text {d}}:= m_{\text {off}}/ m_{\text {on}}.\) The total Ab concentration [A] is the sum of the exogenous \([A]_{\text {exo}}\) and endogenous \([A]_0\) Ab concentrations, and the total concentration of binding sites on mucin fibers is denoted [M]. See Table 1 for a comprehensive list of variables.

We model the Ab–virion interactions using a continuous-time Markov Chain (CTMC) assuming linear state transitions. If a given virion has n occupied (Ab-bound) surface binding sites at time t, then the CTMC transition rates are given by:

$$\begin{aligned} n \underset{nk_{\text {off}}}{{\mathop {\rightleftarrows }\limits ^{(N_*-n) k_{\text {on}} [A]}}} n + 1. \end{aligned}$$
(8)

If there are s simultaneously bound Ab cross-linking the virion to mucin fibers at time t and n occupied virion-surface-binding sites, then the conditional Ab–mucin dynamics are modeled by a CTMC with state transition rates

$$\begin{aligned} s \underset{s m_{\text {off}}}{{\mathop {\rightleftarrows }\limits ^{(n-s) g(s) m_{\text {on}} [M]}}} s + 1 \end{aligned}$$
(9)

for \(s \le n\), where

$$\begin{aligned} g(s) = {\left\{ \begin{array}{ll} 1 &{} s < T \\ c &{} s \ge T. \end{array}\right. } \end{aligned}$$
(10)

The function in Eq. 10 quantifies the impact immobilization has on the rate at which additional antibodies cross-link to the mucin fibers, i.e., the binding cascade effect, and results in a nonlinear transition rate when \(c \ne 1\). We note that the transition \((n,s) \rightarrow (n-1, s-1)\) is omitted from our analysis to facilitate with explicit likelihood calculations. This does not qualitatively affect our results because of the timescale separation between Ab–virion and Ab–mucin kinetics. Since we assume that (conditioned on number of Ab–virion bindings N) the stationary distribution of S is achieved rapidly, the initial “error” using \(S = s\) instead of \(S = s-1\) after an \(N = n \rightarrow n-1\) transition does not affect long-term dynamics.

We show the impact immobilization threshold, T, and the cascade factor, c, have on the immobilization process in Fig. 3. Within each frame, it can be seen that a higher immobilization threshold allows for longer freely diffusion periods, while across frames a higher cascade factor leads to longer immobilized periods. In Fig. 3(d)–(f), we simulated realizations of the processes (N(t), S(t)) for various combinations of T and c. The number of bound antibodies, N(t), is displayed by the purple trajectory, and the number of simultaneously bound Ab with a low immobilization threshold, S(t) when \(T = 1\), and with a higher immobilization threshold, S(t) when \(T = 10\), are shown by the green and blue trajectory, respectively. Moving left to right, the factor by which the Ab–mucin binding rate changes after immobilization increases, \(c = 1, 20,\) and 200, respectively. In Fig. 3(a)–(c), we show how these processes dictate the movement of the virion. The virion with process (N(t), S(t) when \( T=1)\) is colored in green, while (N(t), S(t) when \(T=10)\) is colored in blue. For both trajectories, immobilized periods, \(S(t) \ge T\), are colored in red.

When immobilization does not affect the Ab–mucin binding rate (Fig. 3(d)), the process S(t) rapidly crosses the immobilization threshold (dashed line) resulting in a virion transitioning between states faster than the experimental time step (as shown in Fig. 3(a)), for both \(T = 1\) and \(T = 10\). By increasing the cascade factor (as shown in Fig. 3(e)–(f)), S(t) remains above the immobilization threshold, for observable periods. In this case, the simulated virions in Fig. 3(b), (c) change states on the experimental timescale of 20 seconds and longer than 20 s, respectively.

2.4.3 Our Approximation for the Stationary Probability of Being Immobilized

We assume that the antibody–virion dynamics are slow compared to the antibody–mucin dynamics. To approximate a virion’s long-term probability of being immobilized, we use a product of two factors. The first is the steady-state distribution for the number of surface-bound Ab, N(t). Then, we compute the stationary distribution for the number of simultaneously bound Ab, S(t), conditioned on each value \(N(t) = n\) (where \(n \in \{0, \ldots N_*\}\)).

We introduce the notation b(x, n, p) for the binomial probability mass function. That is, if \(X \sim \text {Binom}(n,p)\), then \({\mathbb {P}}\{X=x\} = b(x,n,p)\). Our approximation to the stationary distribution of immobilization can be understood as an average over the transitions of the fast process S(t). Let \(\sigma \) denote the time a particle spends in the immobilized state, and \(\tau \) the time a particle spends in the freely diffusing state. Then our approximation takes the form

$$\begin{aligned} {\widetilde{\pi }}([A]_{\text {exo}}) = q \sum _{n = T}^{N_*} \frac{{\mathbb {E}}(\sigma ; T,c,n) }{{\mathbb {E}}(\sigma ; T,c,n) + {\mathbb {E}}(\tau ; T,n)} \, b\bigg (n; N_*, \frac{[A]_0 + [A]_{\text {exo}}}{k_{\text {d}}+ ([A]_0 + [A]_{\text {exo}})}\bigg )\nonumber \\ \end{aligned}$$
(11)

where

$$\begin{aligned} \begin{aligned} {\mathbb {E}}(\sigma ; T,c,n)&= \frac{1}{ T m_{\text {off}}} \frac{\sum _{s = T}^{n}b\big (s; n, \frac{c m_{\text {on}}[M]}{m_{\text {off}}+ c m_{\text {on}}[M]} \big ) }{b\big (T; n; \frac{c m_{\text {on}}[M] }{m_{\text {off}}+ c m_{\text {on}}[M]}\big ) }, \\ \text {and } {\mathbb {E}}(\tau ; T,c,n)&= \tfrac{1}{(n-T+1) m_{\text {on}}[M]} \frac{\sum _{s = 0}^{T-1} b\big (s; n; \frac{m_{\text {on}}[M]}{m_{\text {off}}+ m_{\text {on}}[M]}\big )}{b\big (T-1; n, \frac{m_{\text {on}}[M] }{m_{\text {off}}+ m_{\text {on}}[M]}\big )}. \end{aligned} \end{aligned}$$
(12)

The derivation of Eqs. 11 and 12 relies on Markov Chain Theory and Renewal Theory and can be found in Appendix B.

It follows from the law of total expectation and the timescale approximation, the expected time immobilized and expected time freely diffusing are, respectively:

$$\begin{aligned} \begin{aligned} {\mathbb {E}}(\sigma )&= \sum _{n = T}^{N_*} {\mathbb {E}}(\sigma ; T, c, n)\ b\Big (n; N_*, \tfrac{[A]}{k_{\text {d}}+ [A]}\Big ); \\ {\mathbb {E}}(\tau )&= \sum _{n = T}^{N_*} {\mathbb {E}}(\tau ; T, c, n)\ b\Big (n; N_*, \tfrac{[A]}{k_{\text {d}}+ [A]}\Big ). \end{aligned} \end{aligned}$$
(13)

We say that a parameter vector is in the Slow Switching Regime if, for all tested exogenous Ab concentrations, the average times spent in the immobilized and diffusing states are more than 20 s. To be precise, we define

$$\begin{aligned} \varTheta _{\text {slow}} := \big \{ \varvec{}\theta \, : \, {\mathbb {E}}(\sigma ;[A]_{\text {exo}}, \varvec{}{\theta })> 20 \text { and } {\mathbb {E}}(\tau ;[A]_{\text {exo}}, \varvec{}{\theta }) > 20 \text { for all } [A]_{\text {exo}}\in [0,1]\big \}.\!\!\!\nonumber \\ \end{aligned}$$
(14)

2.5 Switch-Point Detection

We develop an algorithm for detecting whether there is a single switch from diffusion to immobilization or immobilization to diffusion. The mathematical model presented in Sect. 2.4 (Eq. 7) assumes complete immobilization, but in fact immobilized virions exhibit spatial motion. Bernstein and Fricks (2016) account for this spatial motion by describing the bound state as a diffusing particle trapped in a potential well. Using an expectation–maximization algorithm, they provide an evolving probability for each particle that it is in an immobilized or diffusing state. In contrast to the many-switch paths considered by Bernstein and Fricks, we argue in Sect. 3.1.2 that the virion paths in our data set have at most one or two switches. We therefore developed and implemented a Bayesian algorithm that is designed to identify the presence of a single switch point.

To derive a likelihood function, we extend our SDE model Eq. 7 to include a path-specific trapping potential well, similar to Bernstein and Fricks (2016). Our extended model for a [diffusion \(\rightarrow \) immobilization] switch is

$$\begin{aligned} d{\mathbf {X}}(t)&= {\left\{ \begin{array}{ll} \sqrt{2D}d{\mathbf {W}}(t) &{} 0 \le t \le \tau \\ -{\tilde{\kappa }}({\mathbf {X}}(t) - {\mathbf {X}}(\tau )) dt + \sqrt{2D}d{\mathbf {W}}(t) &{} \tau < t \end{array}\right. } \nonumber \\ {\mathbf {X}}(0 )&= 0 \end{aligned}$$
(15)

and for [immobilization \(\rightarrow \) diffusion], we have

$$\begin{aligned} d{\mathbf {X}}(t)&= {\left\{ \begin{array}{ll} -{\tilde{\kappa }}({\mathbf {X}}(t) - {\mathbf {X}}(0)) dt + \sqrt{2D}d{\mathbf {W}}(t) &{} 0 \le t \le \tau \\ \sqrt{2D}d{\mathbf {W}}(t) &{} \tau < t \end{array}\right. } \nonumber \\ {\mathbf {X}}(0 )&= 0 \end{aligned}$$
(16)

where \({\mathbf {X}}(t) = (X(t), Y(t))^T\) and \({\mathbf {W}}(t)\) is 2d Brownian Motion. These SDEs are derived from the Langevin equation for particles diffusing in a quadratic (Hookean spring) potential well. The constant \({\tilde{\kappa }}= \kappa / \gamma \) where \(\kappa \) is the spring constant and \(\gamma \) is the viscous drag experienced by the particle. Due to the fluctuation–dissipation relationship, \(\gamma \) also appears in the diffusivity constant, which has the form \(D = k_B {\mathcal {T}} / \gamma \), where \(k_B\) is Boltzmann’s constant and \({\mathcal {T}}\) is the temperature of the fluid. To obtain an analytically trackable likelihood function, we introduce simplifying assumptions that (1) the switch occurs exactly at an observation time point, and (2) there is no measurement error. We derive the likelihood function in Appendix C.

We take a Bayesian approach to jointly estimate D, \({\tilde{\kappa }}\), and \(\tau \) under both switching scenarios using a Gibbs sampling algorithm. If the \(95\%\) credible region for \(\tau \) is completely contained within the interval [0.1\({T}_{\text {final}}\), 0.9\({T}_{\text {final}}\)] where \({T}_{\text {final}}\) is the duration of a path, then we say that path is a candidate for switching. For both switching scenarios, we estimated a false discovery rate for this criterion by simulating freely diffusing particles and setting the false discovery rate to the percent of simulated Brownian particles that were labeled as candidates for switching for the given switching model, as shown in Eqs. 16 or 15. Similarly, we estimated the power of criterion through simulation. For both scenarios, we simulated particles that switched states once and set the power to the fraction of paths that were candidates for switching. See Sect. 5 for more details on how these tests were constructed, and the results are presented in Sect. 3.1.2.

2.6 Uncertainty Quantification

The model given by Eq. 11 depends on the parameter vector \( \varvec{}{\theta } = (T, c, N_*, q, [A]_0, k_{\text {d}}, m_{\text {off}}, \alpha ).\) In specifying the model to HSV-IgG data (Sect. 2.1) we set \(k_{\text {d}}= 0.8969\) (McKinley et al. 2014) and \(\alpha = 0.90\) (Olmsted et al. 2001). The Ab–mucin binding and unbinding rates have not been directly estimated. We assume they are fast compared to the experimental timescale and, for example, set \(m_{\text {off}}= 100s^{-1}\). To assess the remaining parameters, \( \varvec{}{\theta } = (T, c, N_*, q, [A]_0)\)—which are the immobilization threshold value, the binding cascade factor, the number of sites on the surface of virions, the virion–Ab interaction probability, and the endogenous Ab concentration—we employed the numerical method of profile likelihoods (Eisenberg and Hayashi 2014; Raue et al. 2009). We used the numerically obtained relationships among parameters to obtain conditions on \(\varvec{}{\theta }\) such that the Switching Diffusion Hypothesis (in the Slow Switching Regime) is consistent with the data.

In order to quantify the model’s error in predicting the immobilized fraction, for each donor i, we partitioned the paths according to exogenous Ab concentration \(\{[A]_j\}_{j=1}^5\) and introduced the following residual function:

$$\begin{aligned} \chi ^2_i(\varvec{}\theta ) = \sum _{j=1}^5 N_{ij} \frac{({\widetilde{\pi }}([A]_j; \varvec{}{\theta }) - {\hat{p}}_{ij})^2}{{\widetilde{\pi }}([A]_j, \varvec{}{\theta })(1 -{\widetilde{\pi }}([A]_j; \varvec{}{\theta }))}, \end{aligned}$$
(17)

where \({\widetilde{\pi }}([A]_j; \varvec{}{\theta })\) denotes the model evaluated at \([A]_j\) with parameters \(\varvec{}{\theta }\) (as defined in Eq. 11), while \(N_{ij}\) and \({\hat{p}}_{ij}\) are, respectively, the number of paths observed and the fraction that are immobilized in the jth subpopulation associated with donor i. Assuming a normal approximation to the binomial distribution, our residual function can be seen as the sum of five independent squared normal random variables, i.e., with a \(\chi ^2\)-distribution with 5 degrees of freedom.

2.6.1 Numerical Method of Profile Likelihoods to Deduce Parameter Identifiability

Because we assume normal approximation to the binomial distribution, working with a residual function is equivalent to using a likelihood function to define confidence intervals (Press et al. 1996; Raue et al. 2009). For ease of notation in this section, we will suppress the dependence on i when considering the residual function \(\chi ^2(\varvec{}{\theta })\) for donor i.

To discuss identifiability of our model parameters, we use the nomenclature introduced by Raue et al. (2009). Our minimum residual estimator is defined to be \({\widehat{\theta }} := \text {argmin}[\chi ^2(\varvec{}\theta )]\). The likelihood-based confidence region of level\(\alpha \) for \(\varvec{}\theta \) is then defined to be

$$\begin{aligned} \varTheta _{\alpha ,df} := \{\varvec{}{\theta } \, : \, \chi ^2(\varvec{}\theta ) - \chi ^2({{{\widehat{\theta }}}}) < \chi ^2(\alpha , df) \}, \end{aligned}$$
(18)

where \(\chi ^2(\alpha , df)\) is the \(\alpha \) quantile of the \(\chi ^2\) distribution with df degrees of freedom. When establishing a confidence interval for one of the parameters, we set \(df = 1\). When establish a confidence region for multiple parameters, we set df equal to the number of parameters (Press et al. 1996).

A parameter \(\theta _k\) is said to be structurally identifiable when there is a unique minimum of \(\chi ^2(\varvec{}{\theta })\) with respect \(\theta _k\), i.e., if there exists a unique \(\theta _k\) such that

$$\begin{aligned} \theta _k = \big (\mathrm {argmin}_{\varvec{}\theta \in {\mathbb {R}}^5} \{\chi (\varvec{}\theta ))\}\big )_k. \end{aligned}$$

Alternatively, \(\theta _k\) can be unidentifiable due to the structure of the model or because the quality and quantity of the data are insufficient in estimating \(\theta _k.\) For the former case, we say \(\theta _k\) is structurally unidentifiable if the set

$$\begin{aligned} \theta _{\min } := \{ \varvec{}\theta \, : \, \chi (\varvec{}\theta ) = \min _{\vartheta \in {\mathbb {R}}} \chi (\vartheta ) \} \end{aligned}$$

is not unique and contains at least two elements whose \(\theta _k\) components are distinct. This often occurs when there is a functional relationship \(\phi \) among \(\theta _k\) and at least one other parameter, say \(\theta _j\) such that \(\chi \) can be expressed directly in terms of \(\phi (\theta _k,\theta _j)\). As for the latter, data-restricted type of unidentifiability, we say \(\theta _k\) is practically unidentifiable when a unique minimum exists of \(\chi ^2(\varvec{}{\theta })\) with respect to \(\theta _k\), but the likelihood-based confidence interval for \(\varvec{}{\theta }\) extends infinitely in increasing and/or decreasing values of \(\theta _k.\)

These definitions can be interpreted graphically using profile likelihoods. For residual function \(\chi ^2(\varvec{}{\theta })\), the profile likelihood of the k-th parameter defined to be

$$\begin{aligned} \chi _{\text {PL}}^2(\theta _k) = \min _{\theta _{j \ne k}}\big [ \chi ^2(\varvec{}\theta ) \big ]. \end{aligned}$$
(19)

If \(\theta _k\) is a structurally identifiable parameter, then \(\chi _{\text {PL}}^2(\theta _k)\) exceeds the threshold \(\varDelta _{\alpha }\) for both increasing and decreasing values of \(\theta _k\) forming a deep valley around \({{\widehat{\theta }}}_k\). If \(\theta _k\) is structurally unidentifiable, the profile likelihood is flat. Lastly, if \(\theta _k\) is practically unidentifiable, \(\chi _{\text {PL}}^2(\theta _k)\) obtains a unique minimum but does not exceed \(\varDelta _{\alpha }\) in increasing and/or decreasing values of \(\theta _k,\) forming a shallow valley around \({{\widehat{\theta }}}_k.\)

We further investigate unidentifiable combinations of parameters by extending Eq. 19 to profile parameter \(\theta _j\) and \(\theta _k\) simultaneously,

$$\begin{aligned} \chi ^2_{\text {PL}}(\theta _j, \theta _k ) : = \min _{\theta _{i \notin \{j, k\}}} \chi ^2(\varvec{}{\theta }). \end{aligned}$$
(20)

Structural relationships between the two profile parameters manifest as flat valleys extending infinitely along the functional relationship in the contour plots of \(\chi ^2_{\text {PL}}(\theta _j, \theta _k )\). We note this flat valley only traces out the functional relationship \(\theta _j\) and \(\theta _k\) when the dimension of the parameter space is larger than 2.

3 Results

3.1 Data Do Not Support the Incremental Knockdown Hypothesis for a 20 s Time Frame

3.1.1 No Evidence of Fast Switching: Ensemble Effective Diffusivities of the Free Subpopulation are the Same Regardless of Exogenous Ab Concentration

For each donor/Ab–concentration combination, the associated sample of virions contained a clear division among the tracked particles’ MSD and ACF behavior. We used the classification scheme described in Sect. 2.3 to label each tracked virion as Immobilized, Freely Diffusing, Subdiffusive or Outlier. The Immobilized class was characterized by low effective diffusivity (\(< 10^{-1}\,\upmu \text {m}^2/\text {s}\)) and either anti-persistent or uncorrelated increment processes. Meanwhile, the Freely Diffusing class had uncorrelated increment processes and effective diffusivities larger than \(0.2\,\upmu \text {m}^2/\text {s}\). The Subdiffusive and Outlier classifications were rare and did not appear in all samples. For this reason, we removed these categories from the analysis but give a description of them in the SI. In Fig. 4(a)–(c), we display the results of the classification for Donor F08 at 0, 0.1, and 1 \(\upmu \text {g}/ \text {mL}\) added anti-HSV IgG in terms of \(D_{\mathrm {eff}}\) and the average of the x- and y-ACF, as defined in Sects. 2.2.3 and 2.2.1, respectively. The clear separation of groups and locations of the clusters were qualitatively similar for the other donors (further figures included in supplementary materials).

Fig. 4
figure 4

a–c The unweighted composition of the tracked virions for Ab concentration 0, 0.1, and \(1.0\,\upmu \text {g}/ \text {mL}\), respectively, for Donor F08. Each point corresponds to a tracked virion with the given estimated diffusivity on a \(\log 10\) scale and average ACF value. The character of points denotes clusters prescribed by the hierarchical clustering algorithm, and color of the point denotes the class of the cluster. d–f The pathwise MSD for all the tracked virions for Donor F08 at \([A]_{\text {exo}}= 0, 0.1, \) and \(1.00\,\upmu \text {g}/ \text {mL}\). The colors, green, red, and blue, denote the final clusters, Freely Diffusing, Immobilized, and Subdiffusive, respectively. Reference line with slope \(=\) 1 is denoted in black. (We note that the relative size of the different classes in this figure is not reweighted by path length as it is in the population counts reported in Fig. 6.)

Full size image

The pathwise MSDs for Donor F08 virions are displayed in Fig. 4(d)–(f), and we note the similarity of the Freely Diffusing category of virions across all three panels. The Incremental Knockdown Hypothesis would predict that freely diffusing virions would be “slower and slower” in the presence of more and more Ab. However, we found that the diffusivities of the Freely Diffusing classes are consistent across all exogenous Ab concentrations. In Fig. 5, we display this fact in two ways. In the left panel, we display the ensemble MSD averaged over the Freely Diffusing (green triangles) and Immobilized (red x’s) populations for each Ab concentration. There is remarkable overlap within each group. Moreover, in the right panel, we display the ensemble effective diffusivity for the Freely Diffusing class at the various exogenous Ab concentrations for all donors. While there is variation in the effective diffusivity, the overlapping \(BC_a\) confidence intervals indicate there is insufficient evidence to conclude the effective diffusivity decreases with antibody concentration. (We provide \(95\%\) weighted bootstrap confidence intervals for each estimate in Fig. 13 of supplementary material )

Fig. 5
figure 5

a Ensemble MSD of the Freely Diffusing class and the immobilized class at various exogenous antibody concentrations represented by the green and red curves, respectively, for Donor F08. The black line refers to the ensemble MSD of Brownian particles, slope equal to 1. b The estimated ensemble effective diffusivity of the free population versus exogenous antibody concentration where the shade and point style of the curve correspond to Donor. See Fig. 13 of supplementary materials for the ensemble effective diffusivity with \(95\%\)\(BC_a\) confidence intervals

Full size image

We can express this finding in terms of a statistical test by comparing the weighted ensemble effective diffusivity for the freely diffusing subpopulation at the two extreme Ab concentrations. We used a one-tailed paired difference hypothesis test:

$$\begin{aligned} H_0: \langle D_{\mathrm {eff}}([A]_1) \rangle - \langle D_{\mathrm {eff}}([A]_5) \rangle = 0, \qquad H_A: \langle D_{\mathrm {eff}}([A]_1) \rangle - \langle D_{\mathrm {eff}}([A]_5) \rangle > 0\nonumber \\ \end{aligned}$$
(21)

for \([A]_1 = 0.0 \upmu \text {g}/ \text {mL}\) and \([A]_5 = 1.0 \upmu \text {g}/ \text {mL}.\) At an \(\alpha = 0.05\) level of significance, we failed to find significant evidence that the ensemble effective diffusivity of the freely diffusing population decreased when exogenous Ab concentration increased from zero exogenous Ab to the highest concentration \((t_{6} = 0.2567,\)p value\( = 0.4030 )\). We report the results of paired difference tests for all other combinations of the tested exogenous Ab concentration in Table 10 of supplementary materials.

3.1.2 Little Evidence of Switching on the Experimental Timescale

We found little evidence that virions switch between states on the experimental timescale of 20 s. If tracked particles were typically experiencing many subtle switches, we expect that their computed effective diffusivities would be diminished by a factor determined by the time spent immobilized. Moreover, because there are distinct behavioral regimes, the distribution of the increment processes is essentially a mixture of two Gaussian distributions (one for the Immobilized state and one for the Freely Diffusing state). This would manifest itself as a violation of linearity in qqnorm plots, which we do not see for the vast majority of HSV virion paths.

While the qqnorm test can identify paths that might experience switches, they do not affirm the presence of a switch. To this end, we developed a Bayesian method for identifying whether there is a single switch point in a given virion path, described in Sect. 2.5. We say a path of duration \({T}_{\text {final}}\) is a candidate for switching if the \(95\%\) credible region for \(\tau \) was completely contained within the interval [0.1\({T}_{\text {final}}\), 0.9\({T}_{\text {final}}\)]. The method was very effective on simulated data. When we applied the method to simulated Brownian motion (Freely Diffusing), we found a 0.0119 and 0.0080 False Discovery Rate of [diffusion \(\rightarrow \) immobilization] switches and [immobilization \(\rightarrow \) diffusion] switches, respectively. On the other hand, 96.38% of the simulated [diffusion \(\rightarrow \) immobilization] paths were correctly identified as [diffusion \(\rightarrow \) immobilization] switches, while 94.37% of the simulated [immobilization \(\rightarrow \) diffusion] paths were identified as [immobilization \(\rightarrow \) diffusion] switches (Table 2). Under this method, we found that 1.12\(\%\) of the Freely Diffusing class (1689 total tracked virions) were identified as [diffusion \(\rightarrow \) immobilization] switch candidates and 1.24\(\%\) of the free populations were [immobilization \(\rightarrow \) diffusion] switch candidates. We therefore concluded that state switches occurred relatively rarely on the experimental time scale.

Table 2 Fraction of freely diffusing virions that possibly switched states once by Donor
Full size table
Fig. 6
figure 6

a The weighted proportion of the 4 classes for Donor F08 at the various tested exogenous Ab concentrations. b Weighted proportions of immobilized virions for each donor. See Fig. 12 of supplementary materials for plots with \(95\%\)\(BC_a\) confidence intervals

Full size image

3.1.3 Fraction Immobilized Increases with Exogenous Antibody Concentration

While Ab concentration did not seem to affect the behavior of virions labeled Freely Diffusing, it did have a significant effect on the fraction of virions that were placed in this class. This is consistent with the findings reported in Wang et al. (2014). We computed the Immobilized fraction for each Donor/Ab–concentration sample using the method discussed in Sect. 2.3.1 and displayed the results in Fig. 6, where each curve in the panel (b) corresponds to a different donor. While there is heterogeneity in the fraction of Immobilized virions across donors, there is a visible overall increase in proportion immobilized from 0 to \(1 \ \upmu \text {g}/ \text {mL}\). This qualitative assessment is supported by statistical evidence provided by non-overlapping \(BC_a\) confidence intervals between the extreme exogenous Ab concentrations given in Fig. 12 of supplementary materials.

For each donor, the fraction of Immobilized virions increased with Ab concentration in the 0–0.333 \(\upmu \text {g}/ \text {mL}\) range and seemed to be saturated at higher Ab concentrations. We tested the significance of this observed trend by fitting a negative exponential growth model with predictors: exogenous antibody concentration and individual effect terms relative to Donor F08. Let \(\chi _k\) be the indicator function that a virion in the kth donor sample is in the immobilized state. Our negative exponential growth model takes the form

$$\begin{aligned} P(\chi _k = 1) = (\beta _0 + \beta _k) - e^{-(\alpha _0+ \alpha _{\text {exo}} [A]_{\text {exo}}) + \alpha _k} \end{aligned}$$
(22)

where \(\alpha _k\) and \(\beta _k\) are the effect terms for the k-th donor. We found the exogenous antibody concentration (\(\alpha _{\text {exo}} = 15.920\), p value\( < 0.001\)), the growth rate due to the baseline donor (\(\alpha _0 = -0.8427\), p value \(=\) 0.0043), and baseline saturation probability (\(\beta _0 = 0.9138\), p value\( <0.0001\)) were statistically significant in predicting the immobilization probability, whereas the constants accounting for deviations from the baseline due to donor sample were not significant. The model was fit using the R command nls() with the minimization algorithm set to Gauss–Newton’s method.

3.2 The Simple Linear Model Predicts Fast Switching

The results from Sect. 3.1 provide evidence against the hypothesis that switching between the diffusing and immobilized states is fast relative to the experimental timescale. Our next goal was to determine whether there is a parameter regime that predicts slow switching while simultaneously being consistent with the exogenous Ab-dependent Immobilization data displayed in Fig. 6. This analysis depends strongly on two assumptions: (1) whether one virion-bound Ab is sufficient to cross-link the virion to mucin and (2) whether Ab–mucin binding rates increase when the virion is immobilized, the so-called cascade effect. We introduced two variables—T, the threshold number, and c, the cascade factor—in our general model to account for these possible effects. In recent works, it has been assumed either that \(T=c=1\) (Chen et al. 2014; Newby et al. 2017) or that \(T=1\) and \(c>1\) (Wessler et al. 2015). We refer to \(T=c=1\) as the simple linear model (SLM) because all the CTMC transition rates are linear. By computing the expected durations of the immobilized and diffusing states (Eq. 13, derivation in Appendix B.2), we were able to show that the data are not consistent with the SLM, or any case where \(T=1\).

We say a model is consistent with the observed data for a specified donor if there exists a parameter vector \(\varvec{}{\theta }\) that is within the 95% confidence region for the Immobilized Fraction data (denoted \(\varTheta _{\alpha , df}\), defined in Eq. 18) and also predicts expected state times larger than 20 s (denoted \(\varTheta _{\text {slow}}\), defined in Eq. 14). In Fig. 7, we demonstrate that the SLM is not consistent with the data for Donor F08. In the left panel, we show a 2d profile likelihood plot for the endogenous Ab concentration \([A]_0\) and number of virion-surface-binding sites \(N_*\). For each \(([A]_0,N_*)\) pair, we calculated the best fit for the remaining parameter q, the virion–Ab interaction probability, and displayed the residual value by the shading (darker means better fits). The black region represents the 95% confidence region for these two parameters. We uniformly sampled this confidence region, \(\varTheta _{0.05, 3}\), and displayed the predicted Immobilized Fraction curves for these parameter samples in panel (b) and the Ab–concentration-dependent expected state durations in panel (c). We note that all parameter combinations in \(\varTheta _{0.05, 3}\) had diffusing states that lasted less than 0.1 s for all values of \([A]_{\text {exo}}\). We repeated this analysis for all donors and in each case found that \(\varTheta _{0.05, 3}\cap \varTheta _{\text {slow}}= \emptyset \).

Fig. 7
figure 7

a, d Profile likelihood contour plots (Donor F08) for \(\chi _{\text {PL}}^2([A]_0, N_*)\) and \(\chi _{\text {PL}}^2(c,N_*)\) when \(T=c=1\) and \(T=19\), respectively. Darker shades correspond to smaller profile likelihood values and the black region corresponds to the 95% confidence regions \(\varTheta _{0.05, 3}\) and \(\varTheta _{0.05, 5}\). b, e Predicted Immobilized Fraction curves (gray lines) for \(\theta \) sampled from \(\varTheta _{0.05, 3}\) and \(\varTheta _{0.05, 5}\). The black curve is the prediction of the best fit in each case for Donor F08. The observed Immobilized Fraction is shown by the purple line with triangles. c, f Expected duration of Immobilized (red curves) and Freely Diffusing (green curves) states for \(\theta \) sampled from \(\varTheta _{0.05, 3}\) and \(\varTheta _{0.05, 5}\). When \(T=c=1\), frame (c), none of predicted state times are above 20 s, horizontal black line. On the other hand, when \(T=19\), frame (f), there are some parameter combinations that do yield slow switching. These are marked in light blue as appropriate in Panels (d)–(f)

Full size image

3.3 Threshold and Binding Cascade Parameters Allow Slow Switching

By allowing the immobilization process to require multiple cross-linking antibodies, \(T > 1\), and for the Ab–mucin dynamics to be state dependent, \(c \ne 1\), we found both that (1) the subset of parameters that lead to slow switching is non-empty (\(\varTheta _{\text {slow}}\ne \emptyset \)) and (2) there is an overlap between slow-switching parameters and parameters that fit the Immobilized Fraction data well (\(\varTheta _{0.05, 5}\cap \varTheta _{\text {slow}}\ne \emptyset \)). For example, in Fig. 7 panels (d)–(f) we demonstrate this fact assuming \(T = 19\) for Donor F08. The 2d profile likelihood plot in panel (d) shows an inverse relationship between \(N_*\) and the cascade factor c. Again the black region corresponds to all \((c,N_*)\) pairs that appear in \(\varTheta _{0.05, 5}\). For a uniform sample of such pairs, in panel (e) we display the Immobilized Fraction predictions, and in panel (f) the corresponding expected immobilization and diffusion state durations. Only a small subset of \(\varTheta _{0.05, 5}\) allows for slow switches. We mark this subset in blue in all three panels. Notably, conditioned on \(T=19\), we have that \(N_*\le 120\), which is somewhat smaller than the typical estimate for \(N_*\). In the next section, we note that assuming higher values for T leads to higher allowable values for \(N_*\). This type of result holds for all donors: for sufficiently high assumed T, the corresponding parameters sets \(\varTheta _{0.05, 5}\) and \(\varTheta _{\text {slow}}\) overlap.

Fig. 8
figure 8

a Parameter combinations of T and c that predict expected immobilized times greater than 20s, red points, and predict expected freely diffusing times greater than 20s, green points, assuming \(N_*= 300\) and \([A]_0 = 0.1 \upmu \text {g}/ \text {mL}\). The overlapping (T, c) combinations (blue points) are those combinations that satisfy slow-switching condition and subset \((T, c_{\min }\)) are denoted by black. b The minimal value of T required for our model to predict slow switching as a function of \(N_*\), orange curve. Given an \(N_*\) and corresponding minimal T pair, the minimal value of c required for our model to predict slow switching is denoted by the brown curve. The endogenous Ab concentration is fixed at \([A]_0 = 0.1 \upmu \text {g}/ \text {mL}\)

Full size image

By testing over a range of \(\varvec{}{\theta }= (T,c, N_*, [A]_0)\), we uncovered some relationships among the components of the parameter vectors \(\varvec{}{\theta }\) that yield slow switching \(\varTheta _{\text {slow}}\). We first investigated the relationship between T and c by fixing \(N_*\) and \([A]_0\). Noting that \(\mathbb {E}_{\theta }(\tau )\) is independent of c and \(\mathbb {E}_{\theta }(\sigma )\) is an increasing function in c, we calculated the minimal c required to satisfy the slow-switching condition, labeling this value \(c_{\text {min}}\). Though we could not obtain an explicit relationship between T and \(c_{\text {min}}\), we found that virions with a large immobilization threshold T can only satisfy the slow-switching condition if there is a corresponding large cascade effect, large \(c_{\min }\). To visualize this, in Fig. 8(a) we display the parameter combinations of \((T, c, N_*= 300, [A]_0= 0.1)\) that yield \(\mathbb {E}_{\theta }(\tau )> 20\) (green) and \(\mathbb {E}_{\theta }(\sigma ) > 20\) (red) for all exogenous antibody concentrations between 0 and \(1 \upmu \text {g}/ \text {mL}\). The overlapping region (blue points) corresponds to \(\theta \in \varTheta _{\text {slow}}\), and the combinations of interest \((T,c_{\text {min}})\) are shown in black.

We draw the conclusion that if \(N_*= 300\), then T must be at least 34 and c must be at least 63. If we increase the assumption about \(N_*\) while keeping \([A]_0\) fixed, then we found that the minimal allowable T and c for slow switching increase and decrease, respectively. We demonstrate this relationship in Fig. 8(b). For \([A]_0= 0.1 \upmu \text {g}/ \text {mL}\), the orange (circles) curve corresponds to the minimal T value (left y-axis) for the given \(N_*\) (x-axis) required such that \(\varvec{}{\theta } \in \varTheta _{\text {slow}}\) where \([A]_0= 0.1 \upmu \text {g}/ \text {mL}\). The brown (triangles) curve denotes the minimal c value (right y axis) required for the given \(N_*\), minimal T, and \([A]_0= 0.1 \upmu \text {g}/ \text {mL}\) to result in expected state times longer than 20 s.

3.4 Model with Threshold and Binding Cascade Parameter is Unidentifiable

As implied by the results in the preceding section, we found that the introduction of \(T > 1\) and \(c \ne 1\) resulted in issues with identifiability. That is to say, it appears that the confidence region \(\varTheta _{0.05, 5}\) is infinite even when restricted to the subspace \(\varTheta _{0.05, 5}\cap \varTheta _{\text {slow}}\). We use the Immobilized Fraction data for Donor F08 to demonstrate this fact but provide information for each Donor in supplementary materials. Throughout this section, we will use the terminology defined in Sect. 2.6.

Over the full parameter space \(\varTheta \), the 1d profile likelihoods revealed that all three of the parameters T, c, and \(N_*\) are practically unidentifiable over the range we tested. The profile likelihoods are displayed in black in Fig. 9(a)–(c). When we profiled the parameters T, c, and \(N_*\) restricted to the Slow Switching Regime, \(\varTheta _{0.05, 5}\cap \varTheta _{\text {slow}}\), we found T is still practically unidentifiable over the range \(T \ge 19\), while c is practically unidentifiable a large range of positive values. The number of binding sites \(N_*\) does seem to be identifiable, with a deep valley centered around the unique minimum at approximately \(N_*= 120\). These profile likelihoods are represented in blue in Fig. 9(a)–(c). The dashed lines correspond to the 95% CI boundaries for each parameter. Since the blue curves are below the confidence interval, we can say that there exist parameter combinations in the Slow Switching Regime that reasonably fit the Immobilized Fraction data in Fig. 9(e).

Fig. 9
figure 9

a–c The 1d profile likelihoods for the parameters: immobilization threshold T, cascade factor c, and number of Ab binding sites on the virion \(N_*\), respectively, over all tested parameter combinations (black curves) and when restricted to the Slow Switching Regime (blue curves). The 95% CI for each parameter consists of those parameter values with profile likelihood values below the dashed line

Full size image

4 Discussion

We have developed mathematical models and statistical methods to analyze the behavior of HSV virions diffusing in CVM in the presence of various concentrations of cross-linking Ab. With a few exceptions, we found that particle paths can be partitioned into two basic categories: Freely Diffusing and Immobilized. While the fraction of Immobilized virions increases with Ab concentration, we found that the mobility of the Freely Diffusing class is not Ab–concentration dependent.

Because we expect all the individual bonds to be reversible, virions should switch between the Freely Diffusing and Immobilized states. Previously, it had been hypothesized that such switches are rapid with respect to the experimental timescale, but our analysis contradicts that assumption. This raises the question of whether or not it is possible for the basic kinetic model to produce “slow-switching” paths where switches occur on a timescale much larger than the experimental time window. We found that this is possible if the model allows for a lower bound on the number of Ab necessary to immobilize a virion and assuming a “cascade effect” in Ab–mucin binding that encourages entanglement.

Introducing these extra features leads to a fundamental issue with unidentifiability in the statistical analysis. We can make claims like “the minimum number of antibodies needed to immobilize a virion must be greater than 20 or so”, but we cannot be more specific. In order to do so, we would need to have access to time series that are much longer than what is currently experimentally feasible.

While we have shown that it is possible for reversible kinetics to be consistent with the path data, it might also be possible to explain the data with a model that assume all binding events are irreversible. Unfortunately the available data cannot distinguish between the two models. One possible resolution is to conduct experiments that explicitly control for the time between the introduction of Ab to the virion population and the observation of virion trajectories. Based on our model, in which we assume the immobilization process is reversible prior to the system reaching stationarity, switching should be more common when the number of antibodies bound to surface epitopes is low. Therefore, starting the tracking immediately enhances the probability of observing state switches before any long-lasting immobilization events occur.

On the other hand, observing virions at different time points long after Ab introduction will help determine whether or not the system reaches a stationary distribution. If so, there should be substantial information in analyzing how (or if) that stationary distribution depends on the Ab concentration, and the rate at which that stationary distribution is achieved.

From a biological point of view, this uncertainty about the true timescale of state switching prevents conclusive answers for the type of model presented by Chen et al. (2014). Namely, in order to characterize the percentage of virions that can pass through a mucosal layer within two hours, we must know whether state switching takes place on the order of minutes, days or longer. The analysis in this work simply served to eliminate the possibility of switching on a timescale of milliseconds or seconds, which was the standing assumption.

The statistical methods and mathematical model introduced here apply to a broad class of biological systems that are composed of distinct subpopulations. Our classification scheme based on path-by-path analysis detects subpopulation dynamics that can be masked when considering only overall ensemble behavior. Clustering and then analyzing subpopulation ensemble statistics provide insight on the way the proportion and dynamics of these subpopulation change in response to the environmental factors. The model proposed in Sect. 2.4 can be modified to describe the general scenario when nanoparticles work to entrap a diffusing pathogen by anchoring the pathogen to the surrounding environment.