Genreal Remarks

The following content provides an excerpt from a larger scientific work in the field of molecular structural biochemistry focusing on post-transcriptional transfer RNA modifications. It highlights the statistical modeling of the interaction between the 4-thiouridine synthetase ThiI and one of its target tRNA substrates in order to unravel the underlying stoichiometry of the complex that is formed. For this reason a strongly reduced introductory section is presented not going into detail with respect to the structural organization of the macromolecules involved, the photochemistry of 4-thiouridine and its biological role, as well as the underlying persulfide chemistry. Instead, the work focuses on model fitting by non-linear least squares regression, model validation based on maximum likelihood profiling, model evaluation procedures such as Kullback-Leibler divergence analysis and cross validation, as well as generic and robust ANOVA approaches.

Statistical computations were accomplished by using R. Code chunks accompanying the respective sections can be switched on if desired. These provide detailed code lines necessary in order to reproduce the analysis steps.

Some methods are explained in more detail to promote understanding especially when it comes to bridge the gap between classical uderstanding of ANOVA and the general linear regression model in the context of zero hypothesis testing. Thus, the strict framework of purely scientific writing and presentation is somewhat relaxed.

1 Introduction

Transfer RNAs are macromolecular adaptor molecules and central key players of the protein biosynthesis machinery. Serving as a physikal link during translation on the ribosome, they accomplish the important task of specifying the primary amino acid sequence of a given protein as determined by the mRNA transcript of its open reading frame. In order to properly fulfil this task with high fidelity the canonical ribonucleoside building blocks of tRNAs are to a large extent post-transriptionally modified. This one of the processing steps that convert premature tRNA transcripts into the mature and biologically active form. Over 80 distinct chemical modifications have already been reported for tRNAs which is unprecedented among ribonucleic acids. This extraordinary broad and diverse set of modifications is introduced by dedicated tRNA modifying enzymes and provides tRNAs with the chemical upgrade necessary so that they can properly perform at the various stages of the translational process (El Yacoubi, Bailly, and Crécy-Lagard 2012). Thus, this finetuning is mandatory inter alia for the structural integrity of tRNAs, codon-anticodon specificity, reading-frame maintenance, and for the tRNA binding to the ribosoe (Shigi 2014). A special subset of modified ribonucleosides which makes important contributions to all these processes comprise the so-called thionucleosides into which a sulfur atom is incorporated. Among the 16 thionucleosides that are currently known (Crain and McCloskey 1996) 4-thiouridine (s4U) is characterized by unique spectroscopic and photochemical properties that affect the structural integrity of tRNAs. s4U is invariantly found at position 8 of cytosolic tRNAs from Archaea and Bacteria (Griffey et al. 1986) and its biological role is to act as a physiological photosensor for near-UV radiation exposure, thereby providing protection against its damaging effects (Traxler et al. 2008; Thiam and Favre 1984). The biosynthesis of s4U8 necessitates reactive persulfide chemistry as well as the accurate recognition of the target base U8 that is burried in a highly structured tRNA molecule. This biochemically demanding task is accomplished by the cooperative activity of the cysteine desulfurase IscS and the 4-thiouridine synthetase ThiI. Recently the crystal structure of T. maritima ThiI in complex with TPHE39A, a truncated version of E. coli tRNAPhe, has been solved indicating that dimeric ThiITm binds two TPHE39A substrates (Neumann et al. 2014). This underlying 1:2 stoichiometry is re-investigated in this work using ThiITm and the full-length tRNAPhe substrate.

1.1 4-thiouridine synthetases

Bacterial 4-thiouridine synthetases are both persulfide acceptors as well as tRNA targeting enzymes which carry out the formation of s4U8 in cytosolic tRNAs. To accomplish this task they catalyze the site-specific sulfuration at position four of U8 as well as its pre-activation via adenylation. ThiI proteins are characterized by a functional modularization according to which a conserved catalytic domain is combined with additional domains that mediate substrate recognition. The majority of eubacterial ThiI proteins contain three functional domains. To this belongs the ThiI protein of Thermotoga maritima (ThiITm,PDB 4KR7) whose structure has recently been solved by X-ray crystallography (Neumann et al. 2014). In this structure ThiITm is present as a homodimer. Each of the monomers has an elongated, tripartite organization composed of an N-terminal ferredoxin-like domain (NFLD), a THUMP domain, and a C-terminal PP-loop ATP pyrophosphatase (PPase) domain. Both the NFLD and THUMP domains together accomplish the recognition of the respective tRNA substrate whereas the catalytic function resides in the C-terminal PPase domain which also binds the co-substrate ATP (shown as balls and sticks in figure 1.1B) at the active site and which includes the flexible loop (yellow) harboring the catalytically essential Cys344 (orange). ThiITm has additionally been successfully crystallized in complex with a truncated version of tRNAPhe. These complexes give a first insight into the mechanisms of specific target base recognition and inter-domain communication between the substrate binding NFLD/THUMP domains and the catalytic PPase domain.

1.2 ThiITm in complex with truncated tRNAPhe (TPHE39A)

ThiI catalyzes the s4U8 modification of all cytosolic tRNAs. For this reason it has to employ a universal mechanism that allows specific target base recognition that is independent from specific RNA sequence motifs. Crystal structures obtained from ThiITm in complex with truncated tRNAPhe (TPHE39A) reveal that ThiITm acts as a molecular ruler which uses the 3’-ACCA end common to all target tRNAs as a reference point in order to properly place the target base in the active site pocket (Neumann et al. 2014).

TPHE39A is a 39 nucleotide tRNA that represents the minimal efficient substrate for ThiI and was determined based on deletion analysis of full-length tRNAPhe (Lauhon, Erwin, and Ton 2004). The most important binding determinant is the 3’-ACCA end as an invariant primary structural element (Neumann et al. 2014). The 3’-ACCA end is bound in a groove formed by distinct secondary structural elements of the THUMP domain which establish an interaction network including hydrogen bonds as well as polar-, and van-der-Waals contacts, respectively. This highly specific, THUMP-mediated binding of the 3’-ACCA overhang which is common to all tRNA substrates is the basis of the structural motif- and sequence independent recognition of the target base U8. That is, the 3’-ACCA end serves as one reference point of a molecular ruler that allows to define and measure the distance to the site of modification (1.1A). In this way, one RNA molecule bound by the THUMP domain of one monomer places the U8-containing bulge region of TPHE39A in proximity to the catalytic center of the PPase domain of the other monomer. Thus, only dimeric ThiITm is capable of binding and modify TPHE39A. Remarkably, alternative active site conformations in the monomers of a given ThiITm•TPHE39A complex are observed.

These structurally asymmetric active sites can be traced back to alternative conformations of four active site loops (1.1B). The most striking changes are observed in the so-called bottom loop (blue) belonging to the NFLD of one monomer and the Cys344-containing upper loop (yellow) of the PPase domain of the opposing monomer. A change of the loop conformations seems to be induced upon proper binding of the RNA molecule. This coincides with a switch of the active site pocket from an open to a closed state which might reflect the transition from a substrate binding state to a base flipping competent state. The underlying inter-domain communication via flexible loops which signal a properly placed RNA molecule is the basis of the functional modularity of ThiITm. In the closed state the target base U8 (wite balls and sticks) and Cys344 are directed towards the active site and placed in the vicinity of ATP which is kept in a reactive conformation. The fact that ThiITm•TPHE39A complexes are functional classifies ThiITm as a group I tRNA modification enzyme according to Lauhon, Erwin, and Ton (2004). Members of this group do not require structurally intact tRNA substrates and upon binding these substrates are significantly remodeled. For a successful sulfuration and activation of U8 such a remodeling event might be necessary in order for ThiITm to gain access to the target base which is buried in a highly structured tRNA when full-length substrates are considered.

The ThiITm•TPHE39A complex is currently the only structure available concerning the s4U8 biosynthesis in T. maritima and in general. That is why additional crystal structures of functional binary and ternary complexes are needed in order to fully elucidate the s4U8 formation from a structural perspective. One critical aspect that would be solved by crystal structures but which can also be accessed differently is that of the stoichiometry of complexes formed between ThiITm and full-length tRNAs. This information can also well be extracted from titration experients as will be shown in this work (2).

**ThiI<sub>Tm</sub> in complex with TPHE39A.** Shown are a schematic represantion of monomeric ThiI<sub>Tm</sub> (A) and a surface representation (B) of dimeric ThiI<sub>Tm</sub> in complex with truncated TPHE39A (PDB-ID 4KR7). **(A):** Principles of U8 recognition. ThiI<sub>Tm</sub> acts as a molecular ruler in order to recognize the target base U8 independent of sequence- and structural motifs of the substrate tRNA. Therefore it uses the 3'-ACCA overhang which is common to all tRNAs as a reference point. **(B):** Domain communication in ThiI<sub>Tm</sub>. Monomers ([A] and [B]) are colored in scales of red and grey. From the point of view of monomer [A], the 3'ACCA-overhang of TPHE39A is recognized by the, respectively, NFLD[A] and THUMP[A] domains. In this way the target base U8 (white balls and sticks) is placed in the vicinity of the active site of the PPase[B] domain. Its presence is directly communicated via the bootom loop (blue) of NFLD[A] to the catalytically active Cys344 (orange) bearing upper loop (yellow) of the PPase[B] domain, thereby by-passing the PPase[A] domain. This coincides with a switch of the active site pocket from an open to a closed conformation which is supported by adjoining flexible loop regions (green and cyan).

Figure 1.1: ThiITm in complex with TPHE39A. Shown are a schematic represantion of monomeric ThiITm (A) and a surface representation (B) of dimeric ThiITm in complex with truncated TPHE39A (PDB-ID 4KR7). (A): Principles of U8 recognition. ThiITm acts as a molecular ruler in order to recognize the target base U8 independent of sequence- and structural motifs of the substrate tRNA. Therefore it uses the 3’-ACCA overhang which is common to all tRNAs as a reference point. (B): Domain communication in ThiITm. Monomers ([A] and [B]) are colored in scales of red and grey. From the point of view of monomer [A], the 3’ACCA-overhang of TPHE39A is recognized by the, respectively, NFLD[A] and THUMP[A] domains. In this way the target base U8 (white balls and sticks) is placed in the vicinity of the active site of the PPase[B] domain. Its presence is directly communicated via the bootom loop (blue) of NFLD[A] to the catalytically active Cys344 (orange) bearing upper loop (yellow) of the PPase[B] domain, thereby by-passing the PPase[A] domain. This coincides with a switch of the active site pocket from an open to a closed conformation which is supported by adjoining flexible loop regions (green and cyan).

2 Goals

In this work the biochemichal characterization of the s4U8 synthesizing enzyme system covers interaction studies concerning the chimeric ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) sub-complex and extends structural investigations on the ThiITm•TPHE39A complex done by Neumann et al. (2014). This crystal structure revealed that each monomer of dimeric ThiITm binds one truncated TPHE39A substrate. However, alternative active site conformations in the monomers could be observed and it remains elusive whether this asymmetry is part of the catalytic mechanism that would also apply for full-length tRNA substrates or whether this might lead to a half-of-the-sites reaction mechanism which has recently been described for the methyltransferase TrmD (Christian et al. 2010), another tRNA modifying enzyme. In addition, a structural model obtained from the superposition of full-length tRNA\(\mathrm{^{Phe}_{Ec}}\) onto the ThiITm•TPHE39A structure (Neumann et al. 2014) revealed that substantial conformational changes of a full-length tRNA substrate can be expected upon binding which might impair simultaneous binding of two substrates due to potential steric clashes. Against this backdrop the leading question of this work is that of the true underlying stoichiometry of complexes between ThiITm and full-length tRNA substrates and is going to be practically addressed in terms of fluorescence spectroscopic binding studies using the interaction between ThiITm and B. subtilis tRNA\(\mathrm{^{Phe}_{Bs}}\) as an example. A comparative multi-model approach is going to be pursued aimed at finding the statistically best approximating binding model of this interaction in a set of alternative binding models all related to underlying 1:1 (ThiITm : tRNA\(\mathrm{^{Phe}_{Bs}}\)) or 1:2 binding equilibria.

3 Methods

Notice that methods describing the expression and purification protocols of the single tRNA and protein components as well as methods designed to validate stable complex formation were omitted. Only methods involving the key experiment and the following statistical analysis workflow are described here.

3.1 Bioanalytical Methods

3.1.1 Tryptophan Fluorescence Spectroscopy

Complex formation of ThiITm and tRNA\(\mathrm{^{Phe}_{Bs}}\) was analyzed based on quenching of the intrinsic tryptophan fluorescence of ThiITm. Five μM of ThiITm was titrated with incrementally increasing amounts of tRNA\(\mathrm{^{Phe}_{Bs}}\) ranging from 0.03 μM to 60 μM in a buffer containing 50 mM Tris/HCl pH 7.5, 300 mM NaCl, 2 mM DTT, and 5% glycerol so that final volumes of 40 μL were obtained. Reaction set ups were incubated for 15 min at 600 rpm and measured in triplicates. Fluorescence spectra were recorded at room temperature on a Victor Nivo Multimode Plate Reader (PerkinElmer) in 96 well round bottom plates using an excitation filter of 280/10 nm and an emission filter of 355/40 nm. Obtained binding isotherms were fitted by non-linear least squares regression as described in 3.2.1. The change in the observed fluorescence intensity \(\Delta F_{obs}\) was transformed to reflect the proportional signal quenching:

\[\begin{equation} \Delta F_{obs,i} = \frac{\Delta F_{obs,max} - \Delta F_{obs,i}}{\Delta F_{obs,max}} \tag{3.1} \end{equation}\]

where \(\Delta F_{obs,max}\) is the maximum achievable signal provided by the negative control when no tRNA induced signal quenching occurs.

3.2 Mathematical and Computational Tools

3.2.1 Non-Linear Least Squares Analysis of Binding Isotherms

Binding isotherms obtained from tryptophan fluorescence measurements (3.1.1) were analyzed via ordinary nonlinear least squares (nls) regression using the Levenberg-Marquardt algorithm implemented in R (R Core Team 2018). Data were comparatively fitted to different binding equilibria describing ThiI:tRNA binding stoichiometries of 1:1 or 1:2, respectively. All following 1:2 binding models ((3.2), (3.4), (3.5), and (3.6)) are based on the Adair-Klotz formulism (Klotz 1946) describing sequential binding to independent binding sites. Among these the global binding model (least degrees of freedom) which is referred to as the bm1to2.coop model throughout this work is described by equation (3.2):

\[\begin{equation} \Delta F_{obs}=\frac{[H]_0\left(F_{\Delta HG}\,K_1[G]+F_{\Delta HG2}\,K_1K_2[G]^2\right)}{1+K_1[G]+K_1K_2[G]^2} \tag{3.2} \end{equation}\]

Where \([H]_0\) = total concentration of Host (ThiI), \([G]\) = concentration of free Guest (tRNA) at equilibrium, \(F_{\Delta HG}\)/\(F_{\Delta HG2}\) = proportionality factors relating the change in fluorescence signal (\(\Delta F_{obs}\)) to the concentration of the complex having one or two guests bound, respectively, and \(K_1\)/\(K_2\) = macroscopic binding constants of the first and second binding step, respectively. This model accounts for cooperativity which can explicitly be described by the interaction parameter \(\alpha\) according to equation (3.3):

\[\begin{equation} \mathrm{\alpha=\frac{K_1}{4K_2}=} \begin{cases} >1 & \text{positive coop.}\\ 1 & \text{no coop.}\\ <1 & \text{negative coop.} \end{cases} \tag{3.3} \end{equation}\]

Depending on the relationship between the binding constants and the proportionality factors equation (3.2) can be re-parameterized in different ways yielding nested 1:2 binding models having more degrees of freedom. If the proportionality factors are assumed to be additive (\(F_{\Delta HG2}\) = 2\(F_{\Delta HG}\)) equation (3.2) reduces to the bm1to2.coop.add (“cooperative/additive model”) described by equation (3.4):

\[\begin{equation} \Delta F_{obs}=\frac{[H]_0\,F_{\Delta HG}\left(K_1[G]+2K_1K_2[G]^2\right)}{1+K_1[G]+K_1K_2[G]^2} \tag{3.4} \end{equation}\]

Assuming non-cooperativity between the binding sites which can also be referred to as degenerate binding sites, then \(K_1\) = 4\(K_2\) according to equation (3.3). Both binding constants can then be combined to one overall stability constant \(\beta_{12}\) = \(K_1K_2\) so that equation (3.2) simplifies to the “degenerate/non-additive model” bm1to2.deg given by equation (3.5):

\[\begin{equation} \Delta F_{obs}=\frac{[H]_0 \left(F_{\Delta HG}\,2\sqrt{\beta_{12}}\,[G]+F_{\Delta HG2}\,\beta_{12}[G]^2\right)}{1+2\sqrt{\beta_{12}}\,[G] +\beta_{12}[G]^2} \tag{3.5} \end{equation}\]

When non-cooperativity between the binding sites and additivity between the proportionality factors are assumed equation (3.2) reduces to the most simple “degenerate/additive model” bm1to2.deg.add described by equation (3.6):

\[\begin{equation} \Delta F_{obs}=\frac{[H]_0\, F_{\Delta HG}\left(2\sqrt{\beta_{12}}\,[G]+2\beta_{12}[G]^2\right)}{1+2\sqrt{\beta_{12}}\,[G] +\beta_{12}[G]^2} \tag{3.6} \end{equation}\]

Every 1:2 binding model depends on free guest concentration \([G]\). This quantity needed to be derived from the corresponding total concentration \([G]_0\) used in each titration step of the binding experiment. That was accomplished by solving the cubic polynomial in equation (3.7) having the coefficients \(A\), \(B\), \(C\), and \(D\) given below:

\[\begin{gather} \mathrm{(A)[G]^3+(B)[G]^2+(C)[G]-D=0} \quad \text{,with}\\ \mathrm{A=K_1K_2} \notag\\ \mathrm{B=K_1(2K_2[H]_0-K_2[G]_0+1)} \notag\\ \mathrm{C=K_1([H]_0-[G]_0)+1} \notag\\ \mathrm{D=[G]_0} \notag \tag{3.7} \end{gather}\]

The cubic was solved by the root.finder function (see code chunk below). Generally, \([H]_0\) was treated as a constant and set to 5 \(\mu\)m for all fitting routines. For a given vector of 12 starting concentrations \(\overrightarrow{G}_0\) and a set of parameter values (\(K_1\), \(K_2\)) this function yields at most three solutions for every entry of the vector \(\overrightarrow{G}_0\) which might be positive, negative, or complex. The trna_eq function was used to extract the smallest positive real roots of the polynomial so that a vector of 12 equilibrium concentrations \(\vec{G}\) was obtained.

The set of model functions was completed by the 1:1 binding model bm1to1 described by equation (3.8): \[\begin{equation} \mathrm{ \Delta F_{obs}=F_{\Delta HG} \left(\frac{1}{2}\left\{ \left([H]_0+[G]_0+\frac{1}{K_1}\right)-\sqrt{\left([H]_0+[G]_0+\frac{1}{K_1}\right)^2-4[H]_0[G]_0} \right\} \right)} \tag{3.8} \end{equation}\]

3.2.2 Implementation of the nls fitting routine

In case of the two-site binding models, both, the predicted values \(\Delta F_{obs}\) as well as the explanatory variable \([G]\) are dependent on the parameter values of \(K_1\) and \(K_2\) which are supposed to be optimized. Theoretically, each iteration of the Levenberg-Marquardt algorithm (the nlsLM() function) would yield a different vector of the equilibrium concentration of the tRNA (\(\vec{G}\)). This parameter dependent change of the explanatory varibale is not accounted for by the generic nlsLM() function. Due to this, two wrapper functions were written, one for the two degenerative binding models (var_nls_deg) and one for the two cooperative binding models (var_nls_coop) whose core structure is schematically depicted in figure 3.1. Briefly, the var_nls_*() functions take the same arguments as the nlsLM() function. During the initialization phase the point estimates of \(K_1\) and \(K_2\) are firstly found based on the vector of starting concentrations \(\overrightarrow{G}_0\) (start.conc). During the while loop structure the trna_eq function is applied to find the corresponding vector of equilibrium tRNA concentrations \(\overrightarrow{G}\) given the initially found point estimates of \(K_1\) and \(K_2\). This vector is then used as the new explanatory variable and the binding models are re-optimized. The obtained \(RMSE\) (\(RMSE_{new}\)) is compared to the \(RMSE\) of the previous iteration (\(RMSE_{curr}\)) until no further improvement is made.

var_nls_deg <- function(..., start.conc){
  nls.model <- nlsLM(...)
  
  g.current <- start.conc
  k1.current <- as.numeric(coef(nls.model)["k1"])
  # max.iter <- 0
  sigma.new <- 2
  sigma.curr <- 1
  sigma.it <- c(sigma(nls.model) + 1, sigma(nls.model))
  results <- list()
  results$sigma.model[[1]] <- sigma(nls.model)   
  results$k1[[1]] <- k1.current
  results$trna.conc[[1]] <- g.current
  
  while(sigma.it[sigma.curr] >= sigma.it[sigma.new]){
    if(length(coef(nls.model)) == 2){ 
      model.loop <- nlsLM(que_fi ~ bm1to2.deg.add(fhg, k1, g.current),
                          data = fi,
                          start = list(fhg = as.numeric(coef(nls.model)["fhg"]),
                                       k1 = k1.current))
    } else {
      model.loop <- nlsLM(que_fi ~ bm1to2.deg(fhg, fhg2, k1, g.current),
                          data = fi,
                          start = list(fhg = as.numeric(coef(nls.model)["fhg"]),
                                       fhg2 = as.numeric(coef(nls.model)["fhg2"]),
                                       k1 = k1.current))
    }
    
    sigma.curr <- sigma.curr +1
    sigma.new <- sigma.new +1
    sigma.it <- c(sigma.it, sigma(model.loop))
    g.current <- trna_eq(constants = c(coef(model.loop)["k1"], 
                                       coef(model.loop)["k1"]/4, 
                                       5),
                         start.conc)
    k1.current <- as.numeric(coef(model.loop)["k1"])
    print(sigma(model.loop))
    
    results$sigma.model[[sigma.curr]] <- sigma(model.loop)
    results$k1[[sigma.curr]] <- k1.current
    results$trna.conc[[sigma.curr]] <- g.current
    
  }
  
  selector <- which(results$sigma.model == min(results$sigma.model)) - 1 
  g.best <- unlist(results$trna.conc[selector])
  k1.best <- unlist(results$k1[selector])
  ###
  results.list <- list()
  if(length(coef(nls.model)) == 2){
    nls.model.best <- nlsLM(
      que_fi ~ bm1to2.deg.add(fhg, k1, g0 = g.best),
      data = fi,
      start = list(fhg = as.numeric(coef(nls.model)["fhg"]),
                   k1 = k1.best))
    
    results.list <- list(bm1to2.deg.add = nls.model.best,
                         g = g.best) 
    
  } else {
    nls.model.best <- nlsLM(que_fi ~ bm1to2.deg(fhg, fhg2, k1, g.best),
                            data = fi,
                            start = list(fhg = as.numeric(coef(nls.model)["fhg"]),
                                         fhg2 = as.numeric(coef(nls.model)["fhg2"]),
                                         k1 = k1.best))
    results.list <- list(bm1to2.deg = nls.model.best,
                         g = g.best) 
    
  }
  
}
var_nls_coop <- function(..., start.conc){
  nls.model <- nlsLM(...)
  
  g.current <- start.conc
  k1.current <- as.numeric(coef(nls.model)["k1"])
  k2.current <- as.numeric(coef(nls.model)["k2"])
  # max.iter <- 0
  sigma.new <- 2
  sigma.curr <- 1
  sigma.it <- c(sigma(nls.model) + 1, sigma(nls.model))
  results <- list()
  results$sigma.model[[1]] <- sigma(nls.model)   
  results$k1[[1]] <- k1.current
  results$k2[[1]] <- k2.current
  results$trna.conc[[1]] <- g.current
  
  while(sigma.it[sigma.curr] >= sigma.it[sigma.new]){
    if(length(coef(nls.model)) == 3){ 
      model.loop <- nlsLM(que_fi ~ bm1to2.coop.add(fhg, k1, k2, g.current),
                          data = fi,
                          start = list(fhg = as.numeric(coef(nls.model)["fhg"]),
                                       k1 = k1.current,
                                       k2 = k2.current))
    } else {
      model.loop <- nlsLM(que_fi ~ bm1to2.coop(fhg, fhg2, k1, k2, g.current),
                          data = fi,
                          start = list(fhg = as.numeric(coef(nls.model)["fhg"]),
                                       fhg2 = as.numeric(coef(nls.model)["fhg2"]),
                                       k1 = k1.current,
                                       k2 = k2.current))
    }
    sigma.curr <- sigma.curr +1
    sigma.new <- sigma.new +1
    sigma.it <- c(sigma.it, sigma(model.loop))
    g.current <- trna_eq(constants = c(coef(model.loop)["k1"], 
                                       coef(model.loop)["k2"], 
                                       5),
                         start.conc)
    k1.current <- as.numeric(coef(model.loop)["k1"])
    k2.current <- as.numeric(coef(model.loop)["k2"])
    
    print(sigma(model.loop))
    results$sigma.model[[sigma.curr]] <- sigma(model.loop)
    results$k1[[sigma.curr]] <- k1.current
    results$k2[[sigma.curr]] <- k2.current
    results$trna.conc[[sigma.curr]] <- g.current
    
  }
  
  selector <- which(results$sigma.model == min(results$sigma.model)) - 1 
  g.best <- unlist(results$trna.conc[selector])
  k1.best <- unlist(results$k1[selector])
  k2.best <- unlist(results$k2[selector])
  results.list <- list()
  
  if(length(coef(nls.model)) == 3){
    nls.model.best <- nlsLM(que_fi ~ bm1to2.coop.add(fhg, k1, k2, g.best),
                            data = fi,
                            start = list(fhg = as.numeric(coef(nls.model)["fhg"]),
                                         k1 = k1.best,
                                         k2 = k2.best))
    results.list <- list(bm1to2.coop.add = nls.model.best,
                         g = g.best) 
  } else {
    nls.model.best <- nlsLM(que_fi ~ bm1to2.coop(fhg, fhg2, k1, k2, g.best),
                            data = fi,
                            start = list(fhg = as.numeric(coef(nls.model)["fhg"]),
                                         fhg2 = as.numeric(coef(nls.model)["fhg2"]),
                                         k1 = k1.best,
                                         k2 = k2.best))
    results.list <- list(bm1to2.coop = nls.model.best,
                         g = g.best) 
  }
}
**Extension of the basic nlsLM() function in order to fit two-site binding models.** Shown is a schematic representation of the core structure of the var_nls_*() functions designed to handle optimizations of models where the explanatory variable itself is changing with the parameter to be optimized. Descriptions are given in the text.

Figure 3.1: Extension of the basic nlsLM() function in order to fit two-site binding models. Shown is a schematic representation of the core structure of the var_nls_*() functions designed to handle optimizations of models where the explanatory variable itself is changing with the parameter to be optimized. Descriptions are given in the text.

3.2.3 Maximum Likelihood Profiling

Point estimates of model-specific parameter obtained by non linear least squares (nls) analysis (3.2.1) were tested for parameter identifiability based on maximum likelihood profiling as proposed by Raue et al. (2010). In case of statistically independent and additive Gaussian noise as well as constant variance (\(\epsilon\) ~ \(N(0,\sigma^2)\)), the logarithmic form of the likelihood function \(\mathcal{L}(*)\) for a given binding model \(g(x_i, \theta)\) with parameter \(\theta\) is given by equaion (3.9):

\[\begin{equation} log(\mathcal{L}(\theta| x,g)) = -\frac{n}{2}log(2\pi\sigma^2)-\frac{1}{2\sigma^2}\sum\limits_{i=1}^n (y_i-g(x_i,\theta))^2 \tag{3.9} \end{equation}\]

When equation (3.9) is evaluated for a given binding model when its parameters are set to the best fit point estimates determined by the nls regression analysis, \(\hat{\theta}_{nls}\), then the maximum of the log-likelihood function (\(log(\mathcal{L}(\hat{\theta}_{nls}))\)) is obtained since maximizing the likelihood function is equivalent to minimizing the residual sum of squares in a nls sense. In classical hypothesis testing (\(log(\mathcal{L}(\hat{\theta}_{nls}))\)) represents the augmented or full model. When a parameter \(\theta_i\) is iteratively changed around its best fit nls point estimate (\(\hat{\theta}_{i,nls}\)) according to values specified by a vector \(\overrightarrow{p}^t = (\theta_{i,1}, ... , \hat{\theta}_{i,nls}, ..., \theta_{i,n})\), and all other parameters \(\theta_{i\neq j}\) are re-optimized based on equation (3.9), then the profile log-likelihood function \(\mathcal{P}\mathcal{L}_i(\overrightarrow{p})\) for parameter \(\theta_i\) can be obtained according to equation (3.10):

\[\begin{equation} \mathcal{P}\mathcal{L}_i(\overrightarrow{p})=\operatorname*{arg\,max}_{\left\{\theta| \theta_i=p_i\right\}} \,log(\mathcal{L}(\theta| x,g)) \tag{3.10} \end{equation}\]

\(\mathcal{P}\mathcal{L}_i(\overrightarrow{p})\) represents a nested version of the augmented model where the parameter which is profiled out is treated as a constant. Hence, the deviances between the augmented model and each of the nested models, \(2log(\mathcal{L}(\hat{\theta}_{nls}))\) - \(2\mathcal{P}\mathcal{L}_i(\overrightarrow{p})\), are focused 1 degrees of freedom (\(df\)) comparisons that thus all belong to the same \(\chi^2\)-distribution with 1 \(df\) (null distribution). Based on this distribution the corresponding 95% profile log-likelihood confidence interval (\(CI\)) for parameter \(\theta_i\) is given by equation (3.11):

\[\begin{equation} CI_{\mathcal{P}\mathcal{L}(\theta_i)=}\left\{ p_i| \mathcal{P}\mathcal{L}_i(p_i) \geq log(\mathcal{L}(\hat{\theta}_{nls}))-\frac{\chi^2_{0.95,1}}{2} \right\} \tag{3.11} \end{equation}\]

The \(CI\) contains all parameter vlaues of \(\theta_i\) for which the corresponding profile log-likelihood function exceeds a threshold (\(\Delta_\alpha\) in figure 3.2) determined by the critical value according to a significance level of .05 (\(\chi^2_{.95,1}\)). In this way the construction of the \(CI\) can be thought of as a reverse likelihood ratio test. If a parameter is identifiable, the corresponding log-likelihood profile (figure 3.2A) exhibits a unique maximum and a finite \(CI\) (red line). On the other hand, a parameter can be unidentifiable in one of two ways: A structurally unidentifiable parameter (B) has a flat log-likelihood profile with an infinite \(CI\), while a practically unidentifiable parameter shows a unique maximum of the log-likelihood profile but a \(CI\) that is only bound on one site (C). This concept can be extended for two parameters so that joint profile log-likelihood contours and confidence regions can be obtained as described by equations (3.12) and (3.13):

\[\begin{equation} \mathcal{P}\mathcal{L}_{ij}(\overrightarrow{p},\overrightarrow{q})=\operatorname*{arg\,max}_{\left\{\theta| \theta_i=p_i,\theta_j=q_i \right\}} log(\mathcal{L}(\theta| x,g) \tag{3.12} \end{equation}\]

\[\begin{equation} CI_{\mathcal{P}\mathcal{L}(\theta_i,\theta_j)=}\left\{ p_i,q_i| \mathcal{P}\mathcal{L}_{ij}(p_i,q_i) \geq log(\mathcal{L}(\hat{\theta}_{nls}))-\frac{\chi^2_{0.95,2}}{2} \right\} \tag{3.13} \end{equation}\]

For optimization the Nelder-Meadt algorithm was used implemented in R (R Core Team 2018).

**Possible outcomes of a parameter identifiability analysis** The log-likelihood profile for a given parameter and binding model is obtained by treating this parameter as a constant set to a specified value while re-optimizing the remaining parameter of that model as part of its likelihood function. The corresponding 95% profile log-likelihood confidence interval (red line) is derived using a threshold ($\Delta_\alpha$) in the likelihood according to a critical value of a $\chi^2_{.95,1}$-distribution with 1 $df$. The shape of the profile with respect to the $CI$ is used to judge parameter identifiability. **A:** The log-likelihood profile of an identifialbe parameter contains a unique maximum and a finite $CI$. **B:** A structurally unidentifiable parameter has a flat log-likelihood profile and an infinite $CI$. **C:** A practically unidentifiable parameter shows a unique maximum of the log-likelihood profile and a $CI$ that is only bound on either lower or uppper values of the parameter under study.

Figure 3.2: Possible outcomes of a parameter identifiability analysis The log-likelihood profile for a given parameter and binding model is obtained by treating this parameter as a constant set to a specified value while re-optimizing the remaining parameter of that model as part of its likelihood function. The corresponding 95% profile log-likelihood confidence interval (red line) is derived using a threshold (\(\Delta_\alpha\)) in the likelihood according to a critical value of a \(\chi^2_{.95,1}\)-distribution with 1 \(df\). The shape of the profile with respect to the \(CI\) is used to judge parameter identifiability. A: The log-likelihood profile of an identifialbe parameter contains a unique maximum and a finite \(CI\). B: A structurally unidentifiable parameter has a flat log-likelihood profile and an infinite \(CI\). C: A practically unidentifiable parameter shows a unique maximum of the log-likelihood profile and a \(CI\) that is only bound on either lower or uppper values of the parameter under study.

3.2.4 Model evaluation based on the Akaike Information Criterion

The performance of each model function in describing the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) titration profile was evaluated based on the Akaike Information Criterion (\(AIC\)) as described in Anderson and Burnham (2004). In this work the \(AIC\) corrected for small sample sizes (\(AIC_c\)) was used. Following this, binding models were ranked based on the \(AIC_c\) accoridng to equation (3.14):

\[\begin{equation} AIC_c = n\,log\left(\frac{\sum\limits_{i=1}^n \hat{\epsilon}^2_i}{n}\right)+2k+\frac{2k(k+1)}{n-k-1} \tag{3.14} \end{equation}\]

Where \(n\) = number of data points, \(k\) = number of parameters and \(\hat{\epsilon}^2_i\) = the estimated residuals. Binding models were further ranked by \(AIC_c\) differences (\(\Delta\)) specified in equation (3.15):

\[\begin{equation} \Delta_i=AIC_{c,i}-AIC_{c,min} \tag{3.15} \end{equation}\]

Where \(AIC_{c,min}\) is the \(AIC_c\) of the top-ranked binding model that received the minimum \(AIC_c\) score according to equation (3.14). Based on \(AIC_c\) differences Akaike weights (weights of evidence) for each binding model (\(\Delta\)i) with respect to the whole model set (\(\Delta\)r) were calculated based on equation (3.16):

\[\begin{equation} \omega_i=\frac{exp(-\frac{1}{2}\Delta_i)}{\sum\limits_{r=1}^R exp(-\frac{1}{2}\Delta_r)} \tag{3.16} \end{equation}\]

Based on individual Akaike weights evidence ratios between two binding models were calculated based on equation (3.17):

\[\begin{equation} \frac{\omega_i}{\omega_j} \tag{3.17} \end{equation}\]

3.2.5 Complete k-Fold Cross-validation

The predictive power of model functions was assessed by complete \(k\)-fold cross-validation (CV). For the special case where \(k\) = 4, the number of data points \(n\) in the titration data set is partitioned according to the following relation \[ train \backslash test := \{x \in train | x \notin test \} \] into two disjoint sub-sets called folds. The test-fold contains \(x\) = \(n/k\) = 3 data points and the complementary train-fold contains the remaining data points. Subsequently, nls regression was performed for each binding model on the train-fold and obtained parameter point estimates were then fixed and used to predict the data points in the test-fold. This was done for each possible combination \(C_3(12) = \binom{n}{n/k} = \binom{12}{3} = 220\) of test-folds containing three data points. For each test-fold combination, the predictive power was assessed by computation of a CV-score, in this case the root mean squared error (\(RMSE\)) between the predicted (\(\hat{y}_i\)) and the observed (\(y_i\)) data points according to equation (3.18). Based on that the grand mean of the \(RMSE\) (\(\overline{RMSE}\)) was determined taking into account all possible train-fold combinations. Differences between those grand means with respect to the binding models were further investigated by analysis of variance (ANOVA) described in subsequent chapters (3.2.6).

\[\begin{gather} \overline{RMSE} = \frac{1}{C} \sum\limits_{i=1}^{C} RMSE_i \quad \text{,where}\\ RMSE = \sqrt{\frac{SSE}{m}} \quad \text{,and}\\ SSE = \sum\limits_{j=1}^m (\hat{y}_j-y_j)^2 \notag \tag{3.18} \end{gather}\]

3.2.6 Analysis of Variance (ANOVA)

Binding model dependent differences in the central tendencies of CV-scores obtained in 3.2.5 were further investigated by ANOVA based on a one-way independent and balanced factorial design. For \(m = 5\) category means representing the binding models the multiple regression model containing \(J = m - 1\) predictor variables is given by equation (3.19):

\[\begin{equation} Y_{ij} = \beta_0 + \beta_1X_{1i} + \beta_2X_{2i} + \beta_3X_{3i} + \beta_4X_{4i} + \epsilon_{ij} \tag{3.19} \end{equation}\]

The set of regression coefficients \(\beta_{X_k}\) is associated with a complete set of a priori specified orthogonal linear contrasts \(\Psi_k\) which are linked via the following relation (3.20):

\[\begin{gather} \beta_{X_k} = \frac{\Psi_k}{\sum\limits_{j=1}^m c_{jk}^2} \quad \text{, where} \\ \Psi_k = \sum\limits_{j=1}^m c_{jk} \bar{Y_j} \notag \tag{3.20} \end{gather}\]

where \(\bar{Y_j}\) is the respective group mean and \(\Psi_k\) is the \(k^{th}\) linear contrast (\(k = \{1,\dots, J\}\)) with contrast coefficients \(c_j\) (\(j = \{1,\dots, m\}\)) which are given by the columns of the contrast coefficient matrix \(K\):

\[\begin{equation} K = \begin{bmatrix} -4/5 & 0 & 0 & 0 \\ 1/5 & 1/4 & -1/3 & -1/2 \\ 1/5 & 1/4 & -1/3 & 1/2 \\ 1/5 & 1/4 & 2/3 & 0 \\ 1/5 & -3/4 & 0 & 0 \\ \end{bmatrix} \tag{3.21} \end{equation}\]

In this way, by combining equations (3.20) and (3.21), each regression coefficient specifies a distinct difference of group means leading to the null hypotheses (\(H_{01} - H_{04}\)) expressed in (3.22):

\[\begin{equation} \begin{aligned} H_{01}: \Psi_1 &= 0 \\ \beta_{X_1} &= \left(\frac{\bar{Y}_{bm1to2.deg.add} + \bar{Y}_{bm1to2.deg} + \bar{Y}_{bm1to2.coop.add} + \bar{Y}_{bm1to2.coop}}{4}\right) - \bar{Y}_{bm1to1} = 0 \\ H_{02}: \Psi_2 &= 0 \\ \beta_{X_2} &= \left(\frac{\bar{Y}_{bm1to2.deg.add} + \bar{Y}_{bm1to2.deg} + \bar{Y}_{bm1to2.coop.add}}{3}\right) - \bar{Y}_{bm1to2.coop} = 0 \\ H_{03}: \Psi_3 &= 0 \\ \beta_{X_3} &= \bar{Y}_{bm1to2.coop.add} - \left( \frac{\bar{Y}_{bm1to2.deg.add} + \bar{Y}_{bm1to2.deg}}{2}\right) = 0 \\ H_{04}: \Psi_4 &= 0 \\ \beta_{X_4} &= \frac{\bar{Y}_{bm1to2.deg} - \bar{Y}_{bm1to2.deg.add}}{2} = 0 \end{aligned} \tag{3.22} \end{equation}\]

In terms of a model comparison approach, those null hypotheses are expressed by the following compact models \(C_1 - C_4\) which are compared to the alternative hypothesis (\(H_A\)) specified by the augmented model (\(ModelA\)) given in (3.19):

\[\begin{equation} \begin{aligned} H_{01}: Model C1: Y_i &= \beta_0 + \beta_2X_{2i} + \beta_3X_{3i} + \beta_4X_{4i} + \epsilon_i \\ H_{02}: Model C2: Y_i &= \beta_0 + \beta_1X_{1i} + \beta_3X_{3i} + \beta_4X_{4i} + \epsilon_i \\ H_{03}: Model C3: Y_i &= \beta_0 + \beta_1X_{1i} + \beta_2X_{2i} + \beta_4X_{4i} + \epsilon_i \\ H_{04}: Model C4: Y_i &= \beta_0 + \beta_1X_{1i} + \beta_2X_{2i} + \beta_3X_{3i} + \epsilon_i \\ H_A: Model A: Y_i &= \beta_0 + \beta_1X_{1i} + \beta_2X_{2i} + \beta_3X_{3i} + \beta_4X_{4i} + \epsilon_i \end{aligned} \tag{3.23} \end{equation}\]

In a first instance, residuals were assumed to be normally and independently distributed with constant variance (\(\epsilon_{ij}\overset{iid}{\sim} N(0,\sigma^2)\)). In this case the \(F\)-statistic testing the null hypotheses while comparing the variance of the augmented model with the variances of the various contrasts can be represented by equation (3.24), whereby the notation of Judd, McClelland, and Ryan (2017) was followed:

\[\begin{equation} F = \frac{SSR \times \frac{1}{p_A - p_C}}{SSE_A \times \frac{1}{n - p_A}} = \frac{SSR}{MSE_A} \tag{3.24} \end{equation}\]

Here, \(SSR=SSE_{C_i} - SSE_A\) represents the sum of squares that is reduced when a contrast-coded predcitor is included in the model which is always associated with 1 \(df\) determined by the number of parameter in the augmented (\(p_A\)) and compact (\(p_C\)) model, respectively. It can also be derived in closed form via contrasts (equation (3.25)). Since a complete set of \(m-1\) orthogonal contrasts is used, their individual \(SSR\) contributions (\(SSR_{C_i}\)) sum to \(SSR_A\), the sum of squares reduced when the augmented model is used intstead of the simple grand mean model (\(Y_{ij} = \beta_0 + \epsilon{ij}\)). The term \(SSE_A\) captures the sum of squared error of the augmented model (see also equation (3.26)) which leads to the mean square error \(MSE_A\) when divided by the appropriate \(df\).

\[\begin{gather} SSR = SSE_{C_i} - SSE_A = \frac{\Psi_k^2}{\sum\limits_{j=1}^m \frac{c_{jk}^2}{n_k}} \quad \text{, and} \\ SSR_A = \sum\limits_{i=1}^{k} SSR_{C_i} \notag \tag{3.25} \end{gather}\]

\[\begin{gather} SSE_A = \sum\limits_{m=1}^k \sum\limits_{i=1}^{n_m} (Y_{mi} - \bar{Y_i}_.)^2 \quad \text{, and} \\ \bar{Y_i}_. = \frac{1}{n_m} \sum\limits_{i=1}^{n_m} Y_{mi} \notag \tag{3.26} \end{gather}\]

At a given nominal significance level \(\alpha\) the \(F\)-statistic was then compared to the critical \(F^*_{\alpha;df1,df2}\)-statistic having degrees of freedom \(df1\) and \(df2\) in the numerator and denominator, respectively. Based on the \(F^*_{\alpha;df1,df2}\)-statistic 100(\(1 - \alpha\))% confidence intervals (\(CI_{0.95}\)) were additionally computed for individual contrasts according to equation (3.27):

\[\begin{equation} CI_{0.95}: \hat{\Psi}_k \pm \sqrt{F^*_{\alpha;df1,df2}} \times \sqrt{MSE_A \times \sum\limits_{j=1}^m \frac{c_{jk}^2}{n_k}} \tag{3.27} \end{equation}\]

In a second instance, a robust version of the ANOVA was also conducted specifically addressing the assumption of normally distributed residuals.

3.2.7 Robust ANOVA

In order to verify the outcome of the ANOVA as described in 3.2.6 a robust version was additionally applied. In this case the generic ANOVA was made robust against the violation of the assumption of normally distributed residuals. This was achieved by using trimmed means following the description of Wilcox (Wilcox 2011) and using R code provided by the author therein.

When the amount of trimming is denoted by \(\gamma\) and the sample size by \(n\), then \(g\) = \([\gamma n]\), where \([\gamma n]\) is the value of \(\gamma n\) rounded down to the nearest integer. Then, the trimmed mean \(\bar{Y}_t\) of a random sample \(Y_1, \cdots, Y_n\) is computed by removing, respectively, the \(g\) smallest and largest observations of the sample while averaging the remaining ones according to equation (3.28):

\[\begin{equation} \bar{Y}_t = \frac{Y_{(g+1)} + \cdots + Y_{(n-g)}}{n - 2g} \tag{3.28} \end{equation}\]

where \(Y_i\) are the order statistics \(Y_{(1)} \leq Y_{(2)} \leq \cdots \leq Y_{(n)}\). In this work 20% of the observations were trimmed from each tail. The trimmed mean is used in conjunction with the Winsorized sample variance. During Winsorization the values of the \(g\) smallest observations are set equal to \(Y_{(g+1)}\), while the values of the \(g\) largest observations are set equal to \(Y_{(n−g)}\) according to equation (3.29)

\[\begin{equation} W_i = \begin{cases} Y_{(g+1)}, & \text{if} \quad Y_i \leq Y_{(g+1)}\\ Y_i, & \text{if} \quad Y_{(g+1)} < Y_i < Y_{(n-g)}\\ Y_{(n-g)}, & \text{if} \quad Y_i \geq Y_{(n-g)} \end{cases} \tag{3.29} \end{equation}\]

When \(\bar{W}\) represents the Winsorized sample mean then the corresponding Winsorized variance \(s_w^2\) is given by equation (3.30):

\[\begin{equation} s_w^2 = \frac{1}{n-1} \sum\limits_{i=1}^n (W_i - \bar{W})^2 \tag{3.30} \end{equation}\]

The linear contrast associated with a trimmed mean (\(\bar{Y_t}\)) is now given by:

\[\begin{equation} \Psi_k = \sum\limits_{j=1}^J c_{jk} \bar{Y_{jt}} \tag{3.31} \end{equation}\]

An estimate for the squared standard error of the \(J\) linear contrasts first described by Yuen (Yuen 1974) takes the following form:

\[\begin{gather} s_{yu,k}^2 = \sum\limits_{j=1}^J d_j \quad \text{, where} \\ d_j = \frac{c_{jk}^2(n_j - 1) s_{wj}^2}{h_j(h_j - 1)} \tag{3.32} \end{gather}\]

Here, \(h_j\) represents the effective sample size after trimming. The estimated degrees of freedom \(\hat\nu\) associated with the \(k\)th linear contrast is:

\[\begin{gather} \hat\nu = \frac{(s_{yu,k}^2)^2}{D} \quad \text{, where} \\ D = \sum\limits_{j=1}^J\frac{d_j^2}{h_j -1} \tag{3.33} \end{gather}\]

In terms of testing the null hypotheses (e.g., \(H_0: \Psi_k = 0\)) previously described in (3.22) an appropriate test statistic is now the \(T\) statistic according to equation (3.34):

\[\begin{equation} T_k = \frac{\Psi_k}{\sqrt{s_{yu,k}^2}} \tag{3.34} \end{equation}\]

\(T_k\) was considered significant at the nominal \(\alpha\)-level if \(|T_k| > t_k\), where \(t_k\) is the \(1 - \alpha\) percentage point of a C-variate studentized maximum modulus distribution with estimated degrees of freedom of \(\hat\nu\). In this way the familywise error rate (FWE) of at least one Type I error among all J tests was kept at the nominal \(\alpha\)-level. The corresponding \(100(1 - \alpha) \%\) confidence interval was constructed based on equation (3.35):

\[\begin{equation} CI_{0.95}:\quad \Psi_k \pm t_k \sqrt{s_{yu,k}^2} \tag{3.35} \end{equation}\]

smmcrit<-function(nuhat,C){
  #
  #  Determine the .95 quantile of the C-variate Studentized maximum
  #  modulus distribution using linear interpolation on inverse
  #  degrees of freedom
  #  If C=1, this function returns the .975 quantile of Student's t
  #  distribution.
  #
  if(C-round(C)!=0)stop("The number of contrasts, C, must be an  integer")
  if(C>=29)stop("C must be less than or equal to 28")
  if(C<=0)stop("C must be greater than or equal to 1")
  if(nuhat<2)stop("The degrees of freedom must be greater than or equal to 2")
  if(C==1)smmcrit<-qt(.975,nuhat)
  if(C>=2){
    C<-C-1
    m1<-matrix(0,20,27)
    m1[1,]<-c(5.57,6.34,6.89,7.31,7.65,7.93,8.17,8.83,8.57,
              8.74,8.89,9.03,9.16,9.28,9.39,9.49,9.59, 9.68,
              9.77,9.85,9.92,10.00,10.07,10.13,10.20,10.26,10.32)
    m1[2,]<-c(3.96,4.43,4.76,5.02,5.23,5.41,5.56,5.69,5.81,
              5.92,6.01,6.10,6.18,6.26,6.33,6.39,6.45,6.51,
              6.57,6.62,6.67,6.71,6.76,6.80,6.84,6.88, 6.92)
    m1[3,]<-c(3.38,3.74,4.01,4.20,4.37,4.50,4.62,4.72,4.82,
              4.89,4.97,5.04,5.11,5.17,5.22,5.27,5.32, 5.37,
              5.41,5.45,5.49,5.52,5.56,5.59,5.63,5.66,5.69)
    m1[4,]<-c(3.09,3.39,3.62,3.79,3.93,4.04,4.14,4.23,4.31,
              4.38,4.45,4.51,4.56,4.61,4.66,4.70,4.74,4.78,
              4.82,4.85,4.89,4.92,4.95,4.98,5.00,5.03,5.06)
    m1[5,]<-c(2.92,3.19,3.39,3.54,3.66,3.77,3.86,3.94,4.01,
              4.07,4.13,4.18,4.23,4.28,4.32,4.36,4.39,4.43,
              4.46,4.49,4.52,4.55,4.58,4.60,4.63,4.65,4.68)
    m1[6,]<-c(2.80,3.06,3.24,3.38,3.49,3.59,3.67,3.74,3.80,
              3.86,3.92,3.96,4.01,4.05,4.09,4.13,4.16,4.19,
              4.22,4.25,4.28,4.31,4.33,4.35,4.38,4.39,4.42)
    m1[7,]<-c(2.72,2.96,3.13,3.26,3.36,3.45,3.53,3.60,3.66,
              3.71,3.76,3.81,3.85,3.89,3.93,3.96,3.99, 4.02,
              4.05,4.08,4.10,4.13,4.15,4.18,4.19,4.22,4.24)
    m1[8,]<-c(2.66,2.89,3.05,3.17,3.27,3.36,3.43,3.49,3.55,
              3.60,3.65,3.69,3.73,3.77,3.80,3.84,3.87,3.89,
              3.92,3.95,3.97,3.99,4.02,4.04,4.06,4.08,4.09)
    m1[9,]<-c(2.61,2.83,2.98,3.10,3.19,3.28,3.35,3.41,3.47,
              3.52,3.56,3.60,3.64,3.68,3.71,3.74,3.77,3.79,
              3.82,3.85,3.87,3.89,3.91,3.94,3.95, 3.97,3.99)
    m1[10,]<-c(2.57,2.78,2.93,3.05,3.14,3.22,3.29,3.35,3.40,
               3.45,3.49,3.53,3.57,3.60,3.63,3.66,3.69,3.72,
               3.74,3.77,3.79,3.81,3.83,3.85,3.87,3.89,3.91)
    m1[11,]<-c(2.54,2.75,2.89,3.01,3.09,3.17,3.24,3.29,3.35,
               3.39,3.43,3.47,3.51,3.54,3.57,3.60,3.63,3.65,
               3.68,3.70,3.72,3.74,3.76,3.78,3.80,3.82,3.83)
    m1[12,]<-c(2.49,2.69,2.83,2.94,3.02,3.09,3.16,3.21,3.26,
               3.30,3.34,3.38,3.41,3.45,3.48,3.50,3.53,3.55,
               3.58,3.59,3.62,3.64,3.66,3.68,3.69,3.71,3.73)
    m1[13,]<-c(2.46,2.65,2.78,2.89,2.97,3.04,3.09,3.15,3.19,
               3.24,3.28,3.31,3.35,3.38,3.40,3.43,3.46,3.48,
               3.50,3.52,3.54,3.56,3.58,3.59,3.61,3.63,3.64)
    m1[14,]<-c(2.43,2.62,2.75,2.85,2.93,2.99,3.05,3.11,3.15,
               3.19,3.23,3.26,3.29,3.32,3.35,3.38,3.40,3.42,
               3.44,3.46,3.48,3.50,3.52,3.54,3.55,3.57,3.58)
    m1[15,]<-c(2.41,2.59,2.72,2.82,2.89,2.96,3.02,3.07,3.11,
               3.15,3.19,3.22,3.25,3.28,3.31,3.33,3.36,3.38,
               3.39,3.42,3.44,3.46,3.47,3.49,3.50,3.52,3.53)
    m1[16,]<-c(2.38,2.56,2.68,2.77,2.85,2.91,2.97,3.02,3.06,
               3.09,3.13,3.16,3.19,3.22,3.25,3.27,3.29,3.31,
               3.33,3.35,3.37,3.39,3.40,3.42,3.43,3.45,3.46)
    m1[17,]<-c(2.35,2.52,2.64,2.73,2.80,2.87,2.92,2.96,3.01,
               3.04,3.07,3.11,3.13,3.16,3.18,3.21,3.23,3.25,
               3.27,3.29,3.30,3.32,3.33,3.35,3.36,3.37,3.39)
    m1[18,]<-c(2.32,2.49,2.60,2.69,2.76,2.82,2.87,2.91,2.95,
               2.99,3.02,3.05,3.08,3.09,3.12,3.14,3.17, 3.18,
               3.20,3.22,3.24,3.25,3.27,3.28,3.29,3.31,3.32)
    m1[19,]<-c(2.29,2.45,2.56,2.65,2.72,2.77,2.82,2.86,2.90,
               2.93,2.96,2.99,3.02,3.04,3.06,3.08,3.10, 3.12,
               3.14,3.16,3.17,3.19,3.20,3.21,3.23,3.24,3.25)
    m1[20,]<-c(2.24,2.39,2.49,2.57,2.63,2.68,2.73,2.77,2.79,
               2.83,2.86,2.88,2.91,2.93,2.95,2.97,2.98, 3.01,
               3.02,3.03,3.04,3.06,3.07,3.08,3.09,3.11,3.12)
    if(nuhat>=200)smmcrit<-m1[20,C]
    if(nuhat<200){
      nu<-c(2,3,4,5,6,7,8,9,10,11,12,14,16,18,20,24,30,40,60,200)
      temp<-abs(nu-nuhat)
      find<-order(temp)
      if(temp[find[1]]==0)smmcrit<-m1[find[1],C]
      if(temp[find[1]]!=0){
        if(nuhat>nu[find[1]]){
          smmcrit<-m1[find[1],C]-
            (1/nu[find[1]]-1/nuhat)*(m1[find[1],C]-m1[find[1]+1,C])/
            (1/nu[find[1]]-1/nu[find[1]+1])
        }
        if(nuhat<nu[find[1]]){
          smmcrit<-m1[find[1]-1,C]-
            (1/nu[find[1]-1]-1/nuhat)*(m1[find[1]-1,C]-m1[find[1],C])/
            (1/nu[find[1]-1]-1/nu[find[1]])
        }
      }}
  }
  smmcrit
}

qsmm<-function(q, r, nu) {
  #r=number of comparisons
  if (!is.finite(nu)) 
    return(qnorm(1 - 0.5 * (1 - q^(1/r))))
  res = uniroot(function(c, r, nu, q) {
    psmm(c, r = r, nu = nu) - q
  },
  c(0, 100), r = r, nu = nu, q = q)
  res$root
}


psmm = function(x, r, nu) {
  res = integrate(psmm.x, 0, Inf, c = x, r = r, nu = nu)
  res$value
}


psmm.x=function(x, c, r, nu) {
  snu = sqrt(nu)
  sx = snu * x
  lgx = log(snu) - lgamma(nu/2) + (1 - nu/2) * log(2) + 
    (nu - 1) * log(sx) + (-sx^2/2)
  exp(r * log(2 * pnorm(c * x) - 1) + lgx)
}
lincon<-function(x,con=0,tr=.2,alpha=.05,pr=FALSE){
  #
  #  A heteroscedastic test of d linear contrasts using trimmed means.
  #
  #  This version uses an improved method for computing the quantiles of a
  #  Studentized maximum modulus distriburtion
  # 
  #  The data are assumed to be stored in $x$ in list mode, a matrix
  #  or a data frame. If in list mode,
  #  length(x) is assumed to correspond to the total number of groups.
  #  It is assumed all groups are independent.
  #
  #  con is a J by d matrix containing the contrast coefficients that are used.
  #  If con is not specified, all pairwise comparisons are made.
  #
  #  Missing values are automatically removed.
  #
  #  pr=FALSE included to avoid errors using an earlier version of this function when
  #   dealing with two-way and higher designs
  #
  #  Adjusted p-values are based on the Studentized maximum modulus distribution with the 
  #  goal of controlling FWE
  #
  #  To apply the Kaiser-Bowden method, use the function kbcon
  #
  if(tr==.5) stop('Use the R function medpb to compare medians')
  if(is.data.frame(x)) x = as.matrix(x)
  if(is.matrix(x)) x <- listm(x)
  if(!is.list(x))stop('Data must be stored in a matrix or in list mode.')
  con<-as.matrix(con)
  J<-length(x)
  sam=NA
  h<-vector('numeric',J)
  w<-vector('numeric',J)
  xbar<-vector('numeric',J)
  for(j in 1:J){
    xx<-!is.na(x[[j]])
    val<-x[[j]]
    x[[j]]<-val[xx]  # Remove missing values
    sam[j]=length(x[[j]])
    h[j]<-length(x[[j]])-2*floor(tr*length(x[[j]]))
    # h is the number of observations in the jth group after trimming.
    w[j]<-((length(x[[j]])-1)*winvar(x[[j]],tr))/(h[j]*(h[j]-1))
    xbar[j]<-mean(x[[j]],tr)
  }
  if(sum(con^2)==0){
    CC<-(J^2-J)/2
    psihat<-matrix(0,CC,9)
    dimnames(psihat)<-list(NULL,c('Group','Group','psihat','ci.lower','ci.upper',
                                  'p.value','Est.1','Est.2','adj.p.value'))
    test<-matrix(NA,CC,6)
    dimnames(test)<-list(NULL,c('Group','Group','test','crit','se','df'))
    jcom<-0
    for (j in 1:J){
      for (k in 1:J){
        if (j < k){
          jcom<-jcom+1
          test[jcom,3]<-abs(xbar[j]-xbar[k])/sqrt(w[j]+w[k])
          sejk<-sqrt(w[j]+w[k])
          test[jcom,5]<-sejk
          psihat[jcom,1]<-j
          psihat[jcom,2]<-k
          test[jcom,1]<-j
          test[jcom,2]<-k
          psihat[jcom,3]<-(xbar[j]-xbar[k])
          df<-(w[j]+w[k])^2/(w[j]^2/(h[j]-1)+w[k]^2/(h[k]-1))
          test[jcom,6]<-df
          psihat[jcom,6]<-2*(1-pt(test[jcom,3],df))
          psihat[jcom,7]=xbar[j]
          psihat[jcom,8]=xbar[k]
          crit=qsmm(1-alpha,CC,df)
          test[jcom,4]<-crit
          psihat[jcom,4]<-(xbar[j]-xbar[k])-crit*sejk
          psihat[jcom,5]<-(xbar[j]-xbar[k])+crit*sejk
          psihat[jcom,9]=1-psmm(test[jcom,3],CC,df)
        }}}}
  if(sum(con^2)>0){
    if(nrow(con)!=length(x)){
      stop('The number of groups does not match the number of contrast coefficients.')
    }
    CC=ncol(con)
    psihat<-matrix(0,ncol(con),6)
    dimnames(psihat)<-list(NULL,c('con.num','psihat','ci.lower','ci.upper',
                                  'p.value','adj.p.value'))
    test<-matrix(0,ncol(con),5)
    dimnames(test)<-list(NULL,c('con.num','test','crit','se','df'))
    df<-0
    for (d in 1:ncol(con)){
      psihat[d,1]<-d
      psihat[d,2]<-sum(con[,d]*xbar)
      sejk<-sqrt(sum(con[,d]^2*w))
      test[d,1]<-d
      test[d,2]<-sum(con[,d]*xbar)/sejk
      df<-(sum(con[,d]^2*w))^2/sum(con[,d]^4*w^2/(h-1))
      crit=qsmm(1-alpha,CC,df)
      test[d,3]<-crit
      test[d,4]<-sejk
      test[d,5]<-df
      psihat[d,3]<-psihat[d,2]-crit*sejk
      psihat[d,4]<-psihat[d,2]+crit*sejk
      psihat[d,5]<-2*(1-pt(abs(test[d,2]),df))
      psihat[d,6]=1-psmm(abs(test[d,2]),CC,df)
    }
  }
  list(n=sam,test=test,psihat=psihat)
}

4 Results

4.1 Modeling the ThiITm•tRNAPhe interaction

To figure out the stoichiometry that underlies complexes formed between ThiITm and full-length tRNA substrates tryptophan fluorescence spectroscopic binding studies were carried out using the example of the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) interaction. In this context, unveiling the true underlying stoichiometry is inherently linked to finding the most appropriate binding equilibrium describing the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) binding profile. In the following sections an analytic workflow is going to be presented on the basis of a multi-model approach that is capable of objectively finding the statistically best approximating binding model in a set of binding models with varying complexity.

4.1.1 Fluorescence Spectroscopic Binding Studies

ThiITm is characterized by an intrinsic tryptophan fluorescence signal which is quenched in a concentration dependent manner in the presence of tRNA\(\mathrm{^{Phe}_{Bs}}\). This offered the possibility to quantitatively examine ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) complex formation based on a fluorescence spectroscopic titration experiment.

##    replica1 replica2 replica3    avg_fi        trna      que_fi      stdev
## 1    615018   879610   735246  743291.3 30.00000000  0.73054270 0.04908259
## 2   1042800  1127699  1026932 1065810.3 15.00000000  0.61284602 0.03482094
## 3   1525162  1700574  1448364 1558033.3  7.50000000  0.43351491 0.06918158
## 4   1928006  2083840  1806328 1939391.3  3.75000000  0.29479920 0.08162473
## 5   2346354  2283607  1985989 2205316.7  1.87500000  0.19694023 0.11335427
## 6   2514418  2313606  2468219 2432081.0  0.93750000  0.11737962 0.06200229
## 7   2714597  2550381  2565743 2610240.3  0.46875000  0.05183969 0.07894547
## 8   2781468  2665702  2442915 2630028.3  0.23437500  0.04308120 0.11514445
## 9   2753056  2707707  2484951 2648571.3  0.11718750  0.03673332 0.10517550
## 10  2740884  2627076  2543639 2637199.7  0.05859375  0.04156523 0.08830984
## 11  2759639  2739218  2897354 2798737.0  0.02929688 -0.01469239 0.03106489
## 12  2583678  2570650  2857770 2670699.3  0.01464844  0.03281524 0.01490623

Following this, 5 µM of ThiITm were titrated with increasing amounts of tRNA\(\mathrm{^{Phe}_{Bs}}\) ranging from 0.03 μM to 60 μM (3.1.1) and the resulting differences in raw tryptophan fluorescence intensities (replica1 - replica3) were recorded in triplicates (figure 4.1 A). Fluorescence intensities were then averaged across replicates (avg_fi) and for each of the 12 data points the change in fluorescence intensity was subsequently transformed to the proportional signal quenching (que_fi) according to equation (3.1). Binding curves obtained in this way are depicted in figure 4.1B together with standard deviations (stdev) of the mean FI signal highlighted as error bars. The scatter plots exhibit a typical hyperbolic shape as a result of an increasing tryptophan fluorescence quenching with increasing concentrations of tRNA\(\mathrm{^{Phe}_{Bs}}\) probably converging against a saturation level around 0.8. Computed error bars illustrate an increasing measurement uncertainty in the data points with an increasing dilution factor of the tRNA\(\mathrm{^{Phe}_{Bs}}\) concentration. The ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) binding curve was the basis for non-linear least squares (nls) regression analysis of different binding models related to underlying 1:1 (ThiITm:tRNA\(\mathrm{^{Phe}_{Bs}}\)) or 1:2 binding equilibria (3.2.1).

In total a set of five model functions was comparatively analyzed composed of four two-site binding models and one one-site binding model. Best fit regression curves of each model type are shown as line plots in figure 4.1B and corresponding point estimates of the parameter are listed in table 4.1. In order to compare the observed data with the predictions of the binding models the root mean squared error (\(RMSE\)) of the residuals was used as a first model performance metric not taking into account the number of model parameters. All regression curves elevate in a hyperbolic progression in accordance with the titration profile with only subtle differences when phenomenologically judged. The exact underlying fitting procedure thereby was specific for the type of binding models. The set of two-site binding models (bm1to2.*) based on the Adair-Klotz polynomial (equations (3.2), (3.4), (3.5), and (3.6)) describe sequential binding to independent binding sites and are all dependent on free ligand concentrations \([G]\) at equilibrium. Therefore, this quantity was calculated from total ligand concentrations \([G]_0\) used in the titration experiment according to the cubic polynomial described by equation (3.7) which represented an integral part of the least squares fitting routine further detailed in 3.2.2. As a consequence, each binding model of this type is on a unique ligand concentration scale that is shifted in relation to those of the remaining Adair-Klotz type binding models and in relation to the total ligand concentration based on which the one-site binding model was fitted.

The cooperative/non-additive model bm1to2.coop(equation (3.2)) represents the global model of the whole set containing four parameters and the least degrees of freedom (\(df\) = 8). The \(RMSE\) of the residuals for this model amounted to 1.79 \(\times\) 10-2 and best fit values for the macroscopic binding constants \(K_1\) and \(K_2\) were estimated to be 1.38 \(\times\) 106 μM-1 and 1.17 \(\times\) 10-2 μM-1, respectively. Therefore, bm1to2.coop describes tRNA\(\mathrm{^{Phe}_{Bs}}\) binding to ThiITm as an extreme case of positive cooperativity with an unfeasible high interaction parameter \(\alpha\) of 11.8 \(\times\) 106 as determined by equation (3.3). On the basis of the global model a subset of re-parameterized nested 1:2 binding models could be derived by either merging the proportionality factors (bm1to2.coop.add) or the macroscopic binding constants (bm1to2.deg) or both (bm1to2.deg.add). Among these a fit of the cooperative/additive model bm1to2.coop.add resulted in parameter estimations for \(K_1\) and \(K_2\) of 2.69 \(\times\) 10-1 μM-1 and 3.66 \(\times\) 10-2 μM-1 ,respectively, so that the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) interaction would be positively cooperative with an \(\alpha\) value of 7.35. With respect to the \(RMSE\) of the residuals the regression curve of bm1to2.coop.add matched the titration data second best among all models considered (\(RMSE\) = 2.05 \(\times\) 10-2). The degenerate/non-additive model bm1to2.deg describes binding to independent and degenerate binding sites with an overall stability constant \(\beta_{12}\) of 4.58 \(\times\) 10-1 μM-2. The closely related degenerate/additive model bm1to2.deg.add got the same value for its stability constant. Also, both models share the same value for the proportionality factor \(F_{\Delta HG}\) of the first binding step amounting to 8,61 \(\times\) 10-2 μM-1. The least squares regression analysis was completed by the one-site binding model bm1to1 whose fit revealed the highest \(RMSE\) of the residuals in the whole set (2.66 \(\times\) 10-2) and which yielded a \(K_1\) point estimate of 2.7 \(\times\) 10-1 μM-1.

Though each model is adequately descriptive and well capturing the progress of the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) binding curve, the increasing uncertainty in the data points with decreasing tRNA\(\mathrm{^{Phe}_{Bs}}\) concentration as well as the extreme parameter estimates obtained by the global bm1to2.coop model necessitated further investigations prior to final model evaluation steps. In a first step this was achieved by a so-called parameter identifiability analysis based on maximum likelihood profiling of model parameter.

**Intrinsic tryptophan fluorescence quenching asssay**. Five µM of ThiI<sub>Tm</sub> (20 mM Tris/HCl pH 7.5, 300 mM NaCl, 2 mM DTT, 5 mM MGCl2, 5% glycerol) were titrated with increasing amounts of tRNA$^{Phe}_{Bs}$ (0.03 - 60 µM) and the decrease of the intrinsic tryptophan fluoresence of ThiI<sub>Tm</sub> was monitored spectrophotometrically. The titration profile was the basis for nls regression analysis. Different two-site and one-site binding models were fitted all capturing the hyperbolic progress of the binding profile with only subtle differences. **A:** Raw data. Shown is the decrease in the tryptophan fluorescence signal for each of the replicate measurements. The mean FI signal was then transformed to reflect the proportional FI signal quenching. **B:** Two-site binding models were fitted to unique sets of free ligand concentrations $[G]$ individually calculated for each binding model. The remaining one-site binding model was fitted based on total ligand concentration $[G]_0$ as specified in the titration experiment.

Figure 4.1: Intrinsic tryptophan fluorescence quenching asssay. Five µM of ThiITm (20 mM Tris/HCl pH 7.5, 300 mM NaCl, 2 mM DTT, 5 mM MGCl2, 5% glycerol) were titrated with increasing amounts of tRNA\(^{Phe}_{Bs}\) (0.03 - 60 µM) and the decrease of the intrinsic tryptophan fluoresence of ThiITm was monitored spectrophotometrically. The titration profile was the basis for nls regression analysis. Different two-site and one-site binding models were fitted all capturing the hyperbolic progress of the binding profile with only subtle differences. A: Raw data. Shown is the decrease in the tryptophan fluorescence signal for each of the replicate measurements. The mean FI signal was then transformed to reflect the proportional FI signal quenching. B: Two-site binding models were fitted to unique sets of free ligand concentrations \([G]\) individually calculated for each binding model. The remaining one-site binding model was fitted based on total ligand concentration \([G]_0\) as specified in the titration experiment.

Table 4.1: Non-linear least squares fit results. Different two-site and one-site binding models were fit to the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) binding profile obtained from the tryptophan fluorescence quenching assay (figure 23). Listed are best fit point estimates of model-type specific parameter obtained from the regression analysis: \(df\) = degrees of freedom, \(RMSE\) = root mean squared error, \(F_{\Delta HG}\)/\(F_{\Delta HG2}\) = proportionality factors, \(K_1\)/\(K_2\) = macroscopic binding constants, and \(\beta_{12}\) = overall stability constant.
\(Model\) \(df\) \(RMSE\) \(F_{\Delta HG} (µM^{-1})\) \(F_{\Delta HG2} (µM^{-1})\) \(\beta_{12} (µM^{-2})\) \(K_1 (µM^{-1})\) \(K_2 (µM^{-1})\)
bm1to1 10 0.0266 0.16500 2.70e-01
bm1to2.deg.add 10 0.0206 0.08610 0.458
bm1to2.deg 9 0.0217 0.08610 0.172 0.458
bm1to2.coop.add 9 0.0205 0.10200 2.69e-01 0.0366
bm1to2.coop 8 0.0179 0.00362 0.187 1.38e+06 0.1170

4.1.2 Parameter Identifiability Analysis

####..........generic functions and objects....................####

prof.seq <- function(lower, upper, len = 100){
  seq(lower, upper, length = len)
}

single.prof.plot <- function(model.list, mle.obj, theta){
  prof.tbl <- tibble(loglik = unlist(model.list$lik.profile[theta]),
                     theta = unlist(model.list$profile.vals[theta])
  )
  ci95 <- mle.obj$value - (qchisq(0.95,1)/2)
  prof.tbl <- prof.tbl %>% 
    mutate(ci_bound = loglik >= ci95)
  
  ggplot(prof.tbl, aes(y = loglik, x = theta, color = ci_bound))+
    geom_segment(aes(x = theta, 
                     y = 0, 
                     xend = theta, 
                     yend = loglik))+
    geom_point()+
    scale_colour_manual(values = c("grey60", "black"))+
    theme(legend.position = "none")+
    geom_vline(xintercept = mle.obj$par[theta],
               color = "#0033FF", size = 1.1, lty = 2)+
    geom_hline(yintercept = ci95,
               color = "#CC0000", size = 1.1)+
    xlab(names(
      model.list$profile.vals[which(names(model.list$profile.vals) == theta)])
    )
}

contour.plot <- function(model.list, mle.obj, theta){
  conf95 <- mle.obj$value - qchisq(0.95,2)/2
  grid.vals <- data.frame(expand.grid(
    unlist(model.list$profile.vals[theta[1]]), 
    unlist(model.list$profile.vals[theta[2]])))
  prof.vals <- data.frame(model.list$lik.profile[paste(theta[1], 
                                                       theta[2], 
                                                       sep = ".")])
  names(prof.vals) <- "loglik"
  prof.contour <- cbind(grid.vals, prof.vals)
  
  ggplot(prof.contour, aes(y = Var1, 
                           x = Var2, 
                           z = loglik))+
    theme_minimal()+
    geom_tile(aes(fill = loglik))+
    stat_contour(color = "black")+
    geom_point(aes(y = mle.obj$par[theta[1]],
                   x = mle.obj$par[theta[2]]),
               color = "#0033FF")+
    scale_fill_continuous(low="grey70",high="grey15")+
    stat_contour(breaks=c(conf95), color = "#CC0000")+
    ylab(theta[1])+
    xlab(theta[2])
}
profile.vals <- list()
lik.profile <- list()

bm1to2.deg.add_nls$nls.coef <- c(coef(bm1to2.deg.add_nls$bm1to2.deg.add),
                                 sig = sigma(bm1to2.deg.add_nls$bm1to2.deg.add))
bm1to2.deg_nls$nls.coef <- c(coef(bm1to2.deg_nls$bm1to2.deg),
                             sig = sigma(bm1to2.deg_nls$bm1to2.deg))
bm1to2.coop.add_nls$nls.coef <- c(coef(bm1to2.coop.add_nls$bm1to2.coop.add),
                                  sig = sigma(bm1to2.coop.add_nls$bm1to2.coop.add))
bm1to2.coop.add_nls$nls.coef <- c(coef(bm1to2.coop.add_nls$bm1to2.coop.add),
                                  sig = sigma(bm1to2.coop.add_nls$bm1to2.coop.add))

bm1to2.coop_nls$nls.coef <- c(coef(bm1to2.coop_nls$bm1to2.coop),
                              sig = sigma(bm1to2.coop_nls$bm1to2.coop))

mytheme_axes_box =
  theme_bw()+
  theme(axis.ticks.length=unit(.07, "cm"))+
  theme(axis.ticks = element_line(colour = "black", size = 1))+
  theme(axis.text.y = element_text(color="black",size = 10))+
  theme(axis.text.x = element_text(color="black",size = 10))+
  theme(axis.title.y = element_text(size = 12,face = "bold"))+
  theme(axis.title.x = element_text(size = 12,face = "bold"))+
  theme(panel.grid.minor=element_blank(),
        panel.grid.major=element_blank(),
        panel.border = element_rect(size = 1, color = "black"))

The purpose of identifiability analysis is to detect whether the point estimates obtained via nls fitting are uniquely determined given a specific model type and measurement set up. In this work parameter identifiability was investigated via maximum likelihood profiling in a qualitative way. During this process a model parameter under study is fixed at specified values around its best fit nls point estimate and all other model parameter are re-optimized via maximum likelihood estimation based on the underlying likelihood function of that particular model. In this context a paramter is identifiable when its maximum likelihood profile has a unique maximum and a finite confidence interval. Besides being identifiable a parameter might also be found to be either practically or structurally non-identifiable as schematically represented in figure 3.2.

As a first analysis step the logarithmic likelihood functions (log-likelihood) were built for each model as detailed in 3.2.3. Due to the fact that the nls regression analysis was carried out under the assumption of normally distributed residuals with constant variance the respective log-likelihood functions were based on a Gaussian probability density function as specified by equation (3.9). This basic function was then used for maximum likelihood profiling as specified by equation (3.10) which is shown here, respectively, for selected parameter of the cooperative bm1to2.coop and bm1to2.coop.add models as well as for the degenerative bm1to2.deg.add model (figure 4.2). The corresponding log-likelihood and profile log-likelihood functions were implemented in R via the ll.1to2.deg.add function in case of degenerative models and via the ll.1to2.coop function in case of cooperative models (code chunk 2). Regarding the bm1to2.coop model the best fit nls point estimate of the first proportionality factor \(F_{\Delta HG}\) amounted to 3.62 \(\times\) 10-3 μM-1. During the profiling routine \(F_{\Delta HG}\) was iteratively varied around this parameter value between 1 \(\times\) 10-3 μM-1 and 1.7 \(\times\) 10-2 μM-1 while the remaining parameter of the bm1to2.coop model were re-optimized based on the respective log-likelihood function (code chunk 3). The resulting profile of \(F_{\Delta HG}\) (figure 4.2B) contains unique maximized log-likelihood values for each specified parameter value of \(F_{\Delta HG}\) shown as points. The red line represents the threshold of a 0.95 percentile confidence interval (\(CI_{0.95}\)) at 33.7 determined based on an underlying \(\chi^2\)-distribution with one degree of freedom and which was calculated based on equation (3.11). The curve progression of the log-likelihood profile with respect to the \(CI_{0.95}\) was used as a qualitative criterion to judge parameter identifiability. The profile of \(F_{\Delta HG}\) exhibits a smooth parabolic curve typical for an identifiable parameter. It crosses the lower as well as the upper border of the \(CI_{0.95}\) and reaches an expected unique maximum at 33.7 when \(F_{\Delta HG}\) is fixed at its best fit nls point estimate (dashed blue line). In case of the macroscopic binding constant \(K_1\) the analysis revealed a non-identifiability of this parameter (figure 4.2D). The computed maximum likelihood profile is flat over a wide range of the scanned parameter space between 1 \(\times\) 106 μM-1 and 5 \(\times\) 106 μM-1 and it does not possess a finite \(CI_{0.95}\) neither in the scanned range nor beyond. For every specified parameter value of \(K_1\) including that of its best fit nls point estimate (1.38 \(\times\) 106 μM-1) the same maximum of the re-optimized log-likelihood function is found so that the likelihood profile does not reveal a unique maximum. A parameter producing such a flat log-likelihood profile is termed a structurally non-identifiable parameter. To study the impact of a structurally non-identifiable parameter on an identifiable parameter within the same binding model the concept of maximum likelihood profiling was extended and simultaneously applied to both parameter, \(F_{\Delta HG}\) and \(K_1\), within their respective parameter ranges (equation (3.12)). The resulting joint maximum likelihood profile of the two-dimensional parameter space (figure 4.2F) shows parallel likelihood contours with respect to \(F_{\Delta HG}\) and the respective confidence region (equation (3.13)) is reduced to a line that could be infinitely extended. Therefore, changing the value of \(K_1\) can be entirely compensated for by \(F_{\Delta HG}\) during the maximum likelihood re-optimization process. The other type of non-identifiability investigated in this work can be found in one of the parameter of thebm1to2.coop.add model. As highlighted in figure 4.2G, \(K_2\) of the bm1to2.coop.add model can be considered a practically non-identifiable parameter since the \(CI_{0.95}\) around its best fit nls point estimate at 3.2 \(\times\) 10-2 µm-1 does not reveal a lower bound of this parameter in the physically allowed range above zero. When the degrees of freedom are further increased all parameters become identifiable. As exemplarily shown in the left panel of figure 4.2 (A, C, E) for parameters \(F_{\Delta HG}\) and \(\beta_{12}\) of the bm1to2.deg.add model, both log-likelihood profiles peak at the respective parameter values of the best fit nls point estimates at 8.61 \(\times\) 10-2 μM-1 and 4.58 \(\times\) 10-1 μM-2, respectively, and have finite confidence intervals. As a result there exists a closed confidence region in the joint profile log-likelihood space of these parameters (E, red circle).

ll.1to2.deg.add <- function(param, obs, g.curr, profile.par = NULL, ...){
  
  g0 <- g.curr
  if (is.null(profile.par)){
    Fpred <- bm1to2.deg.add(param["fhg"], param["k1"], g0)
  }  else if (all(profile.par == "fhg")){
    Fpred <- bm1to2.deg.add(..., param["k1"], g0)
  } else if (all(profile.par == "k1")){
    Fpred <- bm1to2.deg.add(param["fhg"], ..., g0)
  } else if (all(profile.par == c("fhg", "k1"))){
    Fpred <- bm1to2.deg.add(..., g0)
  }
  
  lnlik <- sum(dnorm(x = obs$que_fi, 
                     mean = Fpred, 
                     sd = param["sig"], 
                     log = TRUE))
  return(lnlik)
}

ll.1to2.coop.add <- function(param, obs, g.curr, profile.par = NULL, ...){
  g0 <- g.curr
  if (is.null(profile.par)){
    Fpred <- bm1to2.coop.add(param["fhg"], param["k1"], param["k2"], g0)
  } else if (all(profile.par == "fhg")){
    Fpred <- bm1to2.coop.add(..., param["k1"], param["k2"], g0)
  } else if (all(profile.par == "k1")){
    Fpred <- bm1to2.coop.add(param["fhg"], ..., param["k2"], g0)
  } else if (all(profile.par == "k2")){
    Fpred <- bm1to2.coop.add(param["fhg"], param["k1"], ..., g0)
  } else if (all(profile.par == c("fhg", "k1"))){
    Fpred <- bm1to2.coop.add(..., param["k2"], g0)
  } else if (all(profile.par == c("fhg", "k2"))){
    Fpred <- bm1to2.coop.add(..., param["k1"], g0)
  } else if (all(profile.par == c("k1", "k2"))){
    Fpred <- bm1to2.coop.add(param["fhg"], ..., g0)
  }
  
  lnlik <- sum(dnorm(x = obs$que_fi, 
                     mean = Fpred, 
                     sd = param["sig"], 
                     log = TRUE))
  return(lnlik)
}

ll.1to2.coop <- function(param, obs, g.curr, profile.par = NULL, ...){
  g0 <- g.curr
  if (is.null(profile.par)){
    Fpred <- bm1to2.coop(param["fhg"], param["fhg2"], param["k1"], 
                         param["k2"], g0)
  } else if (all(profile.par == "fhg")){
    Fpred <- bm1to2.coop(..., param["fhg2"], param["k1"], param["k2"], g0)
  } else if (all(profile.par == "fhg2")){
    Fpred <- bm1to2.coop(param["fhg"], ..., param["k1"], param["k2"], g0)
  } else if (all(profile.par == "k1")){
    Fpred <- bm1to2.coop(param["fhg"], param["fhg2"], ..., param["k2"], g0)
  } else if (all(profile.par == "k2")){
    Fpred <- bm1to2.coop(param["fhg"], param["fhg2"], param["k1"], ..., g0)
  } else if (all(profile.par == c("fhg", "fhg2"))){
    Fpred <- bm1to2.coop(..., param["k1"], param["k2"], g0)
  } else if (all(profile.par == c("fhg", "k1"))){
    Fpred <- bm1to2.coop(..., param["fhg2"], param["k2"], g0)
  } else if (all(profile.par == c("fhg", "k2"))){
    Fpred <- bm1to2.coop(..., param["fhg2"], param["k1"], g0)
  } else if (all(profile.par == c("fhg2", "k1"))){
    Fpred <- bm1to2.coop(param["fhg"], ..., param["k2"], g0)
  } else if (all(profile.par == c("fhg2", "k2"))){
    Fpred <- bm1to2.coop(param["fhg"], ..., param["k1"], g0)
  } else if (all(profile.par == c("k1", "k2"))){
    Fpred <- bm1to2.coop(param["fhg"], param["fhg2"], ..., g0)
  }
  
  lnlik <- sum(dnorm(x = obs$que_fi, 
                     mean = Fpred, 
                     sd = param["sig"], 
                     log = TRUE))
  return(lnlik)
}

mle.deg.add <- optim(bm1to2.deg.add_nls$nls.coef,
                     ll.1to2.deg.add,
                     obs = fi,
                     g.curr = bm1to2.deg.add_nls$g,
                     control = list(fnscale=-1)
)

mle.coop.add <- optim(bm1to2.coop.add_nls$nls.coef,
                      ll.1to2.coop.add,
                      obs = fi,
                      g.curr = bm1to2.coop.add_nls$g,
                      control = list(fnscale=-1))

mle.coop <- optim(bm1to2.coop_nls$nls.coef,
                  ll.1to2.coop,
                  obs = fi,
                  g.curr = bm1to2.coop_nls$g,
                  control = list(fnscale=-1))
####..........Single param profiles bm1to2.deg.add.............####
####k1#####
bm1to2.deg.add_nls$profile.vals$k1 <- prof.seq(0.05, 2)

bm1to2.deg.add_nls$lik.profile$k1 <- sapply(
  bm1to2.deg.add_nls$profile.vals$k1, function(k1) 
    optim(bm1to2.deg.add_nls$nls.coef,
          ll.1to2.deg.add,
          obs = fi,
          g.curr = bm1to2.deg.add_nls$g,
          profile.par = "k1",
          k1 = k1,
          method = "L-BFGS-B",
          lower = c(0.001,0.001),
          control=list(fnscale = -1))$val)

####fhg#####################################################

bm1to2.deg.add_nls$profile.vals$fhg <- prof.seq(0.01, 0.2)

bm1to2.deg.add_nls$lik.profile$fhg <- sapply(
  bm1to2.deg.add_nls$profile.vals$fhg, function(fhg) 
    optim(bm1to2.deg.add_nls$nls.coef,
          ll.1to2.deg.add,
          obs = fi,
          g.curr = bm1to2.deg.add_nls$g,
          profile.par = "fhg",
          fhg = fhg,
          method = "L-BFGS-B",
          lower = c(0.001,0.001),
          control=list(fnscale = -1))$val)

####............Single param profiles bm1to2.coop.add..........####
####k2####
bm1to2.coop.add_nls$profile.vals$k2 <- prof.seq(0.001, 0.12)

bm1to2.coop.add_nls$lik.profile$k2 <- sapply(
  bm1to2.coop.add_nls$profile.vals$k2, function(k2) 
    optim(bm1to2.coop.add_nls$nls.coef,
          ll.1to2.coop.add,
          obs = fi,
          g.curr = bm1to2.coop.add_nls$g,
          profile.par = "k2",
          k2 = k2,
          control=list(fnscale = -1))$val)
####............Single param profiles bm1to2.coop.............####
####fhg####
bm1to2.coop_nls$profile.vals$fhg <- prof.seq(0.001, 0.017)
bm1to2.coop_nls$lik.profile$fhg <- sapply(
  bm1to2.coop_nls$profile.vals$fhg, function(fhg) 
    optim(bm1to2.coop_nls$nls.coef,
          ll.1to2.coop,
          obs = fi,
          g.curr = bm1to2.coop_nls$g,
          profile.par = "fhg",
          fhg = fhg,
          method = "L-BFGS-B",
          lower = c(0.001,0.001),
          control=list(fnscale = -1))$val)

####K1
bm1to2.coop_nls$profile.vals$k1 <- prof.seq(1000000, 5000000)
bm1to2.coop_nls$lik.profile$k1 <- sapply(
  bm1to2.coop_nls$profile.vals$k1, function(k1) 
    optim(bm1to2.coop_nls$nls.coef,
          ll.1to2.coop,
          obs = fi,
          g.curr = bm1to2.coop_nls$g,
          profile.par = "k1",
          k1 = k1,
          method = "L-BFGS-B",
          lower = c(0.001,0.001),
          control=list(fnscale = -1))$val)

####.............profile contour bm1to2.deg.add...............#####
####fhg.k1####
grid.deg.add.fhg.k1 <- data.frame(
  t(expand.grid(unlist(
    bm1to2.deg.add_nls$profile.vals$fhg), 
    unlist(bm1to2.deg.add_nls$profile.vals$k1))))

bm1to2.deg.add_nls$lik.profile$fhg.k1 <- sapply(
  grid.deg.add.fhg.k1, function(v) 
    optim(bm1to2.deg.add_nls$nls.coef,
          ll.1to2.deg.add,
          obs = fi,
          g.curr = bm1to2.deg.add_nls$g,
          profile.par = c("fhg","k1"),
          fhg = v[1],
          k1 = v[2],
          method = "L-BFGS-B",
          lower = c(0.001,0.001),
          control=list(fnscale = -1))$val)

####.............profile contour bm1to2.coop...............#####
####fhg.k1####

grid.coop.fhg.k1 <- data.frame(
  t(expand.grid(unlist(
    bm1to2.coop_nls$profile.vals$fhg), 
    unlist(bm1to2.coop_nls$profile.vals$k1))))

bm1to2.coop_nls$lik.profile$fhg.k1 <- sapply(
  grid.coop.fhg.k1, function(v) 
    optim(bm1to2.coop_nls$nls.coef,
          ll.1to2.coop,
          obs = fi,
          g.curr = bm1to2.coop_nls$g,
          profile.par = c("fhg","k1"),
          fhg = v[1],
          k1 = v[2],
          method = "L-BFGS-B",
          lower = c(0.001,0.001),
          control=list(fnscale = -1))$val)
####.............Profile Plot bm1to2.deg.add...............####
####K1
profile_deg1 <- single.prof.plot(bm1to2.deg.add_nls, mle.deg.add, "k1")+
  mytheme_axes_box+
  ylab(TeX("$LogLik$",
           bold = TRUE))+
  xlab(TeX("$\\beta_{12}$ $\\[ µm^{-2} \\]$",
           bold = TRUE))+
  theme(legend.position = "none")

####fhg
profile_deg2 <- single.prof.plot(bm1to2.deg.add_nls, mle.deg.add, "fhg")+
  mytheme_axes_box+
  xlab(TeX("$F_{\\Delta HG}$ $\\[ µm^{-1} \\]$",
           bold = TRUE))+
  ylab(TeX("LogLik",
           bold = TRUE))+
  theme(legend.position = "none")

####.............Profile Plot bm1to2.coop.add............####
####k2
profile_coop.add <- single.prof.plot(bm1to2.coop.add_nls, mle.coop.add, "k2")+
  mytheme_axes_box+
  xlab(TeX("$K_2$ $\\[ µm^{-1} \\]$",
           bold = TRUE))+
  ylab(TeX("LogLik",
           bold = TRUE))+
  theme(legend.position = "none")

####.............Profile Plot bm1to2.coop...............####
####fhg
profile_coop1 <- single.prof.plot(bm1to2.coop_nls, mle.coop, "fhg")+
  mytheme_axes_box+
  xlab(TeX("$F_{\\Delta HG}$ $\\[ µm^{-1} \\]$",
           bold = TRUE))+
  ylab(TeX("LogLik",
           bold = TRUE))+
  theme(legend.position = "none")

####K1
profile_coop2 <- single.prof.plot(bm1to2.coop_nls, mle.coop, "k1")+
  mytheme_axes_box+
  xlab(TeX("$K_1$ $\\[ µm^{-1} \\]$",
           bold = TRUE))+
  ylab(TeX("LogLik",
           bold = TRUE))+
  theme(axis.text.x = element_text(angle = 25))+
  theme(legend.position = "none")

####.............Profile Contour Plots bm1to2.deg.add...............####
####fhg.k1
profile_deg3 <- contour.plot(bm1to2.deg.add_nls, 
                             mle.deg.add, 
                             theta = c("fhg", "k1"))+
  ylab(TeX("$F_{\\Delta HG}$ $\\[ µm^{-1} \\]$",
           bold = TRUE))+
  xlab(TeX("$\\beta_{12}$ $\\[ µm^{-2} \\]$",
           bold = TRUE))+
  mytheme_axes_box+
  theme(legend.position = "none")

####.............Profile Contour Plots bm1to2.coop......................####
####fhg.k1
profile_coop3 <- contour.plot(bm1to2.coop_nls, 
                              mle.coop, 
                              theta = c("fhg", "k1"))+
  mytheme_axes_box+
  ylab(TeX("$F_{\\Delta HG}$ $\\[ µm^{-1} \\]$",
           bold = TRUE))+
  xlab(TeX("$K_1$ $\\[ µm^{-1} \\]$",
           bold = TRUE))+
  theme(axis.text.x = element_text(angle = 25))+
  theme(legend.position = "none")

ggarrange(profile_deg2, profile_coop1, 
          profile_deg1, profile_coop2, 
          profile_deg3, profile_coop3,
          profile_coop.add,
          labels = c("A", "B", "C", "D", "E", "F", "G"),
          ncol = 2, nrow = 4)
**Maximum likelihood profiling of model parameter**. To assess the extent of confidence in the model-specific parameter, a parameter identifiability analysis was performed based on maximum likelihood profiling. During the profiling process a parameter under study is fixed to specified values (abscissa) around its best fit nls point estimate (dashed line) while the remaining parameter of the log-likelihood function were re-optimized (ordinate). The shape of the profile around its maximum was used to qualitatively judge parameter identifiability. **B, D, F:** Maximum likelihood profiles of $F_{\Delta HG}$ and $K_1$ from the bm1to2.coop model. $F_{\Delta HG}$ (B) is an identifiable parameter with a unique maximum of the profile at its best fit nls estimate (3.62 $\times$ 10<sup>-3</sup> µM<sup>-1</sup>) and a closed $CI_{0.95}$ (red line). $K_1$ (D) is a structurally non-identifiable parameter yielding a profile lacking a unique maximum and having an infinite $CI_{0.95}$. The joint likelihood profile (F) of both parameters exhibits parallel likelihood contours. Shades from black to white correspond to high and low values of the likelihood function, respectively. Likelihood contours are shown as black lines and best fit parameter values are indicated by a blue dot. **A, C, E:** Maximum likelihood profiles of $F_{\Delta HG}$ and $\beta_{12}$ from the bm1to2.deg.add model. Both parameter are identifiable resulting in a closed $CI_{0.95}$ contour in the two-dimensional profile likelihood landscape. **G** Maximum likelihood profile of $K_2$ from the bm1to2.coop.add model. $K_2$ is a practically non-identifiable parameter. Its profile exhibits a unique maximum but the $CI_{0.95}$ only has an upper bound.

Figure 4.2: Maximum likelihood profiling of model parameter. To assess the extent of confidence in the model-specific parameter, a parameter identifiability analysis was performed based on maximum likelihood profiling. During the profiling process a parameter under study is fixed to specified values (abscissa) around its best fit nls point estimate (dashed line) while the remaining parameter of the log-likelihood function were re-optimized (ordinate). The shape of the profile around its maximum was used to qualitatively judge parameter identifiability. B, D, F: Maximum likelihood profiles of \(F_{\Delta HG}\) and \(K_1\) from the bm1to2.coop model. \(F_{\Delta HG}\) (B) is an identifiable parameter with a unique maximum of the profile at its best fit nls estimate (3.62 \(\times\) 10-3 µM-1) and a closed \(CI_{0.95}\) (red line). \(K_1\) (D) is a structurally non-identifiable parameter yielding a profile lacking a unique maximum and having an infinite \(CI_{0.95}\). The joint likelihood profile (F) of both parameters exhibits parallel likelihood contours. Shades from black to white correspond to high and low values of the likelihood function, respectively. Likelihood contours are shown as black lines and best fit parameter values are indicated by a blue dot. A, C, E: Maximum likelihood profiles of \(F_{\Delta HG}\) and \(\beta_{12}\) from the bm1to2.deg.add model. Both parameter are identifiable resulting in a closed \(CI_{0.95}\) contour in the two-dimensional profile likelihood landscape. G Maximum likelihood profile of \(K_2\) from the bm1to2.coop.add model. \(K_2\) is a practically non-identifiable parameter. Its profile exhibits a unique maximum but the \(CI_{0.95}\) only has an upper bound.

Maximum likelihood profiling revealed that the cooperative binding models each contain a non-identifiable parameter, whereby first conclusions about the cause could be drawn from the type of non-identifiability. Practical non-identifiability as seen in the bm1to2.coop.add model is linked to the noise level of the data while structural non-identifiability as seen in the global bm1to2.coop model is linked to over-parameterization. This has to be taken into account for model evaluation steps. Following this, considered binding models are not on an equal footing with respect to the number of parameters so that the \(RMSE\) of the residuals alone is not sufficient to evaluate model performances since it does not account for model complexity. For this reason the so-called Akaike Information Criterion (\(AIC\)) was used in the first stage of binding model evaluation.

4.1.3 Kullback-Leibler Divergence Analysis

In information theory, the information lost when a model is used to approximate the true data generating process is described by the Kullback-Leibler (KL) distance. Hence, the KL divergence between a model and the data generating process can be used as a measure of performance for that particular model. The truly information theoretic KL divergence can be approximated by the \(AIC\) which can be used in the context of nls regression analysis. It thereby takes into account the goodness of fit as determined by the \(RMSE\) of the residuals as well as model complexity in terms of the number of parameter. Calculated \(AIC\) scores for the binding models were used as a figure of merit in order to evaluate model performance in describing the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) binding profile.

In this work the second order \(AIC\) corrected for small sample sizes (\(AIC_c\)) was used (3.2.4). As a first step of the model evaluation procedure the \(AIC_c\) scores for every binding model were calculated based on the respective residual sum of squares obtained from the nls fits (4.1.1) according to equation (3.14). These scores shown in table 4.2 and figure 4.3A represent approximated KL divergences on a negative ordinal scale between the binding models and the true but unknown data-generating process (“reality”, red square) which can be treated as a constant across all models. The lower the \(AIC_c\) score of a model (radial grid lines), the smaller the KL divergence and the better its performance in describing the data. In this way, calculated \(AIC_c\) scores point towards the bm1to2.deg.add model (blue diamond) to be the KL best model in the set since it received the lowest \(AIC_c\) score amounting to -52.4. Competing binding models were further ranked by calculation of \(AIC_c\) differences (\(\Delta_i\)) for each binding model with respect to the top-ranked bm1to2.deg.add model (equation (3.15)). In this way the unknown “reality” is canceled out and \(\Delta_i\) values now reflect increments of KL divergences between the bm1to2.deg.add model (\(\Delta_i\) = 0) and all other competing models (figure 4.3B). Those \(\Delta_i\) values were used as a basis to quantify the plausibility of each model as being the KL best model in the set instead of the actual top-ranked bm1to2.deg.add model. This was achieved by normalization of the \(\Delta_i\) values to be a set of Akaike weights \(\omega_i\) (evidence weights) adding to 1 as described by equation (3.16). Additionally, evidence ratios \(\frac{\omega_i}{\omega_j}\) (equation (3.17)) between the bm1to2.deg.add model and the remaining models in the set were computed. As a result of these further calculations, the top-ranked bm1to2.deg.add model received 70.8% of the total weight of evidence of being the actual KL best model considering all models in the set. According to a plausibility key introduced by Anderson and Burnham (2004) which is based on the \(\Delta_i\) value, the bm1to2.deg.add model was the only one found to have substantial evidence of being the KL best model since no other model received a \(\Delta_i\) value less than 2. The bm1to2.coop.add model is ranked the second best in the set having an \(AIC_c\) score of -49.0 corresponding to a \(\Delta_i\) value of 3.35 which accounts for a weight of evidence of 13.3%. This and the remaining models all having \(\Delta_i\) values greater than 2 but less than 7 would be considered as having considerable less evidence of being the KL best model in the set. Also considering evidence ratios, the bm1to2.deg.add model is 5.3 times more likely to be the KL best model instead of the bm1to2.coop.add model. The evidence ratio for the bm1to2.deg.add model versus the two next best models ranged from 10.6 to 11.9 including the bm1to2.deg and the global bm1to2.coop models, respectively. The model ranked lowest comprises the one-site binding model bm1to1 which got the highest \(AIC_c\) score (-46.2) and a weight of evidence of 3.19% of being the KL best model in the set. Hence, it is 22.2 times more likely that the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) binding profile was generated by the bm1to2.deg.add model than by the bm1to1 model.

Until that point analytic steps involving the evaluation of model performance were all conditional on the data set based on which the nls regression itself was carried out. In order to evaluate how generalizable those performances would be on new data, binding models were additionally subjected to cross-validation.

**Ranking of binding models according to their AIC<sub>c</sub> scores. A:** AIC<sub>c</sub> scores evaluating model performances in describing the ThiI<sub>Tm</sub>•$tRNA^{Phe}_{Bs}$ titration profile are highlighted as radial grid lines in the radar plot. They represent approximated KL distances on an ordinal scale with respect to the true data-generating process ("reality", red diamond). "Reality" represents the origin of the plot and is arbitrarily set to -60 for visualization purposes. Models were ranked based on the AIC<sub>c</sub> scores with the bm1to2.deg.add model (blue diamond) being the top-ranked model which got the smallest AIC<sub>c</sub> score (-52.4) and therefore best describes the binding curve. **B:** Binding models were further ranked based on AIC<sub>c</sub> differences ($\Delta_i$) with respect to the top-ranked bm1to2.deg.add model ($\Delta_i = 0$) so that "reality" treated as a constant across all models is canceled out. Calculated $\Delta_i$ values point towards the bm1to2.coop.add model as being the most competitive model ($\Delta_i = 3.3$) in the set with respect to the top-ranked bm1to2.deg.add model.

Figure 4.3: Ranking of binding models according to their AICc scores. A: AICc scores evaluating model performances in describing the ThiITm•\(tRNA^{Phe}_{Bs}\) titration profile are highlighted as radial grid lines in the radar plot. They represent approximated KL distances on an ordinal scale with respect to the true data-generating process (“reality”, red diamond). “Reality” represents the origin of the plot and is arbitrarily set to -60 for visualization purposes. Models were ranked based on the AICc scores with the bm1to2.deg.add model (blue diamond) being the top-ranked model which got the smallest AICc score (-52.4) and therefore best describes the binding curve. B: Binding models were further ranked based on AICc differences (\(\Delta_i\)) with respect to the top-ranked bm1to2.deg.add model (\(\Delta_i = 0\)) so that “reality” treated as a constant across all models is canceled out. Calculated \(\Delta_i\) values point towards the bm1to2.coop.add model as being the most competitive model (\(\Delta_i = 3.3\)) in the set with respect to the top-ranked bm1to2.deg.add model.

Table 4.2: Results of the KL divergence analysis. Performance of binding models in the nls regression of the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) binding isotherm was evaluated based on \(AIC_c\) scores according to which models are ranked in this table. \(RMSE\) values are those transferred from table 3.1. \(\Delta_i\) = AICc differences, \(\omega_i\) = evidence weights, CUSUM = cumulative sum of the evidence weights, k = number of parameter including an additional parameter for estimating the variance.
\(Model\) \(k\) \(RMSE\) \(AIC_c\) \(\Delta_i\) \(\omega_i\) \(CUSUM\)
bm1to2.deg.add 3 0.0205627 -52.35595 0.000000 70.847168 70.84717
bm1to2.coop.add 4 0.0204736 -49.01018 3.345772 13.298334 84.14550
bm1to2.deg 4 0.0216753 -47.64134 4.714610 6.707457 90.85296
bm1to2.coop 5 0.0178691 -47.40347 4.952481 5.955316 96.80828
bm1to1 3 0.0266239 -46.15603 6.199927 3.191725 100.00000

4.1.4 Quantifying Model Generelizability

Cross-validation (CV) offers the possibility to evaluate model performance on future data based on a single data set. In order to assess how generalizable considered binding models are a complete 4-fold cross validation was performed (3.2.5). This was accomplished by partitioning the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) titration data set consisting of 12 data points in total into test-folds containing three data points and into disjoint train-folds containing the complementary set of data points. This was repeated for a total of 220 test-folds comprising every possible combination of three data points. Following the partitioning, nls regression analysis was carried out for each binding model on the train-fold combinations and model-specific parameter point estimates obtained in this way were then kept constant and used to predict the data-triplet in the complementary test-fold. The predictive power was iteratively assessed by computation of a CV-score representing the \(RMSE\) between the predicted and the observed data. The distribution of these 220 individual CV-scores are shown for each binding model as grey dots in figure 4.4A. The mean value which hereafter is denoted as the mean fold prediction error is highlighted as a red square and the standard deviation as a black line. The mean fold prediction errors are also highlighted as blue bars in figure 4.4B and are compared to the \(RMSE\)s obtained from the respective nls regressions on the complete titration data set (red bars, cf. table 4.1) which hereafter are denoted as nls prediction errors.

The complete cross-validation revealed that the bm1to1 model is the least predictive model in the set. It received the highest absolute mean fold prediction error amounting to 3.44 \(\times\) 10-2. Also, it has the largest deviation from the nls prediction error (\(\Delta RMSE\)) amounting to 7.76 \(\times\) 10-3. The global bm1to2.coop model which obtained a mean fold prediction error of 1.98 \(\times\) 10-2 that is best balanced with the nls prediction error (\(\Delta RMSE\) = 1.95 \(\times\) 10-3) is the best predictive model in the set. In between are, respectively, the bm1to2.deg.add and bm1to2.deg models sharing a mean fold prediction error of 2.37 \(\times\) 10-2 as well as the bm1to2.coop.add model whose mean fold prediction error amounts to 2.25 \(\times\) 10-2. The fold prediction errors of these models are also well balanced with their respective nls prediction errors comparable in size to the bm1to2.coop model.

fi <- fi[,1:2]
bm1to1_nls <- k_fold_cross(data = fi,
                           model.list = bm1to1_nls,
                           model.fun = c("bm1to1"),
                           nFolds = 4)

bm1to2.deg_nls <- k_fold_cross(
  data = fi, 
  model.list = bm1to2.deg_nls, 
  model.fun = c("bm1to2.deg"),
  nFolds = 4)

bm1to2.deg.add_nls <- k_fold_cross(
  data = fi, 
  model.list = bm1to2.deg.add_nls, 
  model.fun = c("bm1to2.deg.add"),
  nFolds = 4)

bm1to2.coop.add_nls <- k_fold_cross(
  data = fi, 
  model.list = bm1to2.coop.add_nls, 
  model.fun = c("bm1to2.coop.add"),
  nFolds = 4)

bm1to2.coop_nls <- k_fold_cross(
  data = fi, 
  model.list = bm1to2.coop_nls, 
  model.fun = c("bm1to2.coop"),
  nFolds = 4)

cross_valid.ls <- list(bm1to1_nls,
                       bm1to2.deg.add_nls,
                       bm1to2.deg_nls,
                       bm1to2.coop.add_nls,
                       bm1to2.coop_nls)
names(cross_valid.ls) <- c("bm1to1", 
                           "bm1to2.deg.add", 
                           "bm1to2.deg", 
                           "bm1to2.coop.add",
                           "bm1to2.coop")

cross_error <- cross_valid.ls %>% 
  map("cross_valid") %>% 
  flatten()

model_sigma <- cross_valid.ls %>% 
  map(`[[`, 1) %>% 
  map_dbl(~ sigma(.x))

nRow <- 12
nFolds <- 4
fold_comb <- combn(nRow, nRow / nFolds)
sublist_app <- map(seq_len(ncol(fold_comb)), ~fold_comb[,.])
sublist_app_chr <- sublist_app %>% 
  map(`[`, c(1,2,3)) %>%
  map_chr(~paste0(.x, collapse = "-")) 

compCV.tbl <- tibble(model = factor(rep(names(cross_valid.ls), each = 220), 
                                   levels = names(cross_valid.ls)),
                    RMSE = unlist(cross_error),
                    fold = rep(c(1:220), times = 5),
                    fold_points = rep(sublist_app_chr, times = 5))

compCV_bar.tbl <- compCV.tbl %>% 
  group_by(model) %>% 
  summarise(cv_RMSE = mean(RMSE)) %>% 
  ungroup() %>% 
  mutate(nls_RMSE = model_sigma)

###############Plots#####################
posn.d <- position_dodge(0.9)
plotA <- ggplot(compCV.tbl, aes(y = RMSE, 
                                x = model))+
  geom_jitter(width = 0.3, color = "grey50", alpha = 0.6)+
  stat_summary(fun.y = mean, geom = "point", 
               position = posn.d,
               fill = "red", size = 3.5, shape = 23)+
  stat_summary(fun.data = mean_sdl, 
               fun.args = list(mult = 1), 
               geom = "errorbar", width = 0, 
               position = posn.d, size = 1.1)+
  mytheme_axes_grids+
  theme(axis.line.x = element_blank(),
        axis.text.x = element_text(angle = 15))

plotB <- compCV_bar.tbl %>% 
  pivot_longer(cols = - model, 
               names_to = "Prediction",
               values_to = "RMSE") %>% 
  mutate(Prediction = factor(Prediction,
                             levels = c("nls_RMSE", "cv_RMSE"))) %>%
  
  ggplot(aes(x = model, y = RMSE, fill = Prediction))+
  geom_bar(stat = "identity",position = "dodge")+
  scale_fill_brewer(palette = "Set1")+
  mytheme_axes_grids+
  theme(legend.position = c(0.8, 0.9),
        axis.line.x = element_blank(),
        axis.text.x = element_text(angle = 15))

ggarrange(plotA, plotB, 
          labels = c("A", "B"),
          ncol = 2, nrow = 1)
**Predictive power of binding models determined by cross-validation. A:** A complete 4-fold cross-validation was conducted in order to assess binding model generalizability on new data. The data set was iteratively split into test-folds containing three data points and into train-folds containing the complementary set of data points. For each binding model least squares regression was performed on the train-folds and obtained best fit parameter estimates were then used to predict the data points in the test-folds. This was iterated over each possible test-fold combination comprising three data points (220 times). For each combination the $RMSE$ was determined (grey dots) as well as the mean $RMSE$ (red square) and the standard deviation (black line). **B:** Mean $RMSE$ from A (blue) as well as the $RMSE$ obtained from the nls regression (red) are highlighted as bars. The bm1to1 model is the least predictive model while the bm1to2.coop model is the most predictive model in the set.

Figure 4.4: Predictive power of binding models determined by cross-validation. A: A complete 4-fold cross-validation was conducted in order to assess binding model generalizability on new data. The data set was iteratively split into test-folds containing three data points and into train-folds containing the complementary set of data points. For each binding model least squares regression was performed on the train-folds and obtained best fit parameter estimates were then used to predict the data points in the test-folds. This was iterated over each possible test-fold combination comprising three data points (220 times). For each combination the \(RMSE\) was determined (grey dots) as well as the mean \(RMSE\) (red square) and the standard deviation (black line). B: Mean \(RMSE\) from A (blue) as well as the \(RMSE\) obtained from the nls regression (red) are highlighted as bars. The bm1to1 model is the least predictive model while the bm1to2.coop model is the most predictive model in the set.

The effect of the binding models on the respective CV-score distributions was further investigated by means of analysis of variance (ANOVA). However, as indicated by figure 4.4A CV-scores are not normally distributed. This is confirmed by figure 4.5A where the CV-score distributions are summarized by boxplots. These plots highlight outlier with respect to every binding model. Further analysis of the outlier revealed that data points one and two of the tiration data set are especially prevalent amongst the data point triplets of the test-folds based on which CV-scores were computed. Due to the fact that each binding model is affected by outlier it is to be expected that the errors (residuals) of the ANOVA model are also not going to be normally distributed. However, the generic ANOVA approach crucially depends on the assumption of normally distributed residuals. For this reason, a robust ANOVA was additionally conducted by using CV-score distributions whose tails were trimmed by 20% (figure 4.5B).

**CV-score distributions are affected by test- fold combinations. A:** The raw CV-score distributions shown in figure 3.4A are summarized by box plots across levels of binding model. Outlier are shown as dots and labeled by data point triplets comprising the respective test fold. The presence of outlier leads to a bimodal distribution of CV-scores violating the normality assumption. Data points one and two are especially prevalent amongst outlier test folds. Note that some data point triplets were omitted from the plot in order to prevent over-plotting. **B:** Shown are box plots of the CV-score distributions across levels of binding model after having trimmed 20% of the raw data from each side. The trimming has resolved the normality assumption about CV-scores. Generally, outlier are defined as the set of CV-scores according to the following formula: $X_{out} = \{X|X < Q_1 - 1.5 \times IQR \lor X > Q_3 + 1.5 \times IQR \}$, where $Q_1$ and $Q_3$ are, respectively, the first and third quartiles and $IQR$ represents the interquartile range.

Figure 4.5: CV-score distributions are affected by test- fold combinations. A: The raw CV-score distributions shown in figure 3.4A are summarized by box plots across levels of binding model. Outlier are shown as dots and labeled by data point triplets comprising the respective test fold. The presence of outlier leads to a bimodal distribution of CV-scores violating the normality assumption. Data points one and two are especially prevalent amongst outlier test folds. Note that some data point triplets were omitted from the plot in order to prevent over-plotting. B: Shown are box plots of the CV-score distributions across levels of binding model after having trimmed 20% of the raw data from each side. The trimming has resolved the normality assumption about CV-scores. Generally, outlier are defined as the set of CV-scores according to the following formula: \(X_{out} = \{X|X < Q_1 - 1.5 \times IQR \lor X > Q_3 + 1.5 \times IQR \}\), where \(Q_1\) and \(Q_3\) are, respectively, the first and third quartiles and \(IQR\) represents the interquartile range.

4.1.5 CV-ANOVA

In order to test whether the central tendency of the CV-score distributions (4.1.4) is affected by the type of binding model, as suggested by figure 4.4A, an analysis of variance (ANOVA) was conducted (3.2.6). ANOVA is generally used to study the effect that one or more factors might have on the variance in a response variable based on the respective factor means. The type of ANOVA used in this study was specifically designed to investigate the binding model as a factor containing five level in a one-way independent factorial design with the use of a complete set of orthogonal linear contrasts (eq. (3.21)). This design was formalised as a genral linear regression model (\(\overrightarrow{Y} = \overrightarrow{\beta}X + \overrightarrow{\epsilon}\), see eq. (3.19)) where each regression coefficient associated with a contrast coded predictor variable specifies a particular difference of group means on the basis of those linear contrasts. During the ANOVA each contrast is tested against the null hypothesis (\(H_0\), eq. (3.22)) via the appropriate \(F\)-stastic (eq. (3.24)). In a first approach the distributional requirements of the residuals in order to use the theorectic \(F\)-distribution under \(H_0\) were assumed to be true. The summary statistics of the ANOVA are highlighted in table 4.3.

compCV.tbl <- compCV.tbl%>% 
  mutate(fold_class = ifelse(grepl("1-2", fold_points),
                                  "outlier",
                                  "normal"),
         fold_class = factor(fold_class, levels = c("normal", 
                                                    "outlier")))

custom_cont <- matrix(c(-4/5,1/5,1/5,1/5,1/5,
                        0,1/4,1/4,1/4,-3/4,
                        0,-1/3,-1/3,2/3,0,
                        0,-1/2,1/2,0,0), 
                      nrow = 5)
colnames(custom_cont) <- c("k1", "k2", "k3", "k4")
rownames(custom_cont) <- c("bm1to1", "bm1to2.deg.add", 
                      "bm1to2.deg", "bm1to2.coop.add",
                      "bm1to2.coop")
####################LinearReg and ANOVA#####################
contrasts(compCV.tbl$model) <- custom_cont
contrasts(compCV.tbl$fold_class) <- c(1,-1)
compModel.lm <- lm(RMSE ~ 1, data = compCV.tbl)
augModel1.lm <- lm(RMSE ~ model, data = compCV.tbl)
augModel2.lm <- lm(RMSE ~ fold_class + model, data = compCV.tbl)

augModel1.aov <- aov(augModel1.lm)
augModel2.aov <- aov(augModel2.lm)

cont_summary1 <- summary(augModel1.aov,
        split = list(
          model = list(
            k1 = 1, k2 = 2, k3 = 3, k4 = 4)
        ))

####Robust linear Reg and ANOVA #######################
#### Linear Regression
trim <- 0.2
n <- 220
cutoff <- floor(n * trim) 

compCV.trim.tbl<- compCV.tbl %>% 
  group_by(model) %>% 
  arrange(model,RMSE) %>% 
  mutate(index = c(1:220)) %>% 
  filter(index > cutoff & index <= 220 - cutoff) %>%
  dplyr::select(-index) %>% 
  ungroup()

augModel1.trim.lm <- lm(RMSE ~ model, data = compCV.trim.tbl)

#### ANOVA
# bring tidy tibble into list mode
compCV.ls <- compCV.tbl %>% 
  dplyr::select(c(model, RMSE)) %>% 
  mutate(index = rep(c(1:220), times = 5)) %>% 
  pivot_wider(names_from = model, values_from = RMSE) %>% 
  dplyr::select(-index) %>% 
  as.list()

robustComp <- lincon(compCV.ls, custom_cont)
options(knitr.kable.NA = '')

pre <- function(sum.aov){
  nRow <- nrow(sum.aov[[1]][2])
  ssr <- sum.aov[[1]][2][-nRow,]
  sseA <- sum.aov[[1]][2][nRow,]
  pre <- ssr / (ssr + sseA)
  return(pre)
}

#####defining ANOVA summary table####
linconts <- diag(apply(custom_cont, 2, 
                  function(x) augModel1.lm$coefficients[-1] * sum(x^2)))
aov.compModel <- tidy(anova(compModel.lm))
aov.df <- data.frame(Source = c("AugModel", "psi1", "psi2", "psi3", "psi4",
                                "Residuals", "Total"),
                     psihat = c(NA, linconts, NA, NA),
                     df = unlist(c(cont_summary1[[1]][1], 
                                   aov.compModel[1,"df"])),
                     SS = unlist(c(cont_summary1[[1]][2], 
                                       aov.compModel[1,"sumsq"])),
                     MS = unlist(c(cont_summary1[[1]][3], 
                                   aov.compModel[1,"meansq"])),
                     Fstat = unlist(c(cont_summary1[[1]][4],NA)),
                     p = unlist(c(cont_summary1[[1]][5],NA)),
                     PRE = c(pre(cont_summary1), NA, NA))
aov.df <- aov.df %>% 
  mutate(across(where(is.numeric), ~ round(., digits = 4)))
rownames(aov.df) <- c("$Y_{ij} = \\beta_0 + \\beta_1X_1 + \\beta_2X_2 + \\beta_3X_3 + \\beta_4X_4 + \\epsilon_{ij}$",
                      "$Y_{ij} = \\beta_0 + 0 + \\beta_2X_2 + \\beta_3X_3 + \\beta_4X_4 + \\epsilon_{ij}$",
                      "$Y_{ij} = \\beta_0 + \\beta_1X_1 + 0 + \\beta_3X_3 + \\beta_4X_4 + \\epsilon_{ij}$",
                      "$Y_{ij} = \\beta_0 + \\beta_1X_1 + \\beta_2X_2 + 0 + \\beta_4X_4 + \\epsilon_{ij}$",
                      "$Y_{ij} = \\beta_0 + \\beta_1X_1 + \\beta_2X_2 + \\beta_3X_3 + 0 + \\epsilon_{ij}$",
                      "---",
                      "$Y_{ij} = \\beta_0  + \\epsilon_{ij}$")

aov.df %>% 
  kable(caption = "**ANOVA summary table.** Shown is a composite table including summary statistics of the ANOVA but also the estimated contrasts ($\\hat{\\Psi}$) and effect sizes ($\\hat{\\eta}^2$). The $SS$ column represents sums of squares. When the source is the augmented model (\"AugModel\") the $SS$ value represents a reduction in sums of squares ($SSR$) on 4 degrees of freedom ($df$). This is also true for the partitioning of the $SSR$ of the augmented model into $SSR$ components of individual contrasts (\"psi1\"-\"psi4\") based on 1 $df$. When the source is \"Residuals\" the $SS$ value refers to the sum of squared error ($SSE$) of the augmented model and when the source is \" Total\" the $SS$ value refers to the $SSE$ of the grand mean model. In the $MS$ and $F_{stat}$ columns, respectively, corresponding mean squared errors and observed $F$-statistics are listed. The $p$ value represents the conditional probability of observing an $F$-statistic as extreme as computed for this instance or even more extreme. The $\\hat{\\eta}^2$ value represents the proportional reduction in error computed as $\\hat{\\eta}^2 = \\frac{SSR_i}{SSR_i + SSE_A}$, where $SSE_A$ represents the sum of suqred error of the augmented model.",      
        col.names = c("$Source$",
                  "$\\hat{\\Psi}$",
                    "$df$",
                    "$SS$",
                    "$MS$",
                    "$F_{stat}$",
                    "$p(\\geq F_{stat})$",
                    "$\\hat{\\eta}^2$")) %>% 
  kable_classic() %>% 
  row_spec(0, background = "grey")
Table 4.3: ANOVA summary table. Shown is a composite table including summary statistics of the ANOVA but also the estimated contrasts (\(\hat{\Psi}\)) and effect sizes (\(\hat{\eta}^2\)). The \(SS\) column represents sums of squares. When the source is the augmented model (“AugModel”) the \(SS\) value represents a reduction in sums of squares (\(SSR\)) on 4 degrees of freedom (\(df\)). This is also true for the partitioning of the \(SSR\) of the augmented model into \(SSR\) components of individual contrasts (“psi1”-“psi4”) based on 1 \(df\). When the source is “Residuals” the \(SS\) value refers to the sum of squared error (\(SSE\)) of the augmented model and when the source is " Total" the \(SS\) value refers to the \(SSE\) of the grand mean model. In the \(MS\) and \(F_{stat}\) columns, respectively, corresponding mean squared errors and observed \(F\)-statistics are listed. The \(p\) value represents the conditional probability of observing an \(F\)-statistic as extreme as computed for this instance or even more extreme. The \(\hat{\eta}^2\) value represents the proportional reduction in error computed as \(\hat{\eta}^2 = \frac{SSR_i}{SSR_i + SSE_A}\), where \(SSE_A\) represents the sum of suqred error of the augmented model.
\(Source\) \(\hat{\Psi}\) \(df\) \(SS\) \(MS\) \(F_{stat}\) \(p(\geq F_{stat})\) \(\hat{\eta}^2\)
\(Y_{ij} = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3 + \beta_4X_4 + \epsilon_{ij}\) AugModel 4 0.0274 0.0068 30.9648 0.0000 0.1016
\(Y_{ij} = \beta_0 + 0 + \beta_2X_2 + \beta_3X_3 + \beta_4X_4 + \epsilon_{ij}\) psi1 -0.0096 1 0.0252 0.0252 113.9305 0.0000 0.0942
\(Y_{ij} = \beta_0 + \beta_1X_1 + 0 + \beta_3X_3 + \beta_4X_4 + \epsilon_{ij}\) psi2 0.0026 1 0.0020 0.0020 9.0324 0.0027 0.0082
\(Y_{ij} = \beta_0 + \beta_1X_1 + \beta_2X_2 + 0 + \beta_4X_4 + \epsilon_{ij}\) psi3 -0.0008 1 0.0002 0.0002 0.8964 0.3439 0.0008
\(Y_{ij} = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3 + 0 + \epsilon_{ij}\) psi4 0.0000 1 0.0000 0.0000 0.0000 0.9996 0.0000
— Residuals 1095 0.2420 0.0002
\(Y_{ij} = \beta_0 + \epsilon_{ij}\) Total 1099 0.2694 0.0002

A regression analysis of the augmented (complete) model revealed that the type of binding model has an overall statistical effect on the CV-score distribution. CV-score predictions conditional on the binding model factor results in a sum of squared error (\(SSE\)) of 2.42 \(\times\) 10-1 (“Residuals” in Source column in tab. 4.3) as opposed to the simple grand mean model the \(SSE\) of which amounts to 2.69 \(\times\) 10-1 (“Total” in tab. 4.3). The resulting reduction in the sums of squares (\(SSR\)) of 2.74 \(\times\) 10-2 (“AugModel” in Source column) leads to a statistically discernible \(F\)-statistic (\(F_{4,1095}\)) of \(\sim\) 31. Conditional on a true \(H_0\), the probability \(p\) of obtaining the observed \(F\)-statistic or hypothetically more extreme ones is less than 2 \(\times\) 10-16. The augmented model was estimated to account for a proportional reduction in error (explained variance) of 10.2% as captured by the \(\hat{\eta}^2\) value. In order to address more focused hypotheses the \(SSR\) of the augmented model was further partitioned into \(SSR\) components (“psi1”-“psi4” in tab. 4.3) attributable to group mean differences as specified by the set of orthogonal linear contrasts.

The first hypothesis (\(H_{01}\) in eq. (3.22)) concerns the first linear contrast \({\Psi}_1\) encoding the difference between the mean CV-score of the bm1to1 model and the pooled mean CV-score of all tow-site binding models. \(H_{01}\) is tested by a focused 1-\(df\) model comparison involving the compact model \(C1\) (see eq. (3.23)) and the augmented model the former of which is obtained by setting \(\beta_{X_1} = 0\). This comparison revealed that including the contrast coded predictor in the regression model results in a \(SSR\) of 2.52 \(\times\) 10-2 which is statistically supported (\(F_{1,1095}\) = 113.9, \(p\) < 2 \(\times\) 10-16). The underlying contrast \(\Psi_1\) was estimated as -9.6 \(\times\) 10-3 (\(\Psi\) column in tab. 4.3) and computed to be situated in a compatible range of values from [-1.42 \(\times\) 10-2, to -9.76 \(\times\) 10-3] at the 95% confidence level (fig. 4.6 A). From the total proportional reduction in error due to the augmanted model 9.42% can be attributed to this contrast.

The second hypothesis \(H_{02}\) tests the difference between the observed mean CV-sore of the global bm1to2.coop model and the average mean CV-score of the remaining two-site binding models via comparison of compact model \(C_2\) (\(\beta_{X_2} = 0\)) with the augmented model. This difference is based on contrast \(\hat{\Psi}_2\) which was estimated as 2.61 \(\times\) 10-3 [\(CI_{0.95}\): 1.21 \(\times\) 10-3 to 5.75 \(\times\) 10-3]. Including \(\hat{\Psi}_2\) as a contrast coded predictor variable in the augmented model results in a \(SSR\) of 2.0 \(\times\) 10-3 that is statistically supported (\(F_{1,1095}\) = 9.03, \(p\) = 0.27%) and which results in a proportional reduction in error of 0.82%.

Hypothesis \(H_{03}\) which is represented by compact model \(C_3\) (\(\beta_{X_3} = 0\)) tests for the differnce in mean CV-scores between the bm1to2.coop.add model and the degenerative models bm1to2.deg and bm1to2.deg.add, respectively. The corresponding contrast \(\Psi_3\) was estimated to be -7.75 \(\times\) 10-4 [\(CI_{95}\): -3.57 \(\times\) 10-3 to 1.25 \(\times\) 10-3]. An inclusion of \(\Psi_3\) as predictor variable in the augmented model leads to a statistically non-discernible \(SSR\) of 2.0 \(\times\) 10-4 (\(F_{1,1095}\) = 8.96 \(\times\) 10-1, \(p\) = 34.4%).

The last hypothesis \(H_{04}\) covered by compact model \(C_4\) (\(\beta_{X_4} = 0\)) tests for the difference in mean CV-scores between the bm1to2.deg and bm1to2.deg.add models, respectively. Here, the contrast was estimated to be basically zero (\(\Psi_4\) = 3.76 \(\times\) 10-7 [\(CI_{95}\): -2.78 \(\times\) 10-3 to 2.78 \(\times\) 10-1 3]). The \(SSR\) amounts to zero when \(\Psi_4\) is included as predictor variable in the augmented model which is statistically not supported (\(F_{1, 1095}\) = 0, \(p\) ~ 100%).

ci.df <- data.frame(lower = contrast.ci(augModel1.lm)[,1],
                    upper = contrast.ci(augModel1.lm)[,2],
                    contrast = c("1","2","3","4"))

ci.A <- ci.df %>%
  pivot_longer(-contrast,names_to = "bound", values_to = "margin") %>% 
  mutate(coefficient = rep(linconts, each = 2)) %>% 
  ggplot(aes(x = margin, y = contrast))+
  geom_point(size = 2)+
  geom_line(size = 1.0)+
  geom_point(aes(x = coefficient), color = "blue", size = 3, shape = 18)+
  geom_vline(xintercept = 0, lty = 2, size = 1.0,
             color = "red")+
  theme_bw()+
  theme(axis.ticks.length=unit(.07, "cm"))+
  theme(axis.ticks = element_line(colour = "black", size = 1))+
  theme(axis.text.y = element_text(color="black",size = 12))+
  theme(axis.text.x = element_text(color="black",size = 12))+
  theme(axis.title.y = element_text(size = 14,face = "bold"))+
  theme(axis.title.x = element_text(size = 14,face = "bold"))+
  theme(panel.grid.minor.x = element_line(color = "grey80"),
        panel.grid.major.x = element_line(color = "grey80"),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.border = element_rect(size = 1.2, color = "black"))+
  xlab(TeX("$\\Delta$CV-Score", bold = TRUE))

#### Residual plot#####

resid.A <- 
  augModel1.lm %>% 
  augment() %>% 
  ggplot(aes(sample = .resid))+
  stat_qq(aes(color = model), alpha = 0.7, size = 1.3)+
  stat_qq_line(color = "black", size = 1.3)+
  scale_color_uchicago()+
  mytheme_axes_box+
  theme(legend.title = element_blank(),
        legend.text = element_text(size = 10,
                                   face = "bold"),
        legend.position = c(0.3,0.72))

ggarrange(resid.A, ci.A,
          labels = c("A","B"),
          ncol = 2, nrow = 1)
**A: ANOVA residual QQ-plots.** The residuals of the augmented model used to fit the CV-scores conditional on the type of binding model were sorted in ascending order (ordinate) and plotted against the quantiles from a standard normal distribution (abscissa). **B: Interval estimates at the 95% confidence level.** For each contrast (1-4) a 95% confidence interval was computed around the respective least squares point estimate (blue diamonds) based on a common variance term ($MSE_A$, see eq. 2.26) amongst each level of the binding model factor.

Figure 4.6: A: ANOVA residual QQ-plots. The residuals of the augmented model used to fit the CV-scores conditional on the type of binding model were sorted in ascending order (ordinate) and plotted against the quantiles from a standard normal distribution (abscissa). B: Interval estimates at the 95% confidence level. For each contrast (1-4) a 95% confidence interval was computed around the respective least squares point estimate (blue diamonds) based on a common variance term (\(MSE_A\), see eq. 2.26) amongst each level of the binding model factor.

The ANOVA conducted in this way is inter alia dependent on certain distributional assumptions about the residuals of the augmented model such as normality and constant variance (\(\epsilon_{ij}\overset{iid}{\sim} N(0,\sigma^2)\)). The raw CV-score distributions across the binding models as summarized by the box plots in figure 4.5A already highlight a clear violation of the normality assumption. This is due to the presence of outlier causing a bimodal distribution of CV-scores. That departure from normality is confirmed by the QQ-plot of the residuals of the ANOVA model (fig. 4.6) especially for the upper tail of the theoretic standard normal distribution (\(N(0,1)\)). For this reason a robust version of ANOVA was also applied in order to validate the current findings.

4.1.6 Robust ANOVA

The ANOVA approach as done in 4.1.5 is based on specific assumptions about the non-deterministic part of the general linear model such as a homoscedastic and normal distribution of its residuals. The ANOVA model can be made robust against deviations from those assumptions in various ways. In this work, the ANOVA model was specifically made robust against the assumption of normally distributed residuals by using trimmed means as outlined in 3.2.7.

The basic strategy was to replace, respectively, the mean and variance of the CV-score distributions of the binding models with 20% trimmed means (eq. (3.28)) and the corresponding Winsorized variances (eq. (3.30)). This replacement necessitated adjustments for the computation of the mean squared error (eq. (3.32)) and the degrees of freedom (eq. (3.33)) which were in part due to the use of effective sample sizes (\(h\)) after trimming. The \(F\)-statistic was replaced by an appropriate \(T\)-statistic(eq. (3.34)). Table 4.4 summarizes the robust ANOVA.

Table 4.4: Robust ANOVA summary table. In the \(\Psi\) column estimated contrast point estimates are presented. \(df\)-values represent estimated degrees of freedom based on effective sample sizes after trimming and taking into account a trimmed t-distribution. Respective mean squared errors (\(s_{yu}^2\)) were computed according to Yuen. The \(p\)-value is adjusted for multiple testing by referring to the C-variate studentized maximum modulus distribution.
\(Source\) \(\hat{\Psi}\) \(df\) \(s_{yu}^2\) \(T_{stat}\) \(p_{adj}(\geq T_{stat})\)
1 -0.0079780 183.7343 0.0004666 -17.0964097 0.0000000
2 0.0024198 165.8233 0.0005072 4.7707657 0.0000160
3 0.0004436 269.1250 0.0003191 1.3899299 0.5141157
4 0.0000003 262.0000 0.0002819 0.0010785 1.0000000

During the robust ANOVA approach, the same hypotheses as in 4.1.5 were tested which were then involving linear contrasts based on trimmed CV-score distribution means. For instance, the first hypothesis \(H_{01}\) tests the difference between the trimmed mean CV-score of the bm1to1 model and the average trimmed mean CV-score of all two-site binding models. This difference is based on linear contrast \(\Psi_1\) which was estimated as -7.98 \(\times\) 10-3 and which leads to a statistically discernible \(T\)-statistic of -17.1 as indicated by the corresponding adjusted \(p\) value (\(p_{adj}\)). The \(CI_{95}\) (eq. (3.35)) of \(\Psi_1\) and those for all other contrasts are shown in figure 4.7B and are in support of those determined in the non-robust ANOVA approach. Due to the trimming, the observed distribution of the residuals of the ANOVA model is more consistent with that of a theoretic standard normal distribution as highlighted in figure 4.7A.

ci.rob.df <- data.frame(lower = robustComp$psihat[,3],
                        upper = robustComp$psihat[,4],
                        contrast = c("1","2","3","4"))

ci.B <- ci.rob.df %>% 
  pivot_longer(-contrast,names_to = "bound", values_to = "margin") %>% 
  mutate(coefficient = rep(robustComp$psihat[,2], each = 2)) %>% 
  ggplot(aes(x = margin, y = contrast))+
  geom_point(size = 2)+
  geom_line(size = 1.0)+
  geom_point(aes(x = coefficient), color = "blue", size = 3, shape = 18)+
  geom_vline(xintercept = 0, lty = 2, size = 1.0,
             color = "red")+
  theme_bw()+
  theme(axis.ticks.length=unit(.07, "cm"))+
  theme(axis.ticks = element_line(colour = "black", size = 1))+
  theme(axis.text.y = element_text(color="black",size = 12))+
  theme(axis.text.x = element_text(color="black",size = 12))+
  theme(axis.title.y = element_text(size = 14,face = "bold"))+
  theme(axis.title.x = element_text(size = 14,face = "bold"))+
  theme(panel.grid.minor.x = element_line(color = "grey80"),
        panel.grid.major.x = element_line(color = "grey80"),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.border = element_rect(size = 1.2, color = "black"))+
  xlab(TeX("$\\Delta$CV-Score", bold = TRUE))

####residual plot####
resid.B <- augModel1.trim.lm %>% 
  augment() %>% 
  ggplot(aes(sample = .resid))+
  stat_qq(aes(color = model), alpha = 0.7, size = 1.3)+
  stat_qq_line(color = "black", size = 1.3)+
  scale_color_uchicago()+
  mytheme_axes_box+
  theme(legend.title = element_blank(),
        legend.text = element_text(size = 10,
                                   face = "bold"),
        legend.position = "none")

ggarrange(resid.B, ci.B,
          labels = c("A","B"),
          ncol = 2, nrow = 1)
**A: Robust ANOVA residual QQ-Plots.** The residuals of the augmented model used to fit the 20% trimmed CV-scores conditional on the type of binding model were sorted in ascending order (ordinate) and plotted against the quantiles from a standard normal distribution. **B: Interval estimates at the 95% confidence level.** For each contrast (1-4) a 95% confidence interval was determined based on Yuen's mean squared error ($s_{yu}^2$).

Figure 4.7: A: Robust ANOVA residual QQ-Plots. The residuals of the augmented model used to fit the 20% trimmed CV-scores conditional on the type of binding model were sorted in ascending order (ordinate) and plotted against the quantiles from a standard normal distribution. B: Interval estimates at the 95% confidence level. For each contrast (1-4) a 95% confidence interval was determined based on Yuen’s mean squared error (\(s_{yu}^2\)).

5 Discussion

In this excerpt special focus was on the chimeric ThiITm•\(\mathrm{tRNA^{Phe}_{Bs}}\) complex which was used as an example in order to dissect the true stoichiometry that underlies complexes that are formed between ThiITm and full-length tRNA substrates. For this reason tryptophan fluorescence spectroscopic binding studies were applied. An analytic work flow was established including ordinary nls regression, maximum likelihood profiling, model evaluation procedures such as KL divergence analysis and cross validation as well as generic and robust ANOVA approaches. On the basis of this analytic work flow it is possible to objectively infer this information from ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) binding profiles by statistical comparison of variably parameterized binding models all related to underlying 1:1 or 1:2 binding equilibria.

Applied model evaluation procedures together provide preliminary evidence that the underlying stoichiometry of ThiITm and its full-length tRNA substrates as determined in place of the ThiITm•\(\mathrm{tRNA^{Phe}_{Bs}}\) complex in the present work is indeed a 1:2 stoichiometry. These findings are discussed in more detail in the following sections.

5.1 Fluorescence Spectroscopic Binding Studies

ThiITm contains three tryptophan residues per monomer and is therefore characterized by a detectable intrinsic tryptophan fluorescence. Structural analysis of the ThiITm•TPHE39A complex (PDB 4kR7) revealed that two of these residues, namely Trp-192 and Trp-301, are situated in non-flexible parts buried inside the PPase domain. The remaining Trp-44, however, is located in the tRNA binding NFLD domain and is directly involved in interaction with TPHE39A by formation of a hydrogen bond with U-29 (figure 5.1). From this one can conclude that binding of a tRNA substrate generally leads to a change in the micro-environment of this residue and ultimately to a decrease of its fluorescence signal (Akbar, Sreeramulu, and Sharma 2016). This circumstance was taken advantage of by establishing a tryptophan fluorescence quenching assay in order to quantitatively characterize the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) interaction (4.1.1).

Remarkably, the recorded binding profile (cf. figure 4.1, 4.1.1) shows an increasing measurement uncertainty in relation to an increasing dilution factor of titrated \(\mathrm{tRNA^{Phe}_{Bs}}\) . This could most probably be attributed to the inherent dependency of the tryptophan signal on its environment. At high concentrations of \(\mathrm{tRNA^{Phe}_{Bs}}\) ThiITm is conformationally stabilized and so is the tryptophan fluorescence signal. At low concentrations of \(\mathrm{tRNA^{Phe}_{Bs}}\), however, conformational fluctuations of ThiITm become more pronounced, thereby destabilizing the tryptophan signal which might cause a higher measurement uncertainty. Therefore, the signal/noise ratio could be improved by increasing the conformational homogeneity of ThiITm. In general, cofactors are known to stabilize proteins (Mezzasalma et al. 2007). So far, the apo form of ThiITm has been used in the quenching assays so that further stabilization of ThiITm might be achieved by also including its cofactor ATP in the assay buffer. Additionally, the Thermofluor set up might be extended to include a large-scale additive screen in order to improve the buffer composition as described by Reinhard et al. (2013). Still, the ThiITm•\(\mathrm{tRNA^{Phe}_{Bs}}\) binding profile measured in triplicates already was of adequate quality for most binding model specific parameter sets as was revealed by subsequent parameter identifiability analysis.

Unveiling the true underlying stoichiometry of the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) complex is inherently linked to finding the most appropriate binding model that could describe the titration data and from which reliable inferences about the parameter can be made that describe the binding process. For this reason, a multi-model approach was pursued by comparatively analyzing a set of model functions that are related to underlying 1:1 and 1:2 binding equilibria. Generally, model functions were derived from the point of view that the quenching of the tryptophan fluorescence signal is of an entirely static nature due to a specific ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) interaction rather than due to unspecific collisions. For this reason, parameter describing dynamic quenching such as the Stern-Vollmer constant were not included in the model functions. Furthermore, focus was on those models that do neither approximate free ligand concentration nor assume simultaneous ligand binding. That is, every model function was derived based on the appropriate binding equilibrium by taking into account the law of mass conservation. That is why for the nls regression analysis of Adair-Klotz-type binding models computations of free ligand concentrations based on a cubic polynomial had to be integrated in the regression routine. Best fit regression curves obtained for each binding model revealed that all models are capturing the progress of the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) binding profile so that no differences in model performance could be detected based on visual inspection. Additionally computed \(RMSE\) of the residuals as a first goodness-of-fit indicator show the expected result that the binding curve is best fitted by the global bm1to2.coop model (cf. table 4.1, 4.1.1). This was the most complex model in the set and the number of parameter included in a model is usually negatively correlated with the \(RMSE\) of the residuals (Anderson and Burnham 2004). At the parameter level, however, the bm1to2.coop model yielded bio-physically impossible parameter estimates. The macroscopic binding constant K1, for instance, was estimated to be in the kilo molar range (1.38 \(\times\) 103 M-1) so that this model ultimately would describe tRNA\(\mathrm{^{Phe}_{Bs}}\) binding to ThiITm as extemely positively cooperative (11.8 \(\times\) 106). This points to the direction that the global model is over-parameterized which could not be determined by taking into account the \(RMSE\) only. In contrast, when model complexity is increased such as from the bm1to2.deg.add to the bm1to2.deg model which is accompanied by a concomittant increase of the \(RMSE\) from 2.06 \(\times\) 10-2 to 2.17 \(\times\) 10-2, then the \(RMSE\) is an already sufficient model performance metric. It indicates that there is not enough information content in the titration data for two proportionality factors as included in the bm1to2.deg model. This finding is supported by the similarity of the parameter values of these two models.

The nls regression analysis gave a first impression concerning model performances. However, based on the \(RMSE\) of the residuals it was not possible to unequivocally determine a best data describing binding model since considered models are not on an equal footing with respect to the number of parameter. Therefore, model evaluation was done by applying the \(AIC_c\) which is capable of trading off goodness-of-fit based on the \(RMSE\) and model complexity. Yet, the outcome of this analysis is dependent on the quality of the model set. For this reason a parameter identifiability analysis was carried out prior to model evaluation. Given the noise level of the data and that the global binding model yielded unusual parameter estimates this analysis should reveal whether in the context of the applied measurement type and set up especially the more complex 1:2 binding models were adequate at all.

**Hydrogen bonding interaction between NFLD Trp-44 and U-29.** Shown is a cartoon representation of the ThiI<sub>Tm</sub>•TPHE39A complex (PDB 4KR7) in a close-up view highlighting Trp-44 of the NFLD domain (red) and U-29 as balls and sticks. The carbonyl oxygen of Trp-44 forms a hydrogen bond with the C2’-hydroxy group of U-29. Electrostatic interactions like this hydrogen bonding cause a detectable quenching of the Trp-44 fluorescence signal.

Figure 5.1: Hydrogen bonding interaction between NFLD Trp-44 and U-29. Shown is a cartoon representation of the ThiITm•TPHE39A complex (PDB 4KR7) in a close-up view highlighting Trp-44 of the NFLD domain (red) and U-29 as balls and sticks. The carbonyl oxygen of Trp-44 forms a hydrogen bond with the C2’-hydroxy group of U-29. Electrostatic interactions like this hydrogen bonding cause a detectable quenching of the Trp-44 fluorescence signal.

5.2 Parameter Identifiability Analysis

All two-site binding models under study consisted of the global bm1to2.coop model and a subset of nested models being special cases of it. The bm1to2.coop model represents the fluorescence spectroscopic equivalent of the “sequential independent binding sites” model frequently used to analyze binding isotherms obtained from ITC measurements (Thordarson 2012). While in the context of ITC there are more constraints on the measurements reducing the noise level such as the control over the temperature and the usage of high-accuracy titration syringes, a relatively higher noise level was expected to be in the data of the quenching assay conducted in this work. The nls regression analysis already revealed that the bm1to2.coop model yielded unjustifiable parameter estimates. Against this backdrop, it had to be ensured whether this finding was due to the quality of the data set or whether specifically the more complex two-site binding models were overly complex and therefore too anticipated in order to describe the binding profile or even a combination of both. All this can be simultaneously addressed by a parameter identifiability analysis by means of maximum likelihood profiling of model parameter (4.1.2). Though computationally rather intensive, valuable information could be gained based on obtained likelihood profiles by just qualitatively analyzing the shape of the profile in the surroundings of its maximum with respect to the specified threshold of the confidence interval. In addition, maximum likelihood profiling has an established framework that is well documented (e.g., (Raue et al. 2010; Kreutz et al. 2013)) and was therefore the method of choice in this work.

When applied to the global bm1to2.coop model maximum likelihood profiling revealed that one of its parameter, namely \(K_1\), is structurally non-identifiable (cf. figure 4.2D, section 4.1.2) which manifests itself in flat likelihood profiles above the 95% confidence threshold lacking a unique maximum. When joint profile likelihoods are considered (cf. figure 4.2F), this results in parallel likelihood contours indicating that changing the value of \(K_1\) can be entirely compensated for by re-optimizing all other parameter of the binding model. That is why currently the titration data set does not provide enough information in order to separately and uniquely estimate this parameter so that this model is over-parameterized with respect to the applied measurement set up. This lack of identifiability can be solved in one of two ways. Firstly, the domain of the binding model can and should be further explored. The concentration range of the tRNA can be shifted towards higher concentrations where the tryptophane fluorescence signal is more stable. This will ultimately increase the signal to noise ratio while preserving the necessary amount of information in the data. The upper limit will be constrained by the absolute concentration of the tRNA that can be achieved after the purification process. In addition, the information content of the data set could be further increased by narrowing the dilution steps of the reactions so that more data points are included. Secondly, the bm1to2.coop model can be re-parameterized to increase the degrees of freedom. For this reason a subset of simpler two-site binding models all representing nested versions of the global bm1to2.coop model were included in the model set. When the proportionality factors are assumed to be additive and merged as in the bm1to2.coop.add model all parameter become structurally identifiable. However, \(K_2\) is still practically non-identifiable. This lack of identifiability can be solved by increasing the quality of the data set, i.e. by increasing the number of independently measured titration profiles. All further re-parameterizations providing yet simpler binding models such as, respectively, the bm1to2.deg and the bm1to2.deg.add models without exception resulted in uniquely identifiable parameter (not shown).

Maximum likelihood profiling revealed that the preliminary set up of the fluorescence spectroscopic titration experiment with an emphasis on the number of independent measurements and the number of data points causes parameter identifiability problems in the binding models describing cooperative binding of tRNA\(\mathrm{^{Phe}_{Bs}}\) substrates to ThiITm. To some extent this might impact the quality of the binding model set with respect to binding model evaluation by the \(AIC_c\) as well as further validation of the statistical robustness of these findings by cross-validation.

5.3 Kuhlback Leibler Divergence Analysis

An inevitable difficulty of the comparative analysis of variably parameterized binding models lies in the differences of the degrees of freedom that are linked to model complexity. As pointed out by Anderson and Burnham (2004) it is always possible to achieve a higher goodness-of-fit by just using a more complex model with more free parameter. For this reason, emphasis was put on finding an appropriate criterion that takes into account goodness-of-fit as quantified by the \(RMSE\) of the residuals but which is capable of neutralizing complexity differences between binding models under investigation. Original considerations included, inter alias, model evaluation on the basis of the extra sum of squares \(F\)-ratio test within the context of classical hypothesis testing. However, as stated in Motulsky and Christopoulos (2004) this approach should be applied to nested models only. With respect to the set of binding models considered in this work this would only be true for the Adair-Klotz-type binding models excluding the one-site binding model from further analysis. Additionally, this approach of binding model comparison requires the models to differ in at least one degree of freedom. Regarding the set of model functions considered here there is no clearly defined simplest binding model since the bm1to2.deg.add and the bm1to1 models, respectively, are mathematically equally simple with respect to the degrees of freedom. Another disadvantage of this approach is that it depends on rather subjective ad hoc decisions regarding the confidence level (\(\alpha\)-value) to choose and that there is currently no mathematical framework that links the \(\alpha\)-value to model complexity (Anderson, Burnham, and Thompson 2000). Therefore, this approach was not followed up on beyond theoretical considerations. An alternative approach to the null hypothesis testing provides model selection based on the \(AIC_c\) which efficiently and reliably trades off goodness-of-fit and model complexity. The great advantage of the \(AIC_c\) lies in its generality. It neither depends on having a clearly defined simplest binding model in the set nor do considered models need to be nested ones. In addition, it does not necessitate the definition of a significance level or threshold. In general, this approach can be thought of as a restricted regression. The idea is to evaluate the set of binding models by computing \(AIC_c\) scores and align the members on an ordinal scale with respect to the best describing top-ranked model based on \(AIC_c\) differences (\(\Delta_i\)). As suggested by Anderson and Burnham (2004) and frequently seen in literature (e.g. Dunlop, Franco, and Murray 2007; Snipes and Taylor 2014) it is convenient to re-express those differences in terms of evidence weights (\(\omega_i\)) as well as evidence ratios (\(\frac{\omega_i}{\omega_j}\)). In this way the relative strength of support for each of the binding models in the set could be assessed.

When the ranking of binding models based on respective \(AIC_c\) scores is compared with the ranking that would be achieved just based on the \(RMSE\) of the regression analysis it becomes readily apparent that the \(AIC_c\) efficiently penalizes the number of parameter (cf. table 4.2, 4.1.3). Interestingly, the bm1to2.deg.add model formerly ranked in third place according to the \(RMSE\) of the residuals is top-ranked based on the \(AIC_c\) score. Conversely, the known over-parameterized bm1to2.coop model is indeed efficiently penalized. It was the top-ranked binding model based on the \(RMSE\) of the residuals, whereas it is just ranked the fourth best binding model when the \(AIC_c\) score is the basis. With respect to the simpler bm1to2.coop.add and bm1to2.deg.add models that are ranked better this finding underpins that it is not justified to use a four-parameter model to describe the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) titration profile under the experimental set up currently applied. This is supported by the relatively large evidence ratio of \(\frac{\omega_{bm1to2.deg.add}}{\omega_{bm1to2.coop}} = 11.9\) stating that currently it is 11.9 times more likely that the data were generated by the bm1to2.deg.add model than by the global bm1to2.coop model. The bm1to2.deg.add model is the actual KL best model that most accurately captures the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) binding profile (\(\omega_i = 70.8 \%\)). As a rough rule of thumb that can be used at this stage of the analysis, competing models having a \(\Delta_i\) within 1-2 of the bm1to2.deg.add model also have substantial support of being the KL best model instead of the bm1to2.deg.add model, whereas binding models having a \(\Delta_i\) within about 4-7 have considerably less support of being the KL best model in the set (Anderson and Burnham 2004). Accodring to this categorization the bm1to2.coop.add model which is the most competing model in the set, is situated in between these ranges having a \(\Delta_i\) of 3.3. Therefore, it might remain elusive whether the KL divergence between this and the top-ranked binding model, respectively, translates into a strength of evidence large enough to convincingly provide support for a single best data describing model in the set when the \(AIC_c\) alone is used as the basis for evalutation. With respect to the binding process this means that no distinction might yet be made whether tRNA\(\mathrm{^{Phe}_{Bs}}\) binding to ThiITm would be best described as being solely statistical or as being positively cooperative. All other binding models received \(\Delta_i\) values at least greater than four and therefore do currently have considerable less support of being the data generating models instead of the top-ranked bm1to2.deg.add model. However, valuable information about the underlying binding stoichiometry of the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) complex could be gained based on consideration of the cumulative weights of evidence (CUSUM). All binding models ranked in places one to four comprise the bm1to2.deg.add, bm1to2.coop.add, bm1to2.deg, and the bm1to2.coop models, respectively. They unite a weight of evidence of 96.8% that the KL best model is among the two-site binding models as opposed to the one-site bm1to1 model having a weight of evidence lagging relatively far behind (3.2%).

In light of total ligand binding experiments as carried out in this work and considering the obtained ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) binding profile the set of binding models together are capturing the most information that can possibly be extracted from this kind of titration experiment, namely whether tRNA\(\mathrm{^{Phe}_{Bs}}\) binding to ThiITm occurs, respectively, in a 1:1 or 1:2 ratio and whether it is purely statistical (degenerate) or cooperative at the macroscopic scale with respect to the binding constants. A deeper analysis involving also the microscopic binding constants might not be justified since there is no change in the tryptophan fluorescence signal that would allow to distinguish the two binding sites of ThiITm in this experimental set up. Though statistically meaningful, obtained findings of the KL divergence analysis rely on residual sum of squares so that they are conditional on the data set based on which the nls regression analysis was carried out. To confirm that these findings are also statistically robust it was assessed how generalizable model performances are with respect to variations in the data as would inevitably be caused by future titration experiments. This was achieved by including a cross-validation as a second model evaluation criterion.

5.4 Quantifying Model Generalizability

As pointed out by Arlot and Celisse (2010) when nls regression analysis as well as model evaluation are conducted just based on the same data set this could lead to results that are overoptimistic with respect to the assessment of model performance. When it comes to model evaluation, just as important as yielding an adequate best fit curve describing the titration data at hand is the capability of binding models to predict yet unseen titration data. For this reason binding models were also evaluated via cross-validation (CV) on the basis of which it is possible to simulate future data sets. In this way, dependencies of binding model performances on variations in the data set can be uncovered. In this work complete 4-fold CV was applied. Firstly, the choice of four folds over five folds (or 10 folds) usually taken was a compromise. On the one hand, a split of a data set into five folds (or 10 folds) was shown to yield prediction error estimates that are statistically most robust (Fushiki 2011). On the other hand, such a partitioning would result into unequally sized folds when applied to the current setting (12 data points) in this work, namely into two folds containing three data points and three folds containing two data points. Having equally sized folds was given higher priority. Secondly, the random fold partitioning causes the estimated CV score to be a random variable whose variance can be quite high especially in small data sets. Repeated CV was shown to reduce this variance and thus making the estimate of the prediction error more stable (Jiang and Wang 2017). Furthermore, repeating it over all possible combinations of test-folds for every binding model, that is fitting 1200 different binding models in total, was computationally affordable in this case. Model selection was done based on a CV-score capturing the prediction accuracy as the average RMSE accross all 220 test-fold combinations. Most importantly, those results confirmed the findings of the KL divergence analysis with respect to the bm1to1 model. This model which is ranked the lowest based on its \(AIC_c\)-score is also the least predictive model in the set (\(\Delta RMSE\) = 3.44 \(\times\) 10-2, cf. figure 4.4B). That is, both evaluation metrics used in this work are in favour for a ThiITm•\(\mathrm{tRNA^{Phe}_{Bs}}\) complex formation that is based on a 1:2 stoichiometry. However, the model evaluation results differ in terms of which two-site binding model best describes the complexation. According to the complete CV analysis both cooperative models are the more predivtive than the degenerative models. The most predictive model is the known over-parameterized bm1to2.coop model (\(\Delta RMSE\) = 1.98 \(\times\) 10-2) which outcompetes the bm1to2.deg.add model which was top ranked based on its \(AIC_c\) score. Notably, the bm1to2.coop.add model containing a practically unidentifiable parameter (\(K_2\)) is the second best in each analysis. The reason for the predictive power of the cooperative binding models despite identifiability problems might be that in each case only one parameter is affected while all other parameters are still uniquely determined. Therefore, a potential noise capturing effect is not detectable. In light of these findings, a clear distinction between, respectively, purely statistical or coopereative binding of \(\mathrm{tRNA^{Phe}_{Bs}}\) to ThiITm has to await further analyses.

5.5 ANOVA

Model evaluation based on a complete 4-fold CV generated CV-score distributions for each binding model. This allowed to add a second dimension to this analysis by further investigating the central tendencies of the CV-score distributions by means of ANOVA. The resulting factorial design was based on binding model type as a factor containing five levels. For this reason the results of an omnibus ANOVA and corresponding effect sizes are based on four \(df\) and therefore are hardly interpretable. That is why a general linear model was formalised that is able to perform contrast analysis. The benefit of contrast analysis is that it allows to formulate focused hypotheses testing differences of groups of means that are based on one \(df\) while controlling the \(\alpha\)-level. In accordance with the CV analysis, a test of \(H_{01}\) which involves the key hypothesis of this work revealed that the one-site binding model bm1to1 has a statistically higher mean fold prediction error than the average two-site binding model (cf. tab 4.3 and fig 4.6B) which is determined by the negative direction of the contrast (\(\Psi_1\) = -9.6 \(\times\) 10-3). \(\Psi_1\) is part of the most important contrast-coded regression coefficient (\(\beta_{X_1}\)) of the ANOVA model due to the fact that the largest PRE is associated with it (\(\hat{\eta}^2\) = 9.42%). Further considering the sub-set of the two-site binding models it was found by testing \(H_{02}\) that the global bm1to2.coop model has a statistically lower mean CV-score than the average remaining two-site binding models as determined by the positive direction of the associated contrast (\(\Psi_2\) = \(\times\) 10-3). However, the corresponding contrast-coded regression coefficient (\(\beta_{X_2}\)) has a relatively low impact on the PRE (\(\hat{\eta}^2\) = 0.82%). Additionally computed interval estimates of the contrasts at the 95% confidence level were based on a common variance term amongst the factor levels (\(MSE_A\)) which reflects the assumtpion about homoscedastic distributed residuals of the general linear ANOVA model. A comparison of the \(CI_{95}\) associated with \(\Psi_1\) and \(\Psi_2\) revealed that the range of possible hypotheses outside these intervals including the respective zero hypotheses \(H_{01}\) and \(H_{02}\), that is, the contrasts equal zero, is larger for \(\Psi_1\) than for \(\Psi_2\). From this one might conclude that currently the data are more consistent with a clearer separation between one-site and two-site binding models (\(H_{01}\)) than with a further distinction of the two-site binding models (\(H_{02}\)). However, this is an indication based on a first experimental analysis. Ultimately, more biological replications are needed which will offer the possibility for meta-analysis of the intervals. Besides the generic ANOVA also a robust ANOVA was applied which should address the major concern that mean CV-scores are biased towards higher scores due to the presence of outlier (cf. figure 4.5A).

Following the generic ANOVA, involved statistics were determined conditional on a true null hypothesis which imposes some distributional assumptions on the residuals of the general linear model. For instance, the residuals should be homoscedastic and normally distributed. Also, there should be no additional source of error besides random measurement error. However, as already discussed in 5.1 there is an increasing measurement uncertainty in the titration data in relation to an increasing dilution factor of the tRNA. As a result, the data in the lowest concentration range including data points one and two introduce systematic error in the CV-ANOVA due to the fact that they are especially prevalent amongst the train fold combinations that lead to CV-score outlier. It is due to these outlier that the CV-score distributions are bimodal. Also, though not statistically investigated but phenomenologically observed, the bm1to1 model seems to be more affected by the presence of outlier than the two-site binding models. As a result of this departure from the normality assumption and due to the fact that the homoscedasticity assumption about the residuals still holds, a robust ANOVA was applied on the basis of trimmed CV-score distributions. Determined robust effect sizes of the contrasts (cf. tab 4.4) as well as additionally computed interval estimates (cf. fig 4.7) are in support of the corresponding estimates previously determined under the generic ANOVA approach.

6 Summary and Conclusions

In the course of this work a chimeric s4U8 synthesizing complex was investigated comprising T. maritima ThiI and B. subtilis tRNAPhe. Currently, the only structural information available about complexes of this system and in general is the crystal structure of the binary complex formed between ThiITm and TPHE39A which represents a truncated version of tRNAPhe still biologically active. These structures reveal that each monomer of dimeric ThiITm binds one TPHE39A substrate. However, it remains elusive whether this 1:2 stoichiometry can be extended to complexes of ThiITm and full-length tRNA substrates. For this reason, the stoichiometry of the complex was re-investigated in this work in terms of fluorescence spectroscopic binding studies using the interaction between ThiITm and full-length tRNA\(\mathrm{^{Phe}_{Bs}}\) as an example. Following this, an analytic work flow was established on the basis of which it is possible to infer the information about the stoichiometry from ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) binding profiles in the context of a multi-model approach by statistical comparison of a set of variably parameterized binding models all related to underlying 1:1 or 1:2 binding equilibria.

Starting with ordinary nls regression analysis the descriptive strength of each binding model in capturing the spectroscopic titration profile could be verified. Ater that, the model set was validated based on maximum likelihood profiling of model parameter. This important analysis step was used to monitor whether the sets of binding model specific parameter are uniquely identified by implicitly taking into account the noise level in the data, the domain restrictions of binding model functions due to lower and upper practical concentration limits, and model complexity. It revealed that the cooperative binding models show, respectively, structural and practical non-identifiability with respect to single parameter so that these models are currently not appropriately parameterized with respect to the information content in the data set. Binding model performance was then firstly evaluated by KL divergence analysis based on \(AIC_c\)-scores which take into account model complexity. The united weight of evidence of all two-site binding models determined during this analysis step provide relatively strong support that the actual KL best model is among this sub-set. Binding model performance was secondly evaluated based on CV-scores taking into account predictive power. The contrast analysis in the context of the generic ANOVA approach revealed that on average the one-site bm1to1 binding model yields statistically higher CV-scores than the two-site binding models which underpins the findings of the KL divergence analysis. The results of the generic ANOVA approach were subsequently supported by a robust ANOVA approach dealing with systematic measurement error induced departure of CV-score distributions from normality by using trimmed distributions.

Overall results from the complete analysis workflow provide evidence that the underlying stoichiometry of complexes between ThiI and its full-length tRNA substrates as determined in place of the chimeric ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) complex in the present work is indeed a 1:2 stoichiometry. The currently best approximating binding model of the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) interaction profile relating this stoichiometry to an appropriate underlying binding equilibrium is the bm1to2.deg.add model. This model describes binding of two tRNA\(\mathrm{^{Phe}_{Bs}}\) ligands to degenerate binding sites of ThiITm based on an overall stability constant \(\beta_{12}\). The bm1to2.deg.add model thereby meets the requirements of three quantifiable prime criteria featuring a good model: It is adequately descriptive and, taking into account model complexity, best fits the ThiITm•tRNA\(\mathrm{^{Phe}_{Bs}}\) binding profile in the simplest possible manner among the two-site binding models of the Adair-Klotz-type as well as it is comparatively generalizable. It therefore provides a solid current working hypothesis.

However, it has to be stressed that the results of this work are just first insights based on a preliminary experimental set up and data set. Ultimately, more biological replications are going to be needed based on an improved experimental set up concerning the applied concentration range, and the number of data points. This improvement is going to allow for a more robust direct comparison of binding models but also for a meta-analytic statistical analysis of CV-score distributions.

References

Akbar, S MD, K Sreeramulu, and Hari C Sharma. 2016. “Tryptophan Fluorescence Quenching as a Binding Assay to Monitor Protein Conformation Changes in the Membrane of Intact Mitochondria.” Journal of Bioenergetics and Biomembranes 48 (3): 241–47.

Allaire, JJ, Yihui Xie, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, Hadley Wickham, Joe Cheng, Winston Chang, and Richard Iannone. 2021. Rmarkdown: Dynamic Documents for R. https://CRAN.R-project.org/package=rmarkdown.

Anderson, David R, Kenneth P Burnham, and William L Thompson. 2000. “Null Hypothesis Testing: Problems, Prevalence, and an Alternative.” The Journal of Wildlife Management, 912–23.

Anderson, D, and K Burnham. 2004. “Model Selection and Multi-Model Inference.” Second. NY: Springer-Verlag 63 (2020): 10.

Arlot, Sylvain, and Alain Celisse. 2010. “A Survey of Cross-Validation Procedures for Model Selection.” Statistics Surveys 4: 40–79.

Christian, Thomas, Georges Lahoud, Cuiping Liu, and Ya-Ming Hou. 2010. “Control of Catalytic Cycle by a Pair of Analogous tRNA Modification Enzymes.” Journal of Molecular Biology 400 (2): 204–17.

Crain, Pamela F, and James A McCloskey. 1996. “The Rna Modification Database.” Nucleic Acids Research 24 (1): 98–99.

Dunlop, Mary J, Elisa Franco, and Richard M Murray. 2007. “A Multi-Model Approach to Identification of Biosynthetic Pathways.” In 2007 American Control Conference, 1600–1605. IEEE.

El Yacoubi, Basma, Marc Bailly, and Valérie de Crécy-Lagard. 2012. “Biosynthesis and Function of Posttranscriptional Modifications of Transfer Rnas.” Annual Review of Genetics 46: 69–95.

Elzhov, Timur V., Katharine M. Mullen, Andrej-Nikolai Spiess, and Ben Bolker. 2016. Minpack.lm: R Interface to the Levenberg-Marquardt Nonlinear Least-Squares Algorithm Found in Minpack, Plus Support for Bounds. https://CRAN.R-project.org/package=minpack.lm.

Fushiki, Tadayoshi. 2011. “Estimation of Prediction Error by Using K-Fold Cross-Validation.” Statistics and Computing 21 (2): 137–46.

Griffey, RH, DR Davis, Z Yamaizumi, S Nishimura, BL Hawkins, and CD Poulter. 1986. “15N-Labeled tRNA. Identification of 4-Thiouridine in Escherichia Coli tRNASer1 and tRNATyr2 by 1H-15N Two-Dimensional Nmr Spectroscopy.” Journal of Biological Chemistry 261 (26): 12074–8.

Henry, Lionel, and Hadley Wickham. 2019. Purrr: Functional Programming Tools. https://CRAN.R-project.org/package=purrr.

Jiang, Gaoxia, and Wenjian Wang. 2017. “Error Estimation Based on Variance Analysis of K-Fold Cross-Validation.” Pattern Recognition 69: 94–106.

Judd, Charles M, Gary H McClelland, and Carey S Ryan. 2017. Data Analysis: A Model Comparison Approach to Regression, Anova, and Beyond. Routledge.

Kassambara, Alboukadel. 2020. Ggpubr: ’Ggplot2’ Based Publication Ready Plots. https://CRAN.R-project.org/package=ggpubr.

Klotz, Irving M. 1946. “The Application of the Law of Mass Action to Binding by Proteins; Interactions with Calcium.” Archives of Biochemistry 9: 109–17.

Kreutz, Clemens, Andreas Raue, Daniel Kaschek, and Jens Timmer. 2013. “Profile Likelihood in Systems Biology.” The FEBS Journal 280 (11): 2564–71.

Lauhon, Charles T, Whitney M Erwin, and Giangthy N Ton. 2004. “Substrate Specificity for 4-Thiouridine Modification in Escherichia Coli.” Journal of Biological Chemistry 279 (22): 23022–9.

Mazerolle, Marc J. 2020. AICcmodavg: Model Selection and Multimodel Inference Based on (Q)AIC(c). https://CRAN.R-project.org/package=AICcmodavg.

Meschiari, Stefano. 2021. Latex2exp: Use Latex Expressions in Plots. https://CRAN.R-project.org/package=latex2exp.

Mezzasalma, Tara M, James K Kranz, Winnie Chan, Geoffrey T Struble, Celine Schalk-Hihi, Ingrid C Deckman, Barry A Springer, and Matthew J Todd. 2007. “Enhancing Recombinant Protein Quality and Yield by Protein Stability Profiling.” Journal of Biomolecular Screening 12 (3): 418–28.

Motulsky, Harvey, and Arthur Christopoulos. 2004. Fitting Models to Biological Data Using Linear and Nonlinear Regression: A Practical Guide to Curve Fitting. Oxford University Press.

Mount, John, and Nina Zumel. 2020a. Vtreat: A Statistically Sound ’Data.frame’ Processor/Conditioner. https://CRAN.R-project.org/package=vtreat.

———. 2020b. Wrapr: Wrap R Tools for Debugging and Parametric Programming. https://CRAN.R-project.org/package=wrapr.

Müller, Kirill, and Hadley Wickham. 2021. Tibble: Simple Data Frames. https://CRAN.R-project.org/package=tibble.

Neumann, Piotr, Kristina Lakomek, Peter-Thomas Naumann, Whitney M Erwin, Charles T Lauhon, and Ralf Ficner. 2014. “Crystal Structure of a 4-Thiouridine Synthetase–Rna Complex Reveals Specificity of tRNA U8 Modification.” Nucleic Acids Research 42 (10): 6673–85.

Raue, A, V Becker, U Klingmüller, and J Timmer. 2010. “Identifiability and Observability Analysis for Experimental Design in Nonlinear Dynamical Models.” Chaos: An Interdisciplinary Journal of Nonlinear Science 20 (4): 045105.

R Core Team. 2018. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Reinhard, Linda, Hubert Mayerhofer, Arie Geerlof, Jochen Mueller-Dieckmann, and Manfred S Weiss. 2013. “Optimization of Protein Buffer Cocktails Using Thermofluor.” Acta Crystallographica Section F: Structural Biology and Crystallization Communications 69 (2): 209–14.

Robinson, David, Alex Hayes, and Simon Couch. 2021. Broom: Convert Statistical Objects into Tidy Tibbles. https://CRAN.R-project.org/package=broom.

Shigi, Naoki. 2014. “Biosynthesis and Functions of Sulfur Modifications in tRNA.” Frontiers in Genetics 5: 67.

Snipes, Michael, and D Christopher Taylor. 2014. “Model Selection and Akaike Information Criteria: An Example from Wine Ratings and Prices.” Wine Economics and Policy 3 (1): 3–9.

Thiam, Khoudia, and Alain Favre. 1984. “Role of the Stringent Response in the Expression and Mechanism of Near-Ultraviolet Induced Growth Delay.” European Journal of Biochemistry 145 (1): 137–42.

Thordarson, Pall. 2012. “Binding Constants and Their Measurement.” Supramolecular Chemistry: From Molecules to Nanomaterials.

Traxler, Matthew F, Sean M Summers, Huyen-Tran Nguyen, Vineetha M Zacharia, G Aaron Hightower, Joel T Smith, and Tyrrell Conway. 2008. “The Global, ppGpp-Mediated Stringent Response to Amino Acid Starvation in Escherichia Coli.” Molecular Microbiology 68 (5): 1128–48.

Wickham, Hadley, and Jennifer Bryan. 2019. Readxl: Read Excel Files. https://CRAN.R-project.org/package=readxl.

Wickham, Hadley, Winston Chang, Lionel Henry, Thomas Lin Pedersen, Kohske Takahashi, Claus Wilke, Kara Woo, Hiroaki Yutani, and Dewey Dunnington. 2020. Ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. https://CRAN.R-project.org/package=ggplot2.

Wickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2021. Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.

Wickham, Hadley, and Lionel Henry. 2020. Tidyr: Tidy Messy Data. https://CRAN.R-project.org/package=tidyr.

Wilcox, Rand R. 2011. Introduction to Robust Estimation and Hypothesis Testing. Academic press.

Xiao, Nan. 2018. Ggsci: Scientific Journal and Sci-Fi Themed Color Palettes for ’Ggplot2’. https://CRAN.R-project.org/package=ggsci.

Xie, Yihui. 2020. Knitr: A General-Purpose Package for Dynamic Report Generation in R. https://CRAN.R-project.org/package=knitr.

———. 2021. Bookdown: Authoring Books and Technical Documents with R Markdown. https://CRAN.R-project.org/package=bookdown.

Yuen, Karen K. 1974. “The Two-Sample Trimmed T for Unequal Population Variances.” Biometrika 61 (1): 165–70.

Zhu, Hao. 2021. KableExtra: Construct Complex Table with ’Kable’ and Pipe Syntax. https://CRAN.R-project.org/package=kableExtra.