Information content in data sets for a nucleated-polymerization model

Abstract

We illustrate the use of statistical tools (asymptotic theories of standard error quantification using appropriate statistical models, bootstrapping, and model comparison techniques) in addition to sensitivity analysis that may be employed to determine the information content in data sets. We do this in the context of recent models [S. Prigent, A. Ballesta, F. Charles, N. Lenuzza, P. Gabriel, L.M. Tine, H. Rezaei, and M. Doumic, An efficient kinetic model for assemblies of amyloid fibrils and its application to polyglutamine aggregation, PLoS ONE 7 (2012), e43273. doi:10.1371/journal.pone.0043273.] for nucleated polymerization in proteins, about which very little is known regarding the underlying mechanisms; thus, the methodology we develop here may be of great help to experimentalists. We conclude that the investigated data sets will support with reasonable levels of uncertainty only the estimation of the parameters related to the early steps of the aggregation process.

Keywords:

• inverse problems
• polyglutamine and aggregation modelling
• information content
• sensitivity
• uncertainty quantification

Mathematics Subject Classification:

• 65M32
• 62P10
• 64B10
• 49Q12

1. Introduction

As mathematical models become more complex with multiple states, increasing numbers of parameters need to be estimated using experimental data. Thus there is a need for critical analysis in model validation related to the reliability of parameter estimates obtained in model fitting. A recent concrete example involves previous HIV models [1,6] with 15 or more parameters to be estimated. In [8], using recently developed parameter selectivity tools [7] based on parameter sensitivity-based scores, it was shown that a number of these parameters could not be estimated with any degree of reliability. Moreover, we found that quantifiable uncertainty varies among patients depending upon the number of treatment interruptions (perturbations of therapy). This leads to a fundamental question: How much information with respect to model validation can be expected in a given data set or collection of data sets?

Here we illustrate the use of other tools: asymptotic theories (as the data sample size increases without bound) of standard error (SE) quantification using appropriate statistical models, bootstrapping (repeated sampling of synthetic data similar to the original data), and statistical (ANOVA type) model comparison techniques – these are discussed in some detail in the recent monograph [11] as well as in numerous statistical texts. Such techniques are employed in addition to sensitivity theory in order to determine the information content in data sets. We pursue this in the context of recent models [25] for nucleated polymerization in proteins.

After briefly outlining the biological context of amyloid formation, we describe the model in Section 2. In Section 3, we investigate the statistical model to be used with our noisy data. This is a necessary step in order to use the correct error model in our generalized least-squares (GLS) minimization. This also reveals information on our experimental observation process. Once we have found parameters which allow a reasonable fit, we determine the confidence we may have in our estimation procedures. We do this in Section 4, using both the condition number of the covariance matrix and a sensitivity analysis. This reveals a smaller number of parameters (than those estimated in [25]) to be reasonably sensitive to the data sets, whereas others do not really affect the quality of the fits to our data. To further support our sensitivity findings, we then apply a bootstrapping analysis in Section 5. We are lead to four main parameters and compare their resulting errors with the asymptotic confidence intervals of Section 4. Finally, in Section 6, we carry out model comparison tests [3,5,11] as used in [10], and these lead us to select three parameters out of the nine original ones estimated in [25] that can be reliably estimated.

1.1. Protein polymerization

It is now known that several neurodegenerative disorders, including Alzheimer's disease, Huntington's disease, and Prion diseases, for example, mad cow, are related to aggregations of proteins presenting an abnormal folding. These protein aggregates are called amyloids and have become a focus of modelling efforts in recent years [12,25,29–31]. One of the main challenges in this field is to understand the key aggregation mechanisms, both qualitatively and quantitatively. In order to test our methodology on a relatively simple case, we focus here on polyglutamine (PolyQ)-containing proteins. This was also the case study chosen to illustrate the fairly general ordinary differential equation (ODE)–partial differential equation (PDE) model proposed in [25]; the reason for our choice is that, as shown in [25], the polymerization mechanisms prove to be simpler for PolyQ aggregation than for other types of proteins ( e.g. PrP) [21]. To understand data sets from experiments carried out by Human Rezaei and his team at INRA (Virologie et Immunologie Moleculaires) [25], we adapt the general model to this context. The data sets (DS1–DS4) of interest to us here are depicted in Figure 1 and record the evolution of normalized total polymerized mass in time.

Figure 1. The data sets of interest from [9,25]. The total polymerized mass is measured by Thioflavin T (ThT), which is one of the most common experimental tools for in vitro protein polymerization [25,29]. (available in colour online)

In [25] and a subsequent effort in [9], the authors sought to investigate several questions including (i) understanding the key polymerization mechanisms, (ii) how to numerically approximate the model, and (iii) how to select parameters and calibrate the model. Here we briefly summarize results related to (ii) and focus primarily on (iii).

2. The model

2.1. Original ODE model

We briefly outline the model that is the same as that of [25]. Let $(V,V^{\ast },c_i)$ be the main variables of interest consisting, respectively, of concentrations of the normal proteins that we will call monomers (basic subunits that are repeated in a chainlike fashion), of the monomeric proteins presenting an abnormal configuration that we will call conformers, and of the i-polymers made of i aggregated abnormal proteins. The following comprise the fundamental dynamics modelled in [25] which the interested reader may consult for further details:

• Monomer–conformer exchange: $V\stackrel{k_I^{+}}{\underset{k_I^{-}}{\rightleftharpoons }}{V}^{\ast }$. This represents the spontaneous formation of an active form of the monomer, denoted here $V^{\ast }$, out of the initially present inert form denoted V. The inert form cannot react and form fibrils, whereas the active conformer may.

• Nucleation: $\underset{i_0}{\underbrace{V^{\ast }+V^{\ast }+\cdots +V^{\ast }}}\stackrel{k_{\mathrm{on}}^N}{\underset{k_{\mathrm{off}}^N}{\rightleftharpoons }}{c}_{i_0}$. This is the spontaneous formation of the smallest stable polymer, formed by the addition of a certain number $i_0$ of active conformers. This smallest stable polymer is called the nucleus.

• Polymerization by conformer addition: $c_i+V^{\ast }\stackrel{k_{\mathrm{on}}^i}{\underset{}{\rightleftharpoons }}{c}_{i+1}$. Once a nucleus is formed, its size grows progressively by addition of an active conformer.

Other reactions like fragmentation and coalescence are negligible for the case of PolyQ-containing proteins (see [25] for experimental justification).

The law of mass action in the deterministic framework (see [5,27] and the numerous references therein) translates $A+B\stackrel{k_I^{+}}{\underset{k_I^{-}}{\rightleftharpoons }}{A}^{\prime }+B^{\prime }$ into the ODE $\mathrm{d}\left[A\right]/\mathrm{d}t=-k^{+}\left[A\right]\left[B\right]+k^{-}\left[A{]}^{\prime }\right[B{]}^{\prime }$. Using these basic ideas we obtain the infinite system of ODEs studied in [25] (1) $\frac{\mathrm{d}V}{\mathrm{d}t}=-k_I^{+}V+k_I^{-}{V}^{\ast },$(1) (2) $\frac{\mathrm{d}{V}^{\ast }}{\mathrm{d}t}=k_I^{+}V-k_I^{-}{V}^{\ast }+i_0{k}_{\mathrm{off}}^N{c}_{i_0}-V^{\ast }\sum _{i\ge i_0}{k}_{\mathrm{on}}^i{c}_i,$(2) (3) $\frac{\mathrm{d}{c}_{i_0}}{\mathrm{d}t}=k_{\mathrm{on}}^N(V^{\ast }{)}^{i_0}-k_{\mathrm{off}}^N{c}_{i_0}-k_{\mathrm{on}}^{i_0}{c}_{i_0}{V}^{\ast },$(3) (4) $\frac{\mathrm{d}{c}_i}{\mathrm{d}t}=V^{\ast }(k_{\mathrm{on}}^{i-1}{c}_{i-1}-k_{\mathrm{on}}^i{c}_i),i=i_0+1,\dots$(4) with initial conditions $V\left(0\right)=c_0,V^{\ast }\left(0\right)=0,c_{i_0}\left(0\right)=c_i\left(0\right)=0$ and the mass balance equation $\frac{\mathrm{d}}{\mathrm{d}t}(V+V^{\ast }+\sum _{i=i_0}^{\mathrm{\infty }}i{c}_i)=0.$

The experiments of interest to us measure the a-dimensional total polymerized mass, or Madim which is given by $M\left(t\right)=\sum _{i\ge i_0}i{c}_i\left(t\right).$

2.2. An approximate PDE system and the associated forward problem

Amyloid formations are characterized by very long polymers (a fibril may contain up to ${10}^6$ monomer units). A PDE version of the standard model, where a continuous variable x approximates the discrete sizes i, is thus a reasonable approximation for large amyloid polymers. However, for small polymer sizes this resulting continuum approximation does not work very well. Thus, we take a ‘hybrid approach’ of retaining the ODE for smaller sizes and using the PDE for larger ones [9].

We define a small parameter $\epsilon =1/i_M$, and let $x_i=i\epsilon$ with $i_M\gg 1$ being the average polymer size defined by $i_M=\frac{\sum _{i\ge i_0}i{c}_i}{\sum c_i}.$

Then after definition of dimensionless quantities $c^{\epsilon }(t,x)=\sum c_i{\mathbb{1}}_{[x_i,x_{i+1}]},$ we may obtain a PDE to replace the infinite ODE system. Rigorous derivations of such continuous integro-PDE models may be found in [22] for coagulation-fragmentation equations, in [15] for the limit of the so-called Becker– Döring system (polymerization and depolymerization reactions are the only reactions considered) towards its continuous limit called the Lifshitz–Slyozov model, and in [19] for the growth-fragmentation ‘Prion Model’. A formal derivation for a full model, also including nucleation, is carried out in [25].

Let $N_0\in \mathbb{N}$. We then use the approximation (5) $\begin{array}{rl}\frac{\mathrm{d}V}{\mathrm{d}t}& =-k_I^{+}V+k_I^{-}{V}^{\ast },\\ \frac{\mathrm{d}{V}^{\ast }}{\mathrm{d}t}& =k_I^{+}V-k_I^{-}{V}^{\ast }+i_0{k}_{\mathrm{off}}^N{c}_{i_0}-V^{\ast }\sum _{i\ge i_0}{k}_{\mathrm{on}}^i{c}_i,\end{array}$(5) (6) $\frac{\mathrm{d}{c}_{i_0}}{\mathrm{d}t}=k_{\mathrm{on}}^N(V^{\ast }{)}^{i_0}-k_{\mathrm{off}}^N{c}_{i_0}-k_{\mathrm{on}}^{i_0}{c}_{i_0}{V}^{\ast },$(6) (7) $\frac{\mathrm{d}{c}_i}{\mathrm{d}t}=V^{\ast }(k_{\mathrm{on}}^{i-1}{c}_{i-1}-k_{\mathrm{on}}^i{c}_i),i\le N_0,$(7) (8) ${\mathrm{\partial }}_t{c}^{\epsilon }(x,t)=-V^{\ast }{\mathrm{\partial }}_x\left(k_{\mathrm{on}}\right(x\left)c^{\epsilon }\right(x,t\left)\right),x\ge N_0,$(8) with initial conditions $V\left(0\right)=c_0,V^{\ast }\left(0\right)=0,c_{i_0}\left(0\right)=c_i\left(0\right)=0,c^{\epsilon }(x,0)=0,$ and the boundary condition $c^{\epsilon }(x=N_0,t)=c_{N_0}\left(t\right).$ Here we have passed to the continuous representations for chain lengths larger than $i=N_0$.

Then an assumed mass balance equation becomes $\frac{\mathrm{d}}{\mathrm{d}t}\left(V+V^{\ast }+\sum _{i=i_0}^{N_0}i{c}_i+{\int }_{N_0}^{\mathrm{\infty }}x{c}^{\epsilon }\left(x\right)\mathrm{d}x\right)=0.$

In [9], we considered requirements for a good discretization scheme including the following: (i) it should conserve the a-dimensional total polymerized mass (Madim), (ii) it should be fast, and most importantly, (iii) it should be accurate.

To ensure the conservation of mass, we replace the ODE for $V^{\ast }$ by the mass conservation equation and obtain (9) $\begin{array}{rl}\frac{\mathrm{d}V}{\mathrm{d}t}& =-k_I^{+}V+k_I^{-}{V}^{\ast },\\ V^{\ast }& =c_0-V-\sum _{i=i_0}^{N_0}i{c}_i-{\int }_{N_0}^{\mathrm{\infty }}x{c}^{\epsilon }\mathrm{d}x,\\ \frac{\mathrm{d}{c}_{i_0}}{\mathrm{d}t}& =k_{\mathrm{on}}^N(V^{\ast }{)}^{i_0}-k_{\mathrm{off}}^N{c}_{i_0}-k_{\mathrm{on}}^{i_0}{c}_{i_0}{V}^{\ast },\\ \frac{\mathrm{d}{c}_i}{\mathrm{d}t}& =V^{\ast }(k_{\mathrm{on}}^{i-1}{c}_{i-1}-k_{\mathrm{on}}^i{c}_i),i\le N_0,\end{array}$(9) (10) ${\mathrm{\partial }}_t{c}^{\epsilon }(x,t)=-V^{\ast }{\mathrm{\partial }}_x\left(k_{\mathrm{on}}\right(x\left)c^{\epsilon }\right(x,t\left)\right),x\ge N_0,$(10) with initial and boundary conditions as before.

We developed methodology for forward solutions in [9]. We first point out that the desired spatial computational domain is very large as determined by the maximum size of observed polymers, with range up to ${10}^6$. The peak in the distribution is at the left side of the domain of interest; for larger polymer sizes, the distribution is almost linearly decreasing.

Based on these and other considerations discussed in [9], the PDE was approximated by the Finite Volume Method (see [23] for discussions of Upwind, Lax–Wendroff, and flux limiter methods) with an adaptive mesh, refined towards the smaller polymer sizes. Furthermore, we kept the ratio between the step size and the corresponding mesh element constant, that is, we used $\mathrm{\Delta }{x}_i/x_i=q<1$ so that $x_i=1/(1-q)x_{i-1}$. This mesh is quasi-linear in the sense of $\mathrm{\Delta }{x}_{i-1}/\mathrm{\Delta }{x}_i=1+O\left(q\right)$. The resulting Upwind and Lax–Wendroff schemes are then consistent on the progressive mesh [23]. For further details on these schemes including examples demonstrating convergence properties, the interested reader may consult [9].

3. The inverse problem

A major question in formulating the model for use in inverse problem scenarios consists of how to best parametrically represent the polymerization parameters $k_{\mathrm{on}}^i$ of Equation (3(3) $\frac{\mathrm{d}{c}_{i_0}}{\mathrm{d}t}=k_{\mathrm{on}}^N(V^{\ast }{)}^{i_0}-k_{\mathrm{off}}^N{c}_{i_0}-k_{\mathrm{on}}^{i_0}{c}_{i_0}{V}^{\ast },$(3) ) and the function $k_{\mathrm{on}}$ of Equation (4(4) $\frac{\mathrm{d}{c}_i}{\mathrm{d}t}=V^{\ast }(k_{\mathrm{on}}^{i-1}{c}_{i-1}-k_{\mathrm{on}}^i{c}_i),i=i_0+1,\dots$(4) ) for our application. We do this in a combined piecewise continuous formulation.

For this continuous polymerization function $k_{\mathrm{on}}\left(x\right)$, we use the piecewise linear representation (Figure 2) $\begin{array}{r}{k}_{\mathrm{on}}\left(x\right)=\left\{\begin{array}{ll}{k}_{\mathrm{on}}^{min}+x\frac{k_{\mathrm{on}}^{max}-k_{\mathrm{on}}^{min}}{x_1{i}_{max}-i_0},& x\le x_1{i}_{max},\\ k_{\mathrm{on}}^{max},& x_1{i}_{max}\le x\le x_2{i}_{max},\\ k_{\mathrm{on}}^{max}-x\frac{k_{\mathrm{on}}^{max}}{i_{max}(1-x_2)},& x_2{i}_{max}\le x\le i_{max},\\ 0,& x\ge i_{max}.\end{array}\right.\end{array}$ In the numerical approximations, we chose $i_0=2,N_0=500$. The discrete polymerization parameters $k_{\mathrm{on}}^i$, $i=i_0,\dots ,N_0$, are then obtained as $k_{\mathrm{on}}^i:=k_{\mathrm{on}}(x=i).$

Figure 2. Parametric representation for $k_{\mathrm{on}}$.

Following [25], we chose to approximate $k_{\mathrm{on}}$ by a function as depicted in Figure 2. (According to our discussions with S. Prigent, H. Rezaei and J. Torrent, other choices like a Gaussian bell curve are also possible, but as we will subsequently conclude, the presently available data will not support estimation of parameters in these representations.) Thus with this parametrization we have five more parameters $k_{\mathrm{on}}^{min},k_{\mathrm{on}}^{max},$ the fractions $x_1,x_2,$ and $i_{max}$ in addition to the four basic parameters $k_I^{+},k_I^{-},k_{\mathrm{on}}^N,k_{\mathrm{off}}^N$ to be estimated using our data sets.

Thus, we seek to estimate (with acceptable quantification of uncertainties) the nine parameters $k_I^{+},k_I^{-},k_{\mathrm{off}}^N,k_{\mathrm{on}}^N$, and $k_{\mathrm{on}}\left(X\right)$ (represented in parametrical form depicted above with the five additional unknowns $k_{\mathrm{on}}^{min},k_{\mathrm{on}}^{max},$ $x_1,x_2,i_{max}$) that fit the data best. To do this we need an efficient discretization method as discussed above for the forward problem as well as a correct assumption on the measurement errors in the inverse problem.

3.1. Estimation of parameters

We make some standard statistical assumptions (see [5,11,17,28]) underlying our inverse problem formulations.

• Assume that there exists a true or nominal set of parameter ${\theta }_0=(k_I^{-},\dots ,i_{max}).$

• Let ${\mathcal{E}}_i$ be i.i.d. with $\mathbb{E}\left({\mathcal{E}}_i\right)=0$ and $\mathrm{cov}({\mathcal{E}}_i,{\mathcal{E}}_i)=\sigma {}_0^2$, where $i=1,\dots ,n$ and n is the number of observations or data points in the given data set.

Denote as $\hat{\theta }$ the estimated parameter for ${\theta }_0$. The inverse problem is based on statistical assumptions on the observation error in the data.

If we assume an absolute error data model, then data points are taken with equal importance. This is represented by observations (here ${ϵ}_i$ is a realization of the error process ${\mathcal{E}}_i$) (11) $y_i=M(t_i,{\theta }_0)+{ϵ}_i.$(11)

On the other hand, if one assumes some type of relative error data model, then the error is proportional in some sense to the measured polymerized mass. This can be represented by observations of the form (12) $y_i=M(t_i,{\theta }_0)+M(t_i,{\theta }_0{)}^{\gamma }{ϵ}_i,\gamma \in (0,1].$(12)

Absolute model error formulations dictate we use ordinary least-squares (OLS) inverse problem [5,11] given by (13) $\hat{\theta }=\arg \underset{\theta }{min}\sum _{i=1}^n(y_i-M(t_i,\theta ){)}^2,$(13) while for relative error models one should use inverse problem formulations with GLS cost functional (14) $\hat{\theta }=\arg \underset{\theta }{min}\sum _{i=1}^n{\left(\frac{y_i-M(t_i,\theta )}{M(t_i,\theta {)}^{\gamma }}\right)}^2,\gamma \in (0,1].$(14)

3.1.1. The residual plots.

To obtain a best statistical model, we used residual plots (see [5,11] for more details) with residuals given by $r_i=\frac{y_i-M(t_i,\hat{\theta })}{M(t_i,\hat{\theta }{)}^{\gamma }},\gamma \in [0,1].$

To illustrate what we are seeking for our data sets, we first used simulated relative error data (simulated data for $\gamma =1$), then carried out the inverse problems for both a relative error cost functional (i.e. $\gamma =1$) and an OLS cost functional (i.e. $\gamma =0$). We then plotted the corresponding residuals vs. time and also residuals vs. the model values. The first plots are related to the correctness of our assumption of independency and identical distributions i.i.d. for the data, whereas the second plots contain information as to the correctness of the form of our proposed statistical model (Figures 3 and 4).

Figure 3. Plots with simulated data: (a) correct cost function vs. time $(\gamma =1)$; (b) incorrect cost function vs. time $(\gamma =0)$.

Figure 4. Plots with simulated data: (a) correct cost function vs. model $(\gamma =1)$; (b) incorrect cost function vs. model $(\gamma =0)$.

3.2. Statistical models of noise

We next carried out similar inverse problems with data set (DS) 4 of our experimental data collection. We first used DS4 on the interval $t\in [0,8]$. Based on some earlier calculations we also chose the nucleation index $i_0=2$ for all our subsequent calculations. The residual plots given in Figures 5 and 6 suggest strongly that neither of the first attempts of assumed statistical models and corresponding cost functionals (absolute error and OLS or relative error with $\gamma =1$ and simple GLS) are correct.

Figure 5. (a) $M\left(t_k\right)$ (Madim) with OLS; (b) residuals vs. model: OLS.

Figure 6. (a) $M\left(t_k\right)$ with GLS, $\gamma =1$; (b) residuals vs. model: GLS.

Based on these initial results and the speculation that early periods of the polymerization process may be somewhat stochastic in nature, we chose to subsequently use all the data sets on the intervals $[t_0,8]$ where $t_0$ is the first time when $M\left(t_0\right)>0.12$ (thus 12% of the a-dimensional total polymerized mass). Moreover, we decided to use other values of γ between 0 and 1 to test DS4. Setting $i_0=2$, we focused on the question of the most appropriate values of γ to use in a GLS approach (again see [11] for further motivation and details). We then obtained the results with DS4 depicted in Figure 7. Analysis of these residuals suggest that either $\gamma =0.6$ or $\gamma =0.7$ might be satisfactory for use in a GLS setting.

Figure 7. Residuals for DS4 using different values of γ.

Motivated by these results, we next investigated the inverse problems for each of the four experimental data sets with initial concentration $c_0=200\text{μ}\mathrm{mol}$ and $i_0=2$. We carried out the optimization over all data points with $M\left(t_k\right)\ge 0.12$ and used the GLS method with $\gamma =0.6$. The resulting graphics depicted in Figure 8 again suggest that $\gamma =0.6$ is a reasonable value to use in our subsequent analysis of the PolyQ data with regard to its information content for inverse problem estimation and parameter uncertainty quantification.

Figure 8. Residuals for the four experimental data sets using $\gamma =0.6$.

4. SEs and asymptotic analysis

4.1. SEs for parameters using GLS

We employed first the asymptotic theory (as $n\to \mathrm{\infty }$) for parameter uncertainty summarized in [5,11,17] and references therein. In the case of GLS, the associated SEs for the estimated parameters $\hat{\theta }=(k_I^{+},\dots ,i_{max})$ (vector length ${\kappa }_{\theta }=9$) are given by the following construction (for details see Chaps 3.2.5 and 3.2.6 of [11]).

We may define the SEs by the formula ${\mathrm{SE}}_k=\sqrt{{\mathrm{\Sigma }}_{kk}\left(\hat{\theta }\right)},k=1,\dots ,9,$ where the covariance matrix Σ is given by $\mathrm{\Sigma }\left(\hat{\theta }\right)={\hat{\sigma }}^2\left({\chi }^{\mathrm{T}}\right(\hat{\theta }\left)W\right(\hat{\theta }\left)\chi \right(\hat{\theta }){)}^{-1}.$ Here $\chi =\frac{\mathrm{\partial }M}{\mathrm{\partial }\theta }=\left(\frac{\mathrm{\partial }M(t_1;\hat{\theta })}{\mathrm{\partial }\theta },\dots ,\frac{\mathrm{\partial }M(t_n;\hat{\theta })}{\mathrm{\partial }\theta }\right)$ is the sensitivity matrix of size $n\times {\kappa }_{\theta }=n\times 9$ (n being the number of data points and ${\kappa }_{\theta }=9$ being the number of estimated parameters) and W is defined by $W^{-1}\left(\hat{\theta }\right)=\mathrm{diag}\left(M\right(t_1;\hat{\theta }{)}^{2\gamma },\dots ,M(t_n;\hat{\theta }{)}^{2\gamma }).$ We use the approximation of the variance $\sigma {}_0^2\approx \hat{\sigma }(\hat{\theta }{)}^2=\frac{1}{n-{\kappa }_{\theta }}\sum _{i=1}^n\frac{1}{M(t_i;\hat{\theta }{)}^{2\gamma }}(M(t_i,\hat{\theta })-y_i{)}^2.$

To obtain a finite SE using asymptotic theory, the ${\kappa }_{\theta }\times {\kappa }_{\theta }=9\times 9$ matrix $F={\chi }^{\mathrm{T}}\left(\hat{\theta }\right)W\left(\hat{\theta }\right)\chi \left(\hat{\theta }\right)$ thus must be invertible. In the above problem, we do indeed obtain a good fit of the curve and good residuals (for the sake of brevity, not depicted here). However, we also found that the condition number of the matrix $F={\chi }^{\mathrm{T}}\left(\hat{\theta }\right)W\left(\hat{\theta }\right)\chi \left(\hat{\theta }\right)$ is $\kappa ={10}^{24}$. Looking more closely at the matrix F reveals a near linear dependence between certain rows, hence the large condition number. We thus quickly reach the following conclusions:

1. We obtain a set of parameters for which the model fits well, but we cannot have any reasonable confidence in them using the asymptotic theories from statistics ( e.g. see the references given above).

2. We suspect that it may not be possible to obtain sufficient information from our data set curves to estimate all nine parameters with a high degree of confidence. This is based on our calculations with the corresponding Fisher matrices as well our prior knowledge in that the graphs depicted in Figure 1 are very similar to Logistic or Gompertz curves. These curves can be quite well fit with parameterized models with only two or three carefully chosen parameters.

To assist in initial understanding of these issues, we consider components of the associated sensitivity matrices $\chi =\mathrm{\partial }M/\mathrm{\partial }\theta$.

4.2. Sensitivity analysis

For the sensitivity analysis, we follow [5,11] and carry out all computations using the differential system of sensitivity equations as detailed in those references. Hereafter all our analysis will be carried out using DS4 and the best estimate $\hat{\theta }$ obtained for the latter. We find that the model is sensitive mainly to four parameters: $k_I^{+},k_I^{-},k_{\mathrm{on}}^N,k_{\mathrm{off}}^N$. The sensitivities for the remaining parameters are on an order of magnitude of ${10}^{-6}$ or less. It also shows some sensitivity with respect to $x_1$. However, the parameter $x_1$ appears in the model only as the factor $x_1{i}_{max}$. The sensitivities depicted in the following use $\hat{\theta }$ for the nine best fit GLS parameters, that is, $\hat{\theta }$ for ${\kappa }_{\theta }=9$. We note that since we use the a-dimensional quantity M in the cost functionals (OLS or GLS), it is the sensitivity of this quantity with respect to the parameters θ, rather than any relative sensitivities, that will determine changes in the cost functionals to be minimized with respect to changes in the parameters (Figures 9–13).

Figure 9. (a) Sensitivity w.r.t. $k_I^{-}$; (b) sensitivity w.r.t. $k_I^{+}$.

Figure 10. (a) Sensitivity w.r.t. $k_{\mathrm{on}}^N$; (b) sensitivity w.r.t. $k_{\mathrm{off}}^N$.

Figure 11. (a) Sensitivity w.r.t. $k_{\mathrm{on}}^{min}$; (b) sensitivity w.r.t. $k_{\mathrm{off}}^{max}$.

Figure 12. (a) Sensitivity w.r.t. $x_1$; (b) sensitivity w.r.t. $x_2$.

Figure 13. (a) Sensitivity w.r.t. $i_{max}$; (b) sensitivity w.r.t. $x_{11}=i_{max}{x}_1$.

5. Sensitivity motivated inverse problems

Based on the sensitivity findings depicted above, we investigated a series of inverse problems in which we attempted to estimate an increasing number of parameters beginning first with the fundamental parameters $k_I^{+}$ and $k_I^{-}$. In each of these inverse problems, we attempted to ascertain uncertainty bounds for the estimated parameters using both the asymptotic theory described above and a GLS version of bootstrapping [13,14,16,18,20].

A quick outline of the appropriate bootstrapping algorithm is given next.

5.1. Bootstrapping algorithm: nonconstant variance data

We suppose now that we are given experimental data $(t_1,y_1),\dots ,(t_n,y_n)$ from the underlying observation process (15) $Y_i=M(t_i;{\theta }_0)+M(t_i;{\theta }_0{)}^{\gamma }\widetilde{{\mathcal{E}}_i},$(15) where $i=1,\dots ,n$ and the ${\tilde{\mathcal{E}}}_i$ are i.i.d. with mean zero and constant variance ${\sigma }_0^2$. Then we see that $\mathbb{E}\left(Y_i\right)=M(t_i;{\theta }_0)$ and $\mathrm{Var}\left(Y_i\right)={\sigma }_0^2{M}^{2\gamma }(t_i,{\theta }_0)$ , with associated corresponding realizations of $Y_i$ given by $y_i=M(t_i;{\theta }_0)+M(t_i;{\theta }_0{)}^{\gamma }\widetilde{{ϵ}_i}.$

A standard algorithm can be used to compute the corresponding bootstrapping estimate ${\hat{\theta }}_{\mathrm{boot}}$ of ${\theta }_0$ and its empirical distribution. We treat the general case for nonlinear dependence of the model output on the parameters θ. The algorithm is given as follows.

1. First obtain the estimate ${\hat{\theta }}^0$ from the entire sample $\left\{y_i\right\}$ using the GLS given in Equation (14(14) $\hat{\theta }=\arg \underset{\theta }{min}\sum _{i=1}^n{\left(\frac{y_i-M(t_i,\theta )}{M(t_i,\theta {)}^{\gamma }}\right)}^2,\gamma \in (0,1].$(14) ) with $\gamma =1$. An estimate ${\hat{\theta }}_{\mathrm{boot}}$ can be solved for iteratively as follows.

2. Define the nonconstant variance standardized residuals ${\overline{s}}_i=\frac{y_i-M\left(t_{;}{\hat{\theta }}^0\right)}{M(t_i;{\hat{\theta }}^0{)}^{\gamma }},i=1,2,\dots ,n.$

   Set $m=0$.

3. Create a bootstrapping sample of size n using random sampling with replacement from the data (realizations) {${\overline{s}}_1$,…,${\overline{s}}_n$} to form a bootstrapping sample $\{s_1^m,\dots ,s_n^m\}$.

4. Create bootstrapping sample points $y_i^m=M(t_i;{\hat{\theta }}^0)+M(t_i;{\hat{\theta }}^0{)}^{\gamma }{s}_i^m,$ where $i=1$,…,n.

5. Obtain a new estimate ${\hat{\theta }}^{m+1}$ from the bootstrapping sample $\left\{y_i^m\right\}$ using GLS.

6. Set $m=m+1$ and repeat steps 3–5 until $m\ge M$ where M is large (e.g. $M=1000$ ).

We then calculate the mean, SE, and confidence intervals using the formulae (16) $\begin{array}{rl}{\hat{\theta }}_{\mathrm{boot}}& =\frac{1}{M}\sum _{m=1}^M{\hat{\theta }}^m,\\ \mathrm{Var}\left({\theta }_{\mathrm{boot}}\right)& =\frac{1}{M-1}\sum _{m=1}^M({\hat{\theta }}^m-{\hat{\theta }}_{\mathrm{boot}}{)}^T({\hat{\theta }}^m-{\hat{\theta }}_{\mathrm{boot}}),\\ {\mathrm{SE}}_k\left({\hat{\theta }}_{\mathrm{boot}}\right)& =\sqrt{\mathrm{Var}({\theta }_{\mathrm{boot}}{)}_{kk}},\end{array}$(16) where ${\theta }_{\mathrm{boot}}$ denotes the bootstrapping estimator.

5.2. Estimation of two parameters

We first carried out estimation for the two parameters $k_I^{+}$ and $k_I^{-}$. We use the GLS formulation with $\gamma =0.6$. We fix globally (based on previous estimations with DS4) the parameter values

Table

$k_{\mathrm{on}}^N$	$k_{\mathrm{off}}^N$	$k_{\mathrm{on}}^{\min }$	$k_{\mathrm{on}}^{\max }$	x1	x2	imax
4616.962	93.332	1684.381	1.5152·10^9	0.0626	0.859	3.542×10^5

and used the initial guesses for the parameters given by

Table

	$k_I^{+}$	$k_I^{-}$
q0	2.1600	10.9270

We then used the bootstrapping algorithm and obtained the following means and SEs for $M=1000$ which, as reported in the following, compare quite well with the asymptotic theory estimates. The corresponding distributions are shown in Figures 14 and 15.

Figure 14. Two parameters estimation ($k_I^{+}$, $k_I^{-}$). Bootstrapping distribution for $k_I^{+}$. We use GLS and $M=1000$ runs.

Figure 15. Two parameters estimation ($k_I^{+}$, $k_I^{-}$). Bootstrapping distribution for $k_I^{-}$. We use GLS and $M=1000$ runs.

Table

	$k_I^{+}$(boot) (GLS)	$k_I^{-}$(boot) (GLS)	$k_I^{+}$(asymp) (GLS)	$k_I^{-}$(asymp) (GLS)
Mean	2.158	10.911	2.157	10.911
SE	0.0044	0.0247	0.00396	0.0225

5.3. GLS estimation of three parameters

We tried next to estimate three parameters using the GLS formulation with $\gamma =0.6$. Once again we fixed all the parameters describing the domain and the polymerization function $k_{\mathrm{on}}$ and also fixed either $k_{\mathrm{off}}^N$ or $k_{\mathrm{on}}^N$ in the corresponding inverse problems.

5.4. GLS estimation for $k_I^{+},k_I^{-},$ and $k_{\mathrm{on}}^N$

We fixed values as follows:

Table

$k_{\mathrm{off}}^N$	$k_{\mathrm{on}}^{\min }$	$k_{\mathrm{on}}^{\max }$	x1	x2	imax
93.33	1684.38	1.5×10^9	0.062	0.859	3.5×10^5

We used as initial parameter values:

Table

	$k_I^{+}$	$k_I^{-}$	$k_{\mathrm{on}}^N$
q0	2.1600	10.9270	4616.962

We obtained the estimated parameters together with the corresponding SEs, variances, and the condition numbers κ of the corresponding sensitivity matrices for the four data sets as reported in the following table. The $95\mathrm{\%}$ confidence results based on the asymptotic theory are also depicted for DS4 in Figure 16.

Table

	$k_I^{+}$	$k_I^{-}$	$k_{\mathrm{on}}^N$	SE	σ^2	κ
DS 1	2.26	13.49	4616.96	(0.012, 0.099, 53.925)	8.52×10^−6	8.89×10^10
DS 2	2.99	16.20	4616.96	(0.021, 0.151, 56.691)	9.67×10^−6	4.37×10^10
DS 3	2.18	15.76	9840.31	(0.011, 0.103, 90.466)	6.45×10^−6	3.94×10^11
DS 4	2.16	10.91	4616.96	(0.0089, 0.0649, 45.262)	6.36×10^−6	7.14×10^10

To compare these asymptotic results with bootstrapping, we carried out bootstrapping with DS4 for the estimation of $k_I^{+}$, $k_I^{-}$, and $k_{\mathrm{on}}^N$ with the same initial values as above. We then obtained the following means and SE for a run with $M=1000$, in comparison to the asymptotic theory.

Figure 16. Confidence intervals.

Table

	$k_I^{+}$(boot)	$k_I^{-}$(boot)	$k_{\mathrm{on}}^N$(boot)	$k_I^{+}$(asymp)	$k_I^{-}$(asymp)	$k_{\mathrm{on}}^N$(asymp)
Mean	2.153	10.887	4616.962	2.157	10.910	4616.962
SE	0.0039	0.0219	0.00003	0.0089	0.0649	45.262

Of particular interest are the values obtained for $k_{\mathrm{on}}^N$ and the bootstrapping SEs for $k_{\mathrm{on}}^N$ which are extremely small. It should be noted that the sensitivity of the model output on $k_{\mathrm{on}}^N$ is also very small. Thus one might conjecture that the iterations in the bootstrapping algorithm do not change the values of $k_{\mathrm{on}}^N$ very much and hence one observes the extremely small SE that are produced for the bootstrapping estimates (Figures 17–19).

Figure 17. Estimation for $k_I^{+}$, $k_I^{-},$ and $k_{\mathrm{on}}^N$: bootstrapping distribution for $k_I^{-}$ for GLS and $M=1000$ runs.

Figure 18. Estimation for $k_I^{+}$, $k_I^{-}$, and $k_{\mathrm{on}}^N$: bootstrapping distribution for $k_I^{+}$ for GLS and $M=1000$ runs.

Figure 19. Estimation for $k_I^{+}$, $k_I^{-}$, and $k_{\mathrm{on}}^N$: bootstrapping distribution for $k_{\mathrm{on}}^N$ for GLS and $M=1000$ runs.

5.5. GLS estimation for $k_I^{+},k_I^{-},$ and $k_{\mathrm{off}}^N$

In another test, we fixed $k_{\mathrm{on}}^N$ and instead estimate $k_{\mathrm{off}}^N$ (along with $k_I^{+}$ and $k_I^{-}$). We use the fixed values:

Table

$k_{\mathrm{on}}^N$	$k_{\mathrm{on}}^{\min }$	$k_{\mathrm{on}}^{\max }$	x1	x2	imax
4616.962	1684.381	1.5152×10^9	0.0626	0.859	3.542×10^5

and the initial guesses for the parameters to be estimated given by

Table

	$k_I^{+}$	$k_I^{-}$	$k_{\mathrm{off}}^N$
q0	2.1600	10.9270	108.256

We obtained the estimated parameters and corresponding SE.

Table

	$k_I^{+}$	$k_I^{-}$	$k_{\mathrm{off}}^N$	SE	σ^2	κ
DS 1	2.203	12.997	99.861	(0.011, 0.091, 1.208)	8.165×10^−6	4.912×10^7
DS 2	2.893	15.474	100.019	(0.019, 0.137, 1.279)	9.323×10^−6	2.486×10^7
DS 3	2.168	15.631	41.935	(0.011, 0.102, 0.424)	6.435×10^−6	9.125×10^6
DS 4	2.181	11.090	90.536	(0.009, 0.066 , 0.936)	6.289×10^−6	3.043×10^7

Also in this case, we carried out bootstrapping for DS4. The bootstrapping distributions for $k_I^{+}$, $k_I^{-}$, and $k_{\mathrm{off}}^N$ are found in Figures 20–22. We then obtained the following means and SEs for a run with $M=1000$ in comparison to the asymptotic theory.

Figure 20. Three parameters estimation ($k_I^{+}$, $k_I^{-},$ and $k_{\mathrm{off}}^N$): bootstrapping distribution for $k_I^{+}$. We used GLS and $M=1000$ runs.

Figure 21. Three parameters estimation ($k_I^{+}$, $k_I^{-},$ and $k_{\mathrm{off}}^N$): bootstrapping distribution for $k_I^{-}$. We used GLS and $M=1000$ runs.

Figure 22. Three parameters estimation ($k_I^{+}$, $k_I^{-}$ and $k_{\mathrm{off}}^N$): bootstrapping distribution for $k_{\mathrm{off}}^N$. We used GLS and $M=1000$ runs.

Table

	$k_I^{+}$(boot)	$k_I^{-}$(boot)	$k_{\mathrm{off}}^N$(boot)	$k_I^{+}$(asymp)	$k_I^{-}$(asymp)	$k_{\mathrm{off}}^N$(asymp)
Mean	2.169	11.013	91.254	2.181	11.090	90.536
SE	0.0094	0.0699	1.0392	0.009	0.066	0.936

5.6. Estimation of four main parameters

Following the sensitivity analysis detailed above, we tried to estimate a combination of the parameters $k_I^{+},k_I^{-},k_{\mathrm{on}}^N,k_{\mathrm{off}}^N$ for the parameter set with ${\kappa }_{\theta }=4$.

Parameters were fixed as follows from the original nine parameter fit:

Table

$k_{\mathrm{on}}^{\min }$	$k_{\mathrm{on}}^{\max }$	x1	x2	imax
1684	1.5×10^9	0.062	0.859	3.5×10^5

We obtained the following result for the estimation of the four parameters using the DS1–DS4. In all of them, the condition number of Fischer's information matrix κ is too large to invert. This along with the sensitivity results above strongly suggests that the data sets do not contain sufficient information to estimate four or more parameters with any degree of certainty attached to the estimates.

Table

	$k_I^{+}$	$k_I^{+}$	$k_{\mathrm{on}}^N$	$k_{\mathrm{off}}^N$	σ^2	κ
DS 1	2.1431	12.4751	4616.962	108.259	8.7219×10^−6	6.1226×10^19
DS 2	2.7995	14.7630	4616.957	108.4308	9.8694×10^−6	1.4442×10^19
DS 3	2.180	15.757	4618.599	41.369	6.4622×10^−6	1.881×10^17
DS 4	2.161	10.9278	4617.3316	93.3265	6.374×10^−6	2.144×10^18

6. Model comparison tests

A type of Residuals Sum of Squares (RSS)-based model selection criterion [3,5,11] can be used as a tool for model comparison for certain classes of models. In particular, this is true for models such as those given in [10] in which potentially extraneous mechanisms can be eliminated from the model by a simple restriction on the underlying parameter space while the form of the mathematical model remains unchanged. In other words, this methodology can be used to compare two nested mathematical models where the parameter set ${\mathrm{\Omega }}_{\theta }^H$ (this notation will be defined explicitly in Section 6.1) for the restricted model can be identified as a linearly restricted subset of the admissible parameter set ${\mathrm{\Omega }}_{\theta }$ of the unrestricted model. Indeed, the RSS-based model selection criterion is a useful tool to determine whether or not certain terms in the mathematical models are important in describing the given experimental data.

6.1. Ordinary least squares

We now turn to the statistical model (11(11) $y_i=M(t_i,{\theta }_0)+{ϵ}_i.$(11) ), where the measurement errors are assumed to be independent and identically distributed with zero mean and constant variance ${\sigma }_0^2$. In addition, we assume that there exists ${\theta }_0$ such that the statistical model (17) $Y_j=M(t_j;{\theta }_0)+{\mathcal{E}}_j,j=1,2,\dots ,n,$(17) correctly describes the observation process. In other words, Equation (17(17) $Y_j=M(t_j;{\theta }_0)+{\mathcal{E}}_j,j=1,2,\dots ,n,$(17) ) is the true model, and ${\theta }_0$ is the true value of the mathematical model parameter θ.

With our assumption on measurement errors, the mathematical model parameter θ can be estimated by using the OLS method; that is, the OLS estimator of θ is obtained by solving ${\theta }^n=\arg \underset{\theta \in {\mathrm{\Omega }}_{\theta }}{min}{J}^n(\theta ;\mathbf{\text{Y}}).$ Here $\mathbf{\text{Y}}=(Y_1,Y_2,\dots ,Y_n{)}^{\mathrm{T}}$ , and the cost function $J^n$ is defined as $J^n(\theta ;\mathbf{\text{Y}})=\frac{1}{n}\sum _{k=1}^n(Y_k-M(t_k;\theta ){)}^2.$ The corresponding realization ${\hat{\theta }}^n$ of ${\theta }^n$ is obtained by solving ${\hat{\theta }}^n=\arg \underset{\theta \in {\mathrm{\Omega }}_{\theta }}{min}{J}^n(\theta ;\mathbf{\text{y}}),$ where $\mathbf{\text{y}}$ is a realization of $\mathbf{\text{Y}}$ (i.e. $\mathbf{\text{y}}=(y_1,y_2,\dots ,y_n{)}^{\mathrm{T}}$).

As alluded to in the introduction, we might also consider a restricted version of the mathematical model in which the unknown true parameter is assumed to lie in a subset ${\mathrm{\Omega }}_{\theta }^H\subset {\mathrm{\Omega }}_{\theta }$ of the admissible parameter space. We assume this restriction can be written as a linear constraint, $\mathcal{H}{\theta }_0=\mathbf{\text{h}}$, where $\mathcal{H}\in {\mathbb{R}}^{{\kappa }_r\times {\kappa }_q}$ is a matrix having rank ${\kappa }_r$ (i.e. ${\kappa }_r$ is the number of constraints imposed), and $\mathbf{\text{h}}$ is a known vector. Thus the restricted parameter space is ${\mathrm{\Omega }}_{\theta }^H=\{\theta \in {\mathrm{\Omega }}_{\theta }:\mathcal{H}\theta =\mathbf{\text{h}}\}.$ Then the null and alternative hypotheses are $\begin{array}{rl}& H_0:{\theta }_0\in {\mathrm{\Omega }}_{\theta }^H,\\ & H_A:{\theta }_0\notin {\mathrm{\Omega }}_{\theta }^H.\end{array}$ We may define the restricted parameter estimator as ${\theta }^{n,H}=\arg \underset{\theta \in {\mathrm{\Omega }}_{\theta }^H}{min}{J}^n(\theta ;\mathbf{\text{Y}}),$ and the corresponding realization is denoted by ${\hat{\theta }}^{n,H}$. Since ${\mathrm{\Omega }}_{\theta }^H\subset {\mathrm{\Omega }}_{\theta }$, it is clear that $J^n({\hat{\theta }}^n;\mathbf{\text{y}})\le J^n({\hat{\theta }}^{n,H};\mathbf{\text{y}}).$ This fact forms the basis for a model selection criterion based upon the residual sum of squares.

Using the standard assumptions (given in detail in [11]), one can establish asymptotic convergence result for the test statistics (which is a function of observations and is used to determine whether or not the null hypothesis is rejected) $U^n=\frac{n\left(J^n\right({\theta }^{n,H};\mathbf{\text{Y}})-J^n({\theta }^n;\mathbf{\text{Y}}\left)\right)}{J^n({\theta }^n;\mathbf{\text{Y}})},$ where the corresponding realization ${\hat{U}}_n$ is defined as ${\hat{U}}^n=\frac{n\left(J^n\right({\hat{\theta }}^{n,H};\mathbf{\text{y}})-J^n({\hat{\theta }}^n;\mathbf{\text{y}}\left)\right)}{J^n({\hat{\theta }}^n;\mathbf{\text{y}})}.$ This asymptotic convergence result is summarized in the following theorem.

Theorem 6.1

Under assumptions detailed in [5,11] and assuming the null hypothesis $H_0$ is true, then $U^n$ converges in distribution $($as $n\to \mathrm{\infty })$ to a random variable U having a chi-square distribution with ${\kappa }_r$ degrees of freedom.

The above theorem suggests that if the sample size n is sufficiently large, then $U^n$ is approximately chi-square distributed with ${\kappa }_r$ degrees of freedom. We use this fact to determine whether or not the null hypothesis $H_0$ is rejected. To do that, we choose a significance level α (usually chosen to be 0.05) and use ${\chi }^2$ tables to obtain the corresponding threshold value τ so that $\mathrm{Prob}(U>\tau )=\alpha$. We next compute ${\hat{U}}^n$ and compare it to τ. If ${\hat{U}}^n>\tau$, then we reject the null hypothesis $H_0$ with confidence level $(1-\alpha )100\mathrm{\%}$; otherwise, we do not reject. We emphasize that care should be taken in stating conclusions: we either reject or do not reject $H_0$ at the specified level of confidence. The following table illustrates the threshold values for ${\chi }^2\left(1\right)$ with the given significance level.

Table

α	τ	Confidence level
.25	1.32	75%
.1	2.71	90%
.05	3.84	95%
.01	6.63	99%
.001	10.83	99.9%

Similar tables can be found in any elementary statistics text or online or calculated by some software package such as Matlab, and is given here for illustrative purposes and also for use in the examples demonstrated in the following.

6.2. Generalized least squares

The model comparison results outlined can be extended to deal with GLS problems in which measurement errors are independent with $\mathbb{E}\left({\mathcal{E}}_k\right)=0$ and $\mathrm{Var}\left({\mathcal{E}}_k\right)={\sigma }_0^2{w}^2(t_k,\hat{\theta })$, $k=1,2,\dots ,n$, where w is some known real-valued function with $w(t,\hat{\theta })\ne 0$ for any t. This is achieved through rescaling the observations in accordance with their variance (as discussed in [11]) so that the resulting (transformed) observations are identically distributed as well as independent.

6.3. Results for PolyQ aggregation models

We then carried out a series of model comparison tests (we again used DS4) for nested models to determine if an added parameter yields a statistically significantly improved model fit. Our null hypothesis in each case was: $H_0$: The restricted model is adequate (i.e. the fit to data is not significantly improved with the model containing the additional parameter as a parameter to be estimated). We obtained the following results.

1. The model with estimation of $\{k_I^{+},k_I^{-}\}$ vs. the model with estimation of $\{k_I^{+},k_I^{-},k_{\mathrm{off}}^N\}$: we find with $n=699$, $J_n({\hat{\theta }}_H^n;\mathbf{\text{Y}})=0.0044192109$, $J_n({\hat{\theta }}^n;\mathbf{\text{Y}})=0.0043709501,$ and ${\hat{U}}_n=7.7178$. Thus, we reject $H_0$ at a 99% confidence level.

2. The model with estimation of $\{k_I^{+},k_I^{-}\}$ vs. the model with estimation of $\{k_I^{+},k_I^{-},k_{\mathrm{on}}^N\}$: we find $J_n({\hat{\theta }}^n;\mathbf{\text{Y}})=0.0044192108$ with ${\hat{U}}_n=7.49\times {10}^{-06}.$ Thus, we do not reject $H_0$ at a 99% confidence level.

3. The model with estimation of $\{k_I^{+},k_I^{-},k_{\mathrm{off}}^N\}$ vs. the model with estimation of $\{k_I^{+},k_I^{-},k_{\mathrm{off}}^N,k_{\mathrm{on}}^N\}$: to the order of computation we find no difference in the cost functions in this case and therefore we do not reject $H_0$ at a confidence level of 99%.

4. The model with estimation of $\{k_I^{+},k_I^{-},k_{\mathrm{on}}^N\}$ vs. the model with estimation of $\{k_I^{+},k_I^{-},k_{\mathrm{on}}^N,k_{\mathrm{off}}^N\}$: we find $J_n({\hat{\theta }}^n;\mathbf{\text{Y}})=0.0043709780$ with ${\hat{U}}_n=7.7133$ and hence we reject $H_0$ with a confidence level of 99%.

From these and the preceding results we conclude that the information content of the typical data set for the dynamics considered here will support at most three parameters estimated with reasonable confidence levels and these are the parameters $\{k_I^{+},k_I^{-},k_{\mathrm{off}}^N\}$.

7. Conclusions and suggested further efforts

For the efforts reported on above we make several conclusions.

For the majority of data sets, the GLS residual plots with $\gamma =0.6$ are random when fitted for data points $M\left(t_k\right)\ge \mathrm{0.12.}$ As conjectured earlier, this may be because the early formation of aggregates is somewhat stochastic in nature which is not well described by either the mathematical and/or statistical models. It appears that one needs special consideration of smaller polymer sizes. Indeed we suspect from additional discussions with our colleagues that perhaps the nucleation step might be dominated by a stochastic rather than deterministic process in the early stages (i.e. for small polymer sizes). This is a possible direction of further investigation.

Based on several different mathematical/statistical methodologies (sensitivities, asymptotic analysis, bootstrapping, and model comparison tests), the data sets we considered do not contain sufficient information for the reliable estimation of all nine parameters of interest. Indeed our findings suggest that at most three parameters can be reliably estimated with the data sets typical of those presented here, and that these parameters are $\{k_I^{+},k_I^{-},k_{\mathrm{off}}^N\}$. Recently related efforts [2] suggest that perhaps there are experimental design questions that could be addressed to collect data that might support the more sophisticated models derived in [25], especially in order to investigate information coming from different initial concentrations. Indeed, we have considered here data sets related to experiments carried out with the same initial concentration. Adapting the previously used techniques to simultaneously or successively use all the information content in data sets carried out for different initial concentration is a challenging problem (see [24] for a discussion of the effect of initial concentration on nucleated polymerization).

Here we concluded that at most three parameters $\{k_I^{+},k_I^{-},k_{\mathrm{off}}^N\}$ can be reliably estimated with the data sets investigated. The two first parameters determine the balance between the normal and abnormal protein concentrations and the third represents the stability of the nucleus against the degradation into monomeric entities. These three parameters are related to the early steps of the aggregation process, and thus we conclude that the model applied to these data sets does not provide any insight into the polymerization of larger polymers. Since this is the case, there is little motivation to modify the polymerization function $k_{\mathrm{on}}$ depicted in Figure 2 until further data collection procedures are pursued. It is difficult to say whether the initiation parameters are biologically more or less interesting than other parameters in the later aggregation process. (Everything interests the biological community since very little is known with certainty.) What we can more prudently conclude is rather a negative conclusion: this type of experiments seems to be more sensitive to the ignition of the reaction than to secondary pathways, which reveals a limit, and argues for the measurement not only of the total polymerized mass but also of size distributions of fibrils if one wants to really select the proper model for self-acceleration, as for instance in [26,30,31].

The methodology we developed here may be carried out for other types of proteins as well as related experiments, and this is of peculiar interest since the ThT measurements as studied here are the most standard method for protein aggregation. For instance, in the seminal paper [29], models are proposed and fit to data; our type of analysis would complement their findings and allow one to assess the quality of these fits and how much confidence the authors can have in their conclusions. These results in turn could be combined with optimal design techniques [2,4,11] to design more informative experiments.

Acknowledgments

The authors are most grateful to a referee whose careful reading of an earlier version of this manuscript led to substantial improvements in the final version.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research was supported in part (MD, CK) by the ERC Starting Grant SKIPPERAD, in part (HTB) by the National Institute of Allergy and Infectious Diseases [grant number NIAID R01AI071915-10], and in part (HTB) by the Air Force Office of Scientific Research [grant number AFOSR FA9550-12-1-0188].