%% BioMed_Central_Tex_Template_v1.06
%%                                      %
%  bmc_article.tex            ver: 1.06 %
%                                       %

%%IMPORTANT: do not delete the first line of this template
%%It must be present to enable the BMC Submission system to
%%recognise this template!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                     %%
%%  LaTeX template for BioMed Central  %%
%%     journal article submissions     %%
%%                                     %%
%%          <8 June 2012>              %%
%%                                     %%
%%                                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                 %%
%% For instructions on how to fill out this Tex template           %%
%% document please refer to Readme.html and the instructions for   %%
%% authors page on the biomed central website                      %%
%% http://www.biomedcentral.com/info/authors/                      %%
%%                                                                 %%
%% Please do not use \input{...} to include other tex files.       %%
%% Submit your LaTeX manuscript as one .tex document.              %%
%%                                                                 %%
%% All additional figures and files should be attached             %%
%% separately and not embedded in the \TeX\ document itself.       %%
%%                                                                 %%
%% BioMed Central currently use the MikTex distribution of         %%
%% TeX for Windows) of TeX and LaTeX.  This is available from      %%
%% http://www.miktex.org                                           %%
%%                                                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% additional documentclass options:
%  [doublespacing]
%  [linenumbers]   - put the line numbers on margins

%%% loading packages, author definitions

%\documentclass[twocolumn]{bmcart}% uncomment this for twocolumn layout and comment line below
\documentclass{bmcart}

%%% Load packages
%\usepackage{amsthm,amsmath}
%\RequirePackage{natbib}
%\RequirePackage[authoryear]{natbib}% uncomment this for author-year bibliography
\usepackage[utf8]{inputenc} %unicode support
%\usepackage[applemac]{inputenc} %applemac support if unicode package fails
%\usepackage[latin1]{inputenc} %UNIX support if unicode package fails
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{booktabs}
\usepackage{float}
\usepackage{color, colortbl}
\usepackage{mathtools}
\usepackage{enumerate}
\usepackage[obeyspaces]{url}
%\usepackage[svgnames]{xcolor}
\usepackage{array}
\usepackage{caption}
\usepackage{subcaption}
\graphicspath{ {img/} }
\usepackage{tabularx}
\usepackage[normalem]{ulem} % for \sout
\usepackage{multirow}
\usepackage{pdflscape}
\usepackage{chngpage}
\usepackage{xr}
\usepackage{adjustbox}
\usepackage{acro}
%\RequirePackage{hyperref}
\usepackage{hyperref}
\usepackage{cleveref}

%\externaldocument{SI}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                             %%
%%  If you wish to display your graphics for   %%
%%  your own use using includegraphic or       %%
%%  includegraphics, then comment out the      %%
%%  following two lines of code.               %%
%%  NB: These line *must* be included when     %%
%%  submitting to BMC.                         %%
%%  All figure files must be submitted as      %%
%%  separate graphics through the BMC          %%
%%  submission process, not included in the    %%
%%  submitted article.                         %%
%%                                             %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\def\includegraphic{}
\def\includegraphics{}

\newcommand{\DG}{\Delta G}
\newcommand{\DDG}{\Delta\Delta G}
\newcommand{\IC}{IC_{50}}
\newcommand{\RF}{RF}


%%% Put your definitions there:
\startlocaldefs
\endlocaldefs

\DeclareAcronym{HIV}{
  short = HIV-1,
  long = human immunodeficiency virus type 1,
  class = abbrev
}
\DeclareAcronym{RAM}{
  short = RAM,
  long = resistance associated mutation,
  class = abbrev
}
\DeclareAcronym{PI}{
  short = PI,
  long = protease inhibitor,
  class = abbrev
}
\DeclareAcronym{MD}{
  short = MD,
  long = molecular dynamics,
  class = abbrev
}
\DeclareAcronym{ART}{
  short = ART,
  long = antiretroviral therapy,
  class = abbrev
}
\DeclareAcronym{AIDS}{
  short = AIDS,
  long = acquired immune deficiency syndrome,
  class = abbrev
}
\DeclareAcronym{APV}{
  short = APV,
  long = Amprenavir,
  class = abbrev
}
\DeclareAcronym{IDV}{
  short = IDV,
  long = Indinavir,
  class = abbrev
}
\DeclareAcronym{DRV}{
  short = DRV,
  long = Darunavir,
  class = abbrev
}
\DeclareAcronym{LPV}{
  short = LPV,
  long = Lopinavir,
  class = abbrev
}
\DeclareAcronym{ATV}{
  short = ATV,
  long = Atazanavir,
  class = abbrev
}
\DeclareAcronym{SQV}{
  short = SQV,
  long = Saquinavir,
  class = abbrev
}
\DeclareAcronym{NFV}{
  short = NFV,
  long = Nelfinavir,
  class = abbrev
}
\DeclareAcronym{FPV}{
  short = FPV,
  long = Fosamprenavir,,
  class = abbrev
}
\DeclareAcronym{RF}{
  short = RF,
  long = resistance factor,
  class = abbrev
}
\DeclareAcronym{FMA}{
  short = FMA,
  long = functional mode analysis,
  class = abbrev
}
\DeclareAcronym{PLS}{
  short = PLS,
  long = partial least-squares,
  class = abbrev
}

%%% Begin ...
\begin{document}

%%% Start of article front matter
\begin{frontmatter}

\begin{fmbox}
\dochead{Research}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the title of your article here     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{Non-Active Site Mutants of HIV-1 Protease Influence Resistance and Sensitization Towards Protease Inhibitors}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors here                   %%
%%                                          %%
%% Specify information, if available,       %%
%% in the form:                             %%
%%   <key>={<id1>,<id2>}                    %%
%%   <key>=                                 %%
%% Comment or delete the keys which are     %%
%% not used. Repeat \author command as much %%
%% as required.                             %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\author[
   addressref={MPII,SBG,AZ},                   % id's of addresses, e.g. {aff1,aff2}
]{\inits{TB}\fnm{Tomas} \snm{Bastys}}
\author[
   addressref={MPIB}
]{\inits{VG}\fnm{Vytautas} \snm{Gapsys}}
\author[
   addressref={ERL}
]{\inits{HW}\fnm{Hauke} \snm{Walter}}
\author[
   addressref={UKK}
]{\inits{EH}\fnm{Eva} \snm{Heger}}
\author[
    addressref={MPII,UK}
]{\inits{NTD}\fnm{Nadezhda T.} \snm{Doncheva}}
\author[
   addressref={UKK}
]{\inits{RK}\fnm{Rolf} \snm{Kaiser}}
\author[
   addressref={MPIB}
]{\inits{BLdG}\fnm{Bert L.} \snm{de Groot}}
\author[
    addressref={HIPS,SU,MPII},
   corref={HIPS},                       % id of corresponding address, if any
   email={Olga.Kalinina@helmholtz-hips.de}
]{\inits{OVK}\fnm{Olga V.} \snm{Kalinina}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors' addresses here        %%
%%                                          %%
%% Repeat \address commands as much as      %%
%% required.                                %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\address[id=MPII]{\orgname{Department for Computational Biology and Applied Algorithmics, Max Planck Institute for Informatics}, \street{D-66123 Saarbr{\"u}cken, Germany}}
\address[id=SBG]{\orgname{Saarbr{\"u}cken Graduate School of Computer Science, University of Saarland} \street{D-66123 Saarbr{\"u}cken, Germany}}
\address[id=AZ]{\orgname{Currently at: Discovery Sciences, AstraZeneca IMED Biotech Unit} \street{SE-431 83 M{\"o}lndal, Sweden}}
\address[id=MPIB]{\orgname{Computational Biomolecular Dynamics Group, Department of Theoretical and Computational Biophysics, Max Planck Institute for Biophysical Chemistry} \street{D-37077 G{\"o}ttingen, Germany}}
\address[id=ERL]{\orgname{Institute for Virology, University of Erlangen-Nuremberg} \street{D-91054 Erlangen, Germany}}
\address[id=UKK]{\orgname{Institute of Virology, University of Cologne} \street{D-50935 K{\"o}ln, Germany}}
\address[id=UK]{\orgname{Faculty of Health and Medical Sciences, University of Copenhagen} \street{DK-2200 Copenhagen, Denmark}}
\address[id=HIPS]{\orgname{Helmholtz Institute for Pharmaceutical Research Saarland, Helmholtz Centre for Infection Research} \street{D-66123 Saarbr{\"u}cken, Germany}}
\address[id=SU]{\orgname{Medical Faculty, Saarland University} \street{D-66421 Homburg, Germany}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter short notes here                   %%
%%                                          %%
%% Short notes will be after addresses      %%
%% on first page.                           %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{fmbox}% comment this for two column layout

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Abstract begins here                 %%
%%                                          %%
%% Please refer to the Instructions for     %%
%% authors on http://www.biomedcentral.com  %%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstractbox}

\begin{abstract}
    \subsection*{\textbf{Background}}

    HIV-1 can develop resistance to antiretroviral drugs, mainly through mutations within the target regions of the drugs.
    In HIV-1 protease, a majority of resistance-associated mutations
    that develop in response to the therapy with protease inhibitors
    are found in the protease's active site that serves also as a binding pocket
    for the protease inhibitors, thus directly impacting the protease-inhibitor
    interactions.
    Some resistance-associated mutations, however, are found in more distant regions, and the exact mechanisms
    how these mutations affect protease-inhibitor interactions are unclear.
    Furthermore, some of these mutations, e.g. N88S and L76V,
    do not only induce resistance to the currently administered drugs,
    but contrarily induce sensitivity towards other drugs.
    In this study, mutations N88S and L76V, along with two other resistance-associated mutations,
    M46I and I84V, are analysed by means of molecular dynamics simulations
    to investigate their role in complexes of the protease with different
    inhibitors and in different background sequence contexts.

    \subsection*{\textbf{Results}}

    Using these simulations for alchemical calculations to estimate the effects of mutations
    M46I, I84V, N88S, and L76V on binding free energies
    shows they are in general in line with the mutations' effect on $IC_{50}$ values.
    For the primary mutation L76V, however, the presence of a background mutation M46I
    in our analysis influences whether the unfavourable effect of L76V on inhibitor binding
    is sufficient to outweigh the accompanying reduction in
    catalytic activity of the protease.
    Finally, we show that L76V and N88S changes the hydrogen bond stability of these residues
    with residues D30/K45 and D30/T31/T74, respectively.

    \subsection*{\textbf{Conclusions}}

    We demonstrate that estimating the effect of both binding pocket and distant mutations
    on inhibitor binding free energy using alchemical calculations can reproduce
    their effect on the experimentally measured $IC_{50}$ values.
    We show that distant site mutations L76V and N88S affect the hydrogen bond network
    in the protease's active site, which offers an explanation for the
    indirect effect of these mutations on inhibitor binding.
    This work thus provides valuable insights on interplay between primary and background mutations
    and mechanisms how they affect inhibitor binding.

\end{abstract}

\begin{keyword}
\kwd{Alchemical binding free energy change calculation}
\kwd{Distant site mutations}
\kwd{HIV-1 protease inhibitors}
\kwd{Hydrogen bond network perturbation}
\kwd{Resistance-associated mutations}
\end{keyword}

\end{abstractbox}
%
%\end{fmbox}% uncomment this for twcolumn layout

\end{frontmatter}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Main Body begins here                %%
%%                                          %%
%% Please refer to the instructions for     %%
%% authors on:                              %%
%% http://www.biomedcentral.com/info/authors%%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%% See the Results and Discussion section   %%
%% for details on how to create sub-sections%%
%%                                          %%
%% use \cite{...} to cite references        %%
%%  \cite{koon} and                         %%
%%  \cite{oreg,khar,zvai,xjon,schn,pond}    %%
%%  \nocite{smith,marg,hunn,advi,koha,mouse}%%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%% start of article main body
% <put your article body there>

%%%%%%%%%%%%%%%%
%% Background %%
%%

\section*{Background}
With around 36.7 million people already infected and 1.8 million being newly infected per year,
the \ac{HIV} (further HIV) remains a global epidemic \cite{UNAIDS17a}. 
Since more than half of the infected individuals receive \ac{ART} \cite{UNAIDS17b},
\ac{AIDS}-related deaths have dropped to 1 million per year \cite{UNAIDS17a}.
For those under treatment, resistance towards drugs is a major cause for the need for switching of the therapy.

HIV protease (Figure \ref{fig:PR_structure}),
a protein responsible for cleaving HIV polyproteins,
is one of the major targets for ART, and protease inhibitors
are currently recommended as second- or third-line ART treatments \cite{WHO16}.
PIs are competitive binders of the protease,
occupying the active site of the protein once bound.
In line with this, most of the major resistance-associated mutations towards PIs
appear in the different structural elements composing the active site pocket,
such as in the active site loop (residues D30, V32, and L33),
the so-called 80s loop (residues V82 and I84)
which together form the sides of the pocket,
and the flap region of the protease
(residues M46, I47, G48, I50, and I54, Figure \ref{fig:PR_structure}).
Yet several RAMs are also found in distant to the binding pocket sites,
e.g.~in the amino acids N88 and L90 in the protease's $\alpha$-helix
or L76 in the protein's hydrophobic core.
The effect of these mutations on inhibitor binding is likely to be not through direct interactions.

Some mutations exhibit opposite effects on binding of some PIs, e.g.~L76V is
associated with resistance towards \ac{APV}, \ac{IDV},
\ac{DRV}, and \ac{LPV},
but increases sensitivity towards \ac{ATV} and \ac{SQV} \cite{Young10,Wiesmann11}.
Similarly, N88S is a RAM towards IDV and \ac{NFV}, but increases
susceptibility towards APV \cite{Ziermann00,Resch02,Vermeiren07}
or its prodrug \ac{FPV} \cite{Rhee10}.

Numerous studies have addressed the molecular effects of RAMs.
Some studies analyse the effects of selected major RAMs 
on binding of different inhibitors \cite{Hou07,Muzammil07,Alcaro09,Shen10,Mittal13,Duan15,Yu15,Bastys18}
and others the effect of different RAMs on binding of the same inhibitor
\cite{Liu05, Kovalevsky06, Kozisek07, Hu09, Cai10, Chen10, Kar12, Meher12, Ragland14, Ragland17}.
Most of the studies are however focused on a single mutation-inhibitor combination,
particularly for major RAMs outside of the binding pocket,
and thus offer only a limited perspective on molecular mechanisms of the protease resistance.
To the best of our knowledge,
the mechanism of action of a mutation on binding of different inhibitors
has not been investigated in the aforementioned cases,
where the same mutation is known to cause resistance to certain inhibitors,
while making the protein sensitive towards other inhibitors \cite{Young10,Wiesmann11,Ziermann00,Resch02,Vermeiren07,Rhee10}.
In clinical practice of HIV treatment,
cases like this potentially provide an opportunity for combined treatment:
a combination of PIs that are associated with opposite effects
of a particular mutation put the virus in a situation where either mutation variant
in this position will render the protease susceptible to one of the drugs
\cite{Larder95, Ziermann00, Wiesmann11}.
Understanding the underlying molecular phenomena could potentially provide
important insights into HIV inhibitor resistance, as well as a possibility
to transfer this knowledge to treatment of other viruses, and for inhibitor design.

In an experimental setting, resistance of mutant proteins towards PIs,
such as in the studies above,
is typically measured in terms of $IC_{50}$ (concentration required to inhibit viral activity by 50\%).
Thus, the ratio between $IC_{50}$ in mutant and the same measurement for the
wildtype protease (typically with the consensus sequence from the strain HXB2),
also called \ac{RF},
is a useful descriptor for resistance of different mutated proteins.
RF is directly related to the free energy of inhibitor binding, $\DG$,
and the protease enzymatic activity, $K_m$ \cite{Cheng73}.

We have previously shown that the effect of mutations in the HIV protease
on inhibitor binding, $\DDG$, can be accurately predicted \emph{in silico}
using alchemical methods based on MD simulations \cite{Bastys18}.
Additionally, analysis of the underlying trajectories from the MD simulations can
also reveal the mechanisms underlying the effects of mutations on protein-inhibitor interactions,
protein structural changes, as well as on their altered dynamics.

In this study, we apply alchemical calculations based on MD simulations
to estimate the effect of RAMs M46I, V84I, N88S, and L76V,
both in the binding pocket and outside of it,
on the binding of PIs APV, IDV, LPV, and SQV:
in total 17 different PI-mutation combinations.
We demonstrate that we can faithfully reproduce the previously reported
effects of the first three mutations on the resistance,
including previously investigated sensitizing and desensitizing effects
of the mutation N88S on binding of APV and IDV \cite{Ziermann00,Vermeiren07}.
We show that changes of the hydrogen bonding network of the protease
that involve D30, located in the active site pocket,
can explain these effects.
For L76V, located in the protease's hydrophobic core,
we investigate its effect on binding of ATV, IDV, LPV, and SQV
in different clinically relevant sequence contexts.
In the phenotypic assays we consistently observe
a resensitizing effect of this mutation towards
ATV and SQV as well as its resistance-associated effect for IDV and LPV,
which was reported previously \cite{Young10,Wiesmann11},
and we generally reproduce these effects computationally.
We suggest a mechanistic explanation of this effect, demonstrating that
these mutations affect the arrangement of residues around the binding
pocket of the protease, which, in turn, similarly to N88S,
affects the hydrogen bonding network of D30,
directly influencing inhibitor binding.

\section*{Results and Discussion}
\subsection*{Estimation of Resistance Factors from the Change in the Inhibitor Binding Free Energy}

We aimed to assess whether we can reproduce the ratios of experimentally measured
resistance factors ($\RF_R = \RF_{mutant}/\RF_{wildtype}$)
between the wildtype and mutant proteins by estimating the change in free energy
of inhibitor binding upon mutations in the protease using MD simulation with
alchemical methods \cite{Gapsys15b}.
For this purpose we selected a dataset of sixteen complexes
from HIVdb database \cite{HIVdb} for which resistance factors towards inhibitors
FPV, IDV, LPV, and SQV had been measured experimentally (Tables \ref{tab:HIVdb_RF} and S1).
These complexes could be paired amongst each other
such that the RF has been measured for the same inhibitor and the same protease strain
with and without the mutation, namely:
IDV and FPV with mutation M46I;
FPV, IDV, LPV, and SQV with mutation I84V;
FPV and IDV with mutation N88S.
This pairing allowed calculating relative resistance factors, $\RF_R$ (Table \ref{tab:RF_ddG}).
Moreover, for all of these mutations
wildtype and mutant protein background sequences were identical
across the different inhibitors,
with the exception of N88S, where complexes with FPV had a background mutation L77I
and complexes with IDV had a R57G background mutation present.
Both of R57G and L77I mutations are found next to each other on two parallel beta-sheets at a distant site of the protease.
Although they are close to the resistance-associated positions 76 and 54, respectively,
unlike for those residues, the side chains of residues 77 and 57 are pointing away from the protease binding pocket.
Nevertheless, the mutation R57G has been suggested to be a protease-inactivating mutation \cite{Ott03}.
L77I on the other hand, has been reported to be a compensatory mutation for I84V, restoring protein stability \cite{Chang11}.
M46I, I84V, and N88S are all considered to be major RAMs
against the corresponding inhibitors as reported in the HIVdb.
The only exception is N88S for which in combination with APV/FPV
sensitivity has been previously reported \cite{Ziermann00,Resch02,Vermeiren07,Rhee10}.
Accordingly, for all of these pairs, protease with a RAM
had a higher RF than the wildtype, while N88S reduced resistance towards FPV (Table \ref{tab:HIVdb_RF}).

We first estimated the effect of the mutation on the free energy
of inhibitor binding, $\DDG$.
For this purpose, we performed MD simulations for the
protease mutants and wildtype in inhibitor-bound and unbound state,
where in the bound state the simulations were performed taking into account both
alternative protonation states of the catalytic residues of the active site,
D25 and D25$^\prime$, to account for asymmetry of this complex.
This allows us to identify which protonation state is more likely for both wildtype and mutant complexes,
as well as to increase the accuracy of the $\DDG$ estimation, as we reported previously \cite{Bastys18}.
The resulting $\DDG$ calculations (Tables \ref{tab:RF_ddG} and S2) overall indicated
a good agreement in discriminating resistant and sensitizing effects of mutations on the protein-ligand binding,
including the opposite effects of N88S towards IDV and APV.
An exception to that is M46I, where the mutation had a modest effect on $\DG$
which was within the estimated error range.
The mutation of this flap residue, whose side-chain points away from the protease binding pocket,
has been associated with resistance towards different PIs,
but it typically appears in combination with other RAMs and has been suggested
to compensate the decreased catalytic activity of mutant proteases
\cite{Ho94,Pazhanisamy96,Schock96,Nijhuis09,Louis11,Henderson12}.

The possibility of mutations having an effect on the catalytic activity of the enzyme, $K_m$,
precludes direct comparison of the $\DDG$ estimates of mutation effects
on inhibitor binding and the $\RF_R$ corresponding to that mutation.
In previous studies of resistance mutations of another enzyme of HIV, reverse transcriptase,
$\DDG$ was considered to approximate changes in $IC_{50}$ \cite{Rizzo00, Wang01, Udier03}.
This is at odds with the fact that the majority of the mutations for which this approximation was used,
such as L100I, V106A, and Y188L, although not located directly in the active site,
have been previously reported to affect the catalytic potential of the enzyme
\cite{Loya94, Maga97, Archer00, Locatelli04, Das05}.
RAMs affecting the catalytic activity have also been reported for HIV protease
\cite{Chen95, Pazhanisamy96, Schock96, Nijhuis99}.
In the present study, similarly to the aforementioned computational studies of the reverse transcriptase,
only the experimentally measured $IC_{50}$ values for mutations are available
for the enzyme and different inhibitors.
To account for the changes in the binding free energy and $K_m$,
we developed a Bayesian method which combines multiple experimental $\RF_R$ measurements
and $\DDG$ estimates to calculate $\RF_R$ (see Methods).

We then compared the estimated $\RF_R$ values to their experimental measurements (Table \ref{tab:RF_ddG}, Figure \ref{fig:RF}).
The increase of resistance towards inhibitors ($\RF_R>1$) was correctly predicted for all
mutations, as was the sensitizing effect of N88S towards APV ($\RF_R<1$),
and the experimental $\RF_R$ values were within the
corresponding calculated distributions based on the $\DDG$ estimates (Figure S1).

While the structural effect of N88S mutation has been previously analysed for NFV,
to the best of our knowledge, the opposite effects of this mutation on susceptibility
towards APV and IDV have not been previously addressed.
It has been suggested that substitution of asparagine with serine
creates a hydrogen bond with the residue D30, which in turn affects
the interaction between D30 and the inhibitor NFV \cite{Ode07}.
A similar observation with regard to the N88S effect on the interaction with D30
has been previously made for L10F/N88S mutant with NFV as well as the unliganded
N88S protease \cite{Vasavi17, Bihani09}.
Another mutation at this site that occurs in patients treated with NFV, N88D,
has been reported to co-occur with mutation D30N,
which coincides with losing water molecules that mediate this site's interactions
with residues T31 and T74 \cite{Mahalingam02}.
Seeking to verify whether these effects extend to complexes of N88S
with APV and IDV, we performed hydrogen bond network analysis of mutant
and wildtype complexes, where we measured the average number of hydrogen bonds over the course of the simulations.
Indeed, similar effects were confirmed:
S88 formed a hydrogen bond with D30 more frequently than N88,
whereas N88 formed a hydrogen bond more frequenly to T31 and T74
compared to S88 (Table S3).

\subsection*{Estimating Resistance Factor for Clinical Samples with the L76V Mutation}

Among patients in Germany undergoing HIV treatment,
individuals who underwent multiple therapy failures against different PI-regimens were identified.
Among the group of patients we described in the manuscript of Wiesman et al.~\cite{Wiesmann11},
there were viral variants which displayed different resistance-determining mutations to ATV, SQV, IDV, and LPV.
The extraordinary observation was that a specific amino acid change L76V increased resistance
to LPV and IDV, while at the same time giving a clinically relevant re-sensitisation to SQV and ATV.
Those variants were observed in the diagnostic procedure, sequenced, and subsequently tested in a phenotypic assay.
We analysed variants which showed some of the highest changes in RF upon mutation
to evaluate whether we can computationally estimate RF values
as measured in the phenotypic assay (Table \ref{tab:pr:clinical_RF}).

Just as for mutations M46I, I84V, and N88S, for the mutation L76V multiple RF measurements
for different inhibitors in the same sequence context were available,
thus enabling us to computationally predict the $\RF_R$ values.
The sequences of the protease complexes analysed
had a large number of background mutations accumulated compared to the reference sequence HXB2,
making it difficult to find complexes in the PDB with sequences matching the studied genotypes.
Thus in the protein modelling stage between 11 and 19 mutations
had to be introduced to create protein models
with sequences corresponding to those for which $\RF_R$ was measured
(see Methods).
Including the target mutation L76V as well meant that up to 20\% of protease residues
had to be modelled \emph{in silico}.

First, we estimated the effect of the mutation L76V on inhibitor binding
in terms of the change of the binding free energy $\DDG$ (Tables \ref{tab:ddG_L76V} and S4).
The increase of the binding free energy, corresponding to the decrease in inhibitor affinity,
was predicted for all complexes where mutations were observed to increase the protease RF
(RU1 with LPV, iZ2 with IDV and LPV).
The decrease of RF, on the other hand,
did not always correspond to a negative value of $\DDG$:
L76V was predicted to increase
the affinity of inhibitor binding for inhibitors ATV and SQV
only for the genotype GH9, but not for the genotype FB15,
nor for inhibitor ATV in the context of the genotypes RU1 or iZ2.
The genotypes FB15 and iZ2 lack the background mutation M46I
(the former being wildtype at that position and the latter having mutation M46L),
which has been suggested to co-occur with L76V to compensate for its compromising
effect on the replication capacity of HIV \cite{Nijhuis09,Louis11}.
This suggests that in this case the dominant effect of the mutation L76V might be
exerted through decreasing the protease's catalytic activity $K_m$.
However for RU1, which has the M46I mutation, we cannot explain the positive $\DDG$ value
for complexes with ATV using this argument.

The $\DDG$ estimates were used to calculate $\RF_R$ in the same fashion
as for mutations M46I, I84V, and N88S (Table \ref{tab:ddG_L76V}, Figure \ref{fig:RF_L76V}).
For most of the complexes we correctly predicted whether the mutation made the protein
more resistant or more sensitive towards the inhibitor.
This included the prediction of sensitizing effect of mutation in the genotype FB15
for both ATV and SQV for which the inhibitor affinity increased based on the $\DDG$ estimates.
Sensitisation towards ATV was, on the other hand, not observed for genotypes RU1 and iZ2.
The experimental $\RF_R$ values were however within the corresponding
calculated distributions based on the $\DDG$ estimates (Figure S2).

\subsection*{Effect of L76V on Molecular Interactions}

Next, we focused on the structural changes in the protease upon the mutation L76V,
for which purpose we first analysed the hydrogen bond network of the protein.
It was consistently seen across all of the different genotypes
that the mutation L76V increases the probability of observing a
hydrogen bond between residues D30 and K45 (Figure \ref{fig:pr:hbond} and Table \ref{tab:hbond_L76V}).
A previous study found a highly significant correlation between mutations in
these sites \cite{Margerison08}, potentially indicating the importance of the
interaction between these two oppositely charged residues.
Seeking whether this was a result of side chain rearrangement,
we performed \ac{FMA} based on \ac{PLS} regression,
a supervised machine learning technique which builds models that
distinguish behaviour in MD trajectories between the wildtype and mutant protease complex
based on their Cartesian atoms coordinates.
These models are interpretable in terms of protein conformational changes associated with the mutation.
In analysing FMA models, we could see that the mutation L76V
caused a tendency of the side chains of residues D30, K45, and Q/E58
to shift towards the binding pocket (Figures \ref{fig:pr:hbond} and S3).
This shift is likely the result of fine rearrangement of residues
in the region as a consequence of a larger side chain
of leucine being replaced by a smaller valine.
The effect of L76V on D30--K45 hydrogen bond is thus similar to
the effect of the other distant site mutation we analysed in this study, N88S,
on hydrogen bonding of D30.
Mutations at site N88 have been reported to be correlated
to mutations of both D30 and K45 \cite{Margerison08}.
Displacement of Q/E58, on the other hand, is in line with the previously reported
co-occurance of mutations L76V and Q58E \cite{Champenois11},
and both of these mutations were found in the patient sample RU1.

We calculated direct protein-inhibitor interaction energies to see whether
the L76V mutation impacted direct interactions of D30, K45,
or other residues with the inhibitors (Figures \ref{fig:energy} and S4).
Indeed, changes in the interaction of D30/D30$^\prime$ with the inhibitors are in general
in line with changes of $\RF_R$:
negative, or favourable, interaction energy values correspond to $\RF_R<1$,
and positive, or unfavourable, interaction energies correspond to $\RF_R>1$.
An exceptions to that is the protease of the iZ2 genotype in complex with IDV,
where a favourable effect on the direct interaction energy of D30/D30$^\prime$
with the inhibitor is observed (Figures \ref{fig:energy} and S4),
and the same enzyme in complex with LPV,
where no notable effect on this interaction is seen
(Figure \ref{fig:energy} and S4).
The effect of L76V on interaction energies between the inhibitors and K45/K45$^\prime$
was, on the other hand, modest (Figure \ref{fig:energy} and S4).
A number of other residues' direct interactions with inhibitor were affected.
These residues are widely distributed across most of the active site pocket,
including the active site loop (including the D30), flap (including the K45), and 80s loop regions.
This is in line with our observations from a previous study of the effect
of mutations I50V, G48V, and L90M on protein-inhibitor interactions,
where interactions of residues in these regions were also affected
by the mutations \cite{Bastys18}.
Interestingly, measurable differences were observed for interactions of the residue at position 76
with the inhibitor for complexes of genotype FB15 with LPV and ATV and genotype iZ2 with IDV.
But given that side chain of residue 76 has minimal exposure to the binding pocket,
those differences are negligible.

Overall, these results indicate a pathway for how the mutation L76V impacts the inhibitor binding
through altering the interactions of other residues with the inhibitor without actual mutations at these sites.
A similar observation has been previously made for another pair of an active
site loop and distant site mutations, namely for mutation L90M
that alters the interactions of the residue at this position with D25,
which in turn affects the interactions of D25 with the inhibitor as well as
with other residues in the binding pocket \cite{Hong00, Ode06, Bastys18}.
Our observations on energetic and structural consequences of the mutation L76V are also in line with
its previously reported effects on the ligand binding affinity
for the inhibitor DRV through both changes in protein-inhibitor interactions and
changes in the inter-residue distances in the binding pocket \cite{Ragland14}.

Recently, a study reported experimentally resolved structures
of the wildtype and the L76V mutant of the HIV protease in complex with inhibitors DRV, LPV, TPV,
as well as with two experimental compounds, GRL-0519 and GRL-5010 \cite{Sam18}.
It has been observed that mutation does not change the backbone structure of the protease,
however residue 76 loses contacts with D30 and T74,
and, for structures with LPV, there is a slight shift of K45 towards the binding pocket
in the mutant structure.
Overall, similar interactions were reported between wildtype and mutant proteins
with different inhibitors, with the exception of GRL-5010,
which interacted with D30$^\prime$ in an altered way.
These results thus partially support the observations made in our study on the effects of the L76V mutation
on the structure and interactions within the HIV protease.

\section*{Conclusions}

In this work, we analysed a set of mutations in the HIV protease associated
with resistance towards PIs in complex with different clinically used inhibitors.
First, we analysed three mutations with resistance factors extracted from the literature,
where resistance factor measurements
for the same sequence and the same inhibitors were available from multiple experiments.
We showed that the effect of the mutation on the resistance factor,
both increasing resistance and sensitizing, was successfully reproduced using
alchemical free energy calculations of affinity of inhibitor binding.
Second, we modelled complexes for sequences based on our own clinical samples
containing the mutation L76V with four PIs.
Even though the sequences in question had a large number of background mutations,
we could in most cases reproduce the resistant and sensitizing effects of L76V.
These calculations gave us insight into whether change in
resistance is predominantly the result of change
in inhibitor binding affinity or a change in the catalytic activity of the protease,
for example for sequences which lacked the compensatory mutation M46I.
Further analysis of L76V in different sequence contexts revealed that the
effect of this mutation on direct protein residue-inhibitor interactions,
including that of D30, is generally line with the changes in the resistance.
Potentially causal to the observed changes
is the favourable effect of the mutation on the
hydrogen bond stability between residues D30 and K45 of the binding pocket.
Analysis of another distant site mutation, N88S,
also revealed changes in hydrogen bonding
of the mutated residue with D30 as well as with T31 and T74,
suggesting changes in hydrogen bonding network
of the protease as a major pathway for how mutations outside
of the binding pocket affect inhibitor binding.

\begin{backmatter}
\section*{Methods}
\subsection*{Computational Studies}

\paragraph{System Preparation}

Crystal structures of protease-inhibitor complexes were obtained from the Protein Data Bank \cite{Berman00}
(PDB IDs 1HPV (APV), 1HXB (SQV), 1K6C (IDV), 1MUI (LPV), 2BPX (IDV), 3EKV (APV), 3EL1 (ATV), 3PWR (SQV)).
Modeller \cite{Sali93} version 9.12 was used to introduce mutations targeted in this study
as well as the background mutations. The following background mutations were introduced in the studied protein from HIVdb dataset:
K7Q, R14K, R57G, T82V, and V84I in 1K6C;
K7Q, R14K, K41R, and V77I in 3EKV;
L10F and S37N in 1MUI;
L10F and S37N in 1HPV (84V);
S37N and A71V in 1HPV (46V);
L10F and S37N in 2BPX (84V);
S37N and A71V in 2BPX (46M).
For the phenotypic assay dataset the following background mutations were introduced:
K7Q, I13V, G16E, K20I, I33F, M36L, S37N, I62V, I63H, A67C, A71V, G73S, I84V, L90M, and 95C in 3PWR for the genotype FB15;
K7Q, I13V, G16E, K20I, I33F, M36L, K41R, I62V, P63H, V64I, A71V, G73S, I84V, and L90M in 3EL1 for the genotype FB15;
K7Q, L10V, I13V, G16E, K20R, I33L, E35D, M36I, S37N, M46I, I54V, Q58E, I62V, I63H, I64V, A67C, V82F, I84V, and A95C
in 3PWR for the genotype GH9;
L10I, I13V, K20M, E35D, M36I, S37N, R41K, M46I, I54V, Q58E, H69K, V82M, and L89I in 1MUI for the genotype RU1;
K7Q, L10V, I13V, R14K, K20M, E35D, M36I, M46I, I54V, Q58E, P63L, V64I, H69K, V82M, and L89I in 3EL1 for the genotype RU1;
L10I, I13V, K20M, E35D, M36I, S37N, R41K, M46I, I54V, Q58E, H69K, V82M, and L89I in 1MUI for the genotype RU1;
K7Q, L10I, I13V, R14K, L24I, L33F, K46L, I62V, A71V, T82A, V84I, and Q92K in 1K6C for the genotype iZ2;
L10I, I13V, L24I, L33F, S37N, M46L, I62V, L63P, A71V, and V82A in 1MUI for the genotype iZ2;
K7Q, L10I, I13V, R14K, L24I, L33F, K41R, M46L, I62V, V64I, A71V, V82A, and Q92K in 3EL1 for the genotype iZ2.

In the following, preparation for simulations of all the structures mentioned above is described
in both holo and apo states,
with the exception of structure 1HXB, 1MUI, and 2BPX for the HIVdb dataset,
for which simulations in apo state were not performed
(see Equilibrium MD Simulations and Free Energy Calculations).

Remaining set up of the system
in this study has been performed in the manner as described previously \cite{Bastys18}.
In short, the Gromacs simulation software package was used to set up (version 4.6.5), carry out,
and analyse the MD simulations (versions 5.0.2 and 5.1.2) \cite{Hess08, Abraham15}.
The $pK_a$ of residues was predicted using Propka \cite{Sondergaard11}
and protease was assigned monoprotonated state on either D25/D25$^\prime$,
here the prime refers to the second subunit of the protein.
The Amber99SB*-ILDN force field was used for parametrisation of the protease.
The Chemaxon Calculator\cite{Chemaxon} was used to determine inhibitor protonation,
while Gaussian09~\cite{g09} was used to optimize inhibitor geometry and calculate electrostatic potential.
Partial charges were assigned by performing restrained electrostatic potential (RESP) fit~\cite{Bayly93}.
The complex was solvated in TIP3P water molecules with 1.4~nm buffer in each dimension
with 0.15~mol/liter concentration of Cl and Na ions \cite{Joung08} to neutralize the system.

\paragraph{Equilibrium MD Simulations and Free Energy Calculations}

The equilibrium MD simulations and the free energy calculations
in this study have been performed in the manner as described previously \cite{Bastys18}.
In short, each system was subjected to steepest descent energy minimisation.
Before equilibrium simulation,
in order to avoid too close contacts between atoms,
simulated annealing in length of 1~ns was performed for the following complexes:
1MUI (76V for the genotype iZ2) in D25$^\prime$ protonation state,
1MUI (76L for the genotype iZ2) in D25 protonation state,
2BPX (46V and 84I variants) in D25 protonation state,
3EKV (88S variant) in D25 protonation state,
3EL1 (76L and 76V for the genotype iZ2) in D25$^\prime$ protonation state.
Ten replicas of 200~ns simulation for each complex were performed at 300~K.
For all of the analyses that followed, the first 20~ns of the simulations
were considered to be a part of the system equilibration process and thus discarded,
with the exception of free energy calculations, where first 10~ns were discarded.
The protocol for free energy calculations was adjusted from the non-equilibrium simulation approach
used in assessing changes in protein thermal stabilities and protein--protein interactions upon amino acid mutation~\cite{Gapsys16}.
For calculating the free energy change upon mutation of apo structures, $\DG_1$,
estimates for structure 1HPV were used in both $\DDG$ estimates for mutation M46I
and all four $\DDG$ estimates for mutation I84V,
as the background sequences for these mutations were the same.

From each of the equilibrium simulations described above,
trajectory frames were extracted equidistantly in time every 10~ns.
After generating hybrid structures for every snapshot using the pmx~\cite{Gapsys15} framework,
short 20~ps simulations were performed to equilibrate velocities,
after which alchemical transitions were carried out in 50~ps.
Identical parameters were used for equilibrium simulations,
equilibration, and alchemical transitions with soft-core potential for non-bonded interactions \cite{Gapsys12}.
The Crooks Fluctuation Theorem~\cite{Crooks99} was used to 
relate the obtained work distributions to the free energy values
by employing maximum likelihood estimator~\cite{Shirts03},
with the error estimates obtained by the bootstrap approach.

\paragraph{Partial Least-Squares Regression}

Partial least-squares regression was performed with the functional mode analysis tool \cite{Krivobokova12}
using the heavy atoms of protein as predictors.
Wildtype and mutant protein simulation trajectories were labelled using constants 0 and 1 as target values, respectively.

For each mutation and inhibitor combination cross validation (CV) was performed to verify the models.
During CV, all trajectories for wildtype and mutant complexes were concatenated,
superimposed to minimize the variance over the ensemble \cite{Gapsys13},
and divided into five equal parts.
For every iteration,
a model was trained on four parts of labelled input in equal proportions from wildtype and mutant simulations,
after which it was used to make prediction for the last part.
The Pearson correlation between the actual and predicted labels was used to evaluate the quality of the model.
Number of components $i = 1, \dots, 25$ was tested in each iteration.
For the final model, the number of components was chosen from the correlation curve in CV
such that choosing a higher number of components only marginally improves the performance of the model.

\paragraph{Estimation of Resistance Factor Change from Free Energy of Inhibitor Binding Change}

Cheng-Prusoff equation \cite{Cheng73} relates inhibitor's $\IC$ and binding affinity $K_i$,
which in turn can be estimated from inhibitor binding free energy $\DG$:
\begin{equation}
    \IC = K_i\bigg(1+\frac{[S]}{K_m}\bigg) = e^\frac{\DG}{k_B T}\bigg(1+\frac{[S]}{K_m}\bigg)\,,
\end{equation}
where [S] is a fixed substrate concentration,
$K_m$ is the concentration of the substrate at which the enzyme is at its half-maximal activity,
$k_B$ is the Boltzmann constant, and T is the temperature.
Thus, given two RF values for two proteases with sequences $A$ and $B$, $\RF_A$ and $\RF_B$, their ratio can be related to $\DDG$:

\begin{equation}
\RF_R = \frac{\RF_A}{\RF_B} = \frac{e^\frac{\DG_A}{k_BT}}{e^\frac{\DG_B}{k_BT}}\Bigg(\frac{1+\frac{[S]}{K_m^A}}{1+\frac{[S]}{K_m^B}}\Bigg) =
e^{\frac{\Delta\Delta G}{k_BT}}\Bigg(\frac{1+\frac{[S]}{K_m^A}}{1+\frac{[S]}{K_m^B}}\Bigg) = e^{\frac{\Delta\Delta G}{k_BT}} C\,.
\end{equation}

We are interested in obtaining a distribution of the $\RF_R$ values after calculating the
double free energy differences $\DDG$:
\begin{equation}
    p(\RF_R |\DDG, C) \propto p(\DDG, C|\RF_R)p(\RF_R)
\end{equation}

When there are multiple $\RF_R$ measurements and $\DDG$ calculations available for the wildtype
and mutant protein complexed with different ligands, $C$ can
be expressed as a function of the available values $C = C_i(\DDG_i , \RF_R^i), i = 1, \dots, n$.
This gives the final posterior distribution:
\begin{equation}
    p(\RF_R |\DDG, \DDG_i , \RF_R^i) \propto p(\DDG, \DDG_i , \RF_R^i | \RF_R) p(\RF_R),
    \label{eq:RF_R_calc}
\end{equation}
\noindent where $C_i = \RF_R^i e^\frac{-\DDG_i}{k_B T}$.
    The $\DDG$ values are sampled from a Gaussian distribution
with the mean and standard deviation corresponding to the calculated double free energy
difference and estimated error, respectively.

\subsection*{Phenotypic Assay for Resistance Factor Value Estimation}

The experimental data on L76V resistance is based on
samples of patients who underwent multiple therapy failures with different PIs.
Those variants were observed in the diagnostic procedure, sequenced, and subsequently tested in a phenotypic assay
as described by Walter et al.~\cite{Walter99}.
The tests were carried out after the patient's variant was cloned into a recombinant derivate of the HIV NL4-3, called pNL4-3-Delta-PRT5.
The L76V mutation was reverted to wildtype by site directed mutagenesis.
This allowed to determine the effect of the genetic background upon the L76V.
These variants were analysed in cell culture experiments where they were exposed to different PIs in different concentrations
to estimate their RF values (Table \ref{tab:pr:clinical_RF}).
Based on these variants, the clones were specifically modified by site-directed mutagenesis
so that different variants of L76V could be tested in different genetic backgrounds. 
For simplicity, regardless of the residue at position 76 of protease as present
in the original clinical samples, in the context of this paper L76 is referred to as wildtype residue and
V76 as the mutant residue as per reference sequence HXB2.


\printacronyms[include-classes=abbrev,name=Abbreviations]

%\section*{Declarations}
\section*{Declarations}

\subsection*{Ethics approval and consent to participate}

The study was approved by the University Clinics Cologne (16-460).

\subsection*{Consent for publication}

All study participants consent to publication.

\subsection*{Availability of data and materials}

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

\subsection*{Competing interests}

The authors declare no competing interests.

\subsection*{Funding}

OVK is grateful to the Klaus Faber Foundation for their financial support.
VG and BLdG were supported by the BioExcel CoE (www.bioexcel.eu),
a project funded by the European Union contracts H2020-INFRAEDI-02-2018-823830, H2020-EINFRA-2015-1-675728.

\subsection*{Author's contributions}

All authors wrote and approved the manuscript.

\subsection*{Acknowledgements}

We are grateful to Alejandro Pironti for providing resistance factor data from HIVdb and 
to Mazen Ahmad for discussions and technical advice on the simulations.
The simulations presented in this work were in part performed on HPC systems of the Max Planck Computing and Data Facility.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  The Bibliography                       %%
%%                                                         %%
%%  Bmc_mathpys.bst  will be used to                       %%
%%  create a .BBL file for submission.                     %%
%%  After submission of the .TEX file,                     %%
%%  you will be prompted to submit your .BBL file.         %%
%%                                                         %%
%%                                                         %%
%%  Note that the displayed Bibliography will not          %%
%%  necessarily be rendered by Latex exactly as specified  %%
%%  in the online Instructions for Authors.                %%
%%                                                         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% if your bibliography is in bibtex format, use those commands:
\bibliographystyle{bmc-mathphys} % Style BST file (bmc-mathphys, vancouver, spbasic).
%\bibliography{manuscript}
\bibliography{bmc_article}      % Bibliography file (usually '*.bib' )
% for author-year bibliography (bmc-mathphys or spbasic)
% a) write to bib file (bmc-mathphys only)
% @settings{label, options="nameyear"}
% b) uncomment next line
%\nocite{label}

% or include bibliography directly:
% \begin{thebibliography}
% \bibitem{b1}
% \end{thebibliography}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Figures                       %%
%%                               %%
%% NB: this is for captions and  %%
%% Titles. All graphics must be  %%
%% submitted separately and NOT  %%
%% included in the Tex document  %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
%% Do not use \listoffigures as most will included as separate files

\section*{Figures}
\begin{figure}[!ht]
    \caption{HIV protease structure. Flap region in cyan, 80s loop in brown, active-site proximate loop in olive colours.
	Mutations analysed in this study (red), catalytic site residue (blue) and bound inhibitor (magenta) are shown in sticks model.}
	\label{fig:PR_structure}
\end{figure}


\begin{figure}[!ht]
    \caption{Predicted and experimental RF measurements.
	Each symbol corresponds to a unique sequence background
    and colour corresponds to inhibitor.
\emph{Nota bene:} in case of APV, $\RF^{exp}_R$ measurements are for its prodrug FPV.}
    \label{fig:RF}
\end{figure}

\begin{figure}[!ht]
    \caption{Predicted and experimental RF measurements.
	Each symbol corresponds to a unique sequence background
    and colours correspond to inhibitor.}
    \label{fig:RF_L76V}
\end{figure}

\begin{figure}[!ht]
    \caption{Interpolation between the extremes of the FMA models for the protease (genotype RU1) in complex with LPV.
Blue-to-magenta bands correspond to the interpolation along the mode as represented
as cartoon for backbone and as sticks for residues 30, 45, and 58,
with blue corresponding to L76 state and magenta to V76 state.
Green dashed line represents a hydrogen bond between residues D30 and K45.
Mutated residue 76, here semi-transparent in yellow, as well as hydrogen atoms, here in gray,
were not part of the FMA models and are here for representational purposes only.}
	\label{fig:pr:hbond}
\end{figure}

\begin{figure}[!ht]
    \caption{Energy differences of non-bonded interaction between protein and
	inhibitor in wildtype and mutant complexes.
	Residues, for which the difference ($E_{MUT}-E_{WT}$) between the wildtype
	and the mutant complexes is higher than the propagated error (SE)
	and its absolute value higher than 0.1~kcal/mol, are represented as a colored circle,
	where the color represents relative interaction energy and the size
	of the circle relates inversely to the standard error of the estimate.
    Residues' 30 and 45 interactions are highlighted in green box.}
    \label{fig:energy}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Tables                        %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Use of \listoftables is discouraged.
%%

\section*{Tables}
\begin{table}[!ht]
    \caption{RF values for different mutant and corresponding wildtype sequences from HIVdb \cite{HIVdb}.
        First column indicates the mutation analysed,
	while the second column indicates background polymorphisms that are present in both wildtype
	and mutant sequences compared to the reference HIV sequence HXB2.
    Multiple RF measurements for the same protein are separated by comma.
	* indicates that the R57G background mutation was found in sequences where RF
	for IDV was measured and L77I in sequences where RF for FPV was measured.
	M46I wildtype and mutant measurements reported from different studies,
but performed using the same susceptibility test method (Phenosense, Monogram,
San Francisco, USA).}
\newcolumntype{B}{>{\bf}c}
        \resizebox{\textwidth}{!}{
        \begin{tabular}{Bccccccc}
	\toprule
    \multicolumn{2}{c}{\multirow{2}{*}{\bf{Mutation}}}	& \bf{Background}	& \multicolumn{4}{c}{\bf{RF}}\\
\cmidrule(lr){4-7}
		&   & \bf{polymorphisms}	& \bf{IDV}	& \bf{SQV} & \bf{LPV} & \bf{FPV}\\
	\midrule
	\multirow{2}{*}{M46I}	& wildtype~\cite{Colonno04}	& \multirow{2}{*}{V3I, S37N, A71V}	& 0.6, 1.0, 1.2 & --	& --	& 0.3, 0.6, 0.7 \\
				& mutant \cite{Petropoulos00}	&	& 4.4	& --	& --	& 2.2 \\
	\midrule
	\multirow{2}{*}{I84V}	& wildtype~\cite{Prado02}	& \multirow{2}{*}{V3I, L10F, S37N}	& 1.5 & 1.3 	& 1.6	& 1.8 \\
				& mutant\cite{Prado02}	&	& 2.1, 3.2	& 2.7, 3.7	& 6.2, 7.7	& 4.6, 8.4 \\
	\midrule
	\multirow{2}{*}{N88S}	& wildtype~\cite{Ziermann00}	& \multirow{2}{*}{V3I, S37N, L63P, R57G/L77I*}	& 1.1	& --	& --	& 1.0 \\
				& mutant~\cite{Ziermann00}	&	& 2.6	& --	& --	& 0.1\\
	\bottomrule
    \end{tabular}
    }
    \label{tab:HIVdb_RF}
\end{table}

\begin{table}[!ht]
    \begin{center}
        \caption{Change of the binding free energy ($\DDG$) of inhibitors upon mutation.
    All values in kcal/mol, and $\pm$ shows bootstrap error estimate.
$\RF^{exp}_R$ indicates value ranges calculated using $\RF_R = \RF_{mutant}/\RF_{wildtype}$
from previously reported experimental RF measurements,
where * stands for measurements for FPV, the prodrug of APV.
$\RF^{calc}_R$ indicates average calculated $\RF_R$ value
from the distribution described in Equation \ref{eq:RF_R_calc}
with subscript and superscript corresponding to lower and upper bound of 95\% credible interval.}
    \setlength{\tabcolsep}{9pt}
    \setlength\extrarowheight{2pt}
    \begin{tabular}{ccr@{~$\pm$~}lr@{--}lc}
	\toprule
    Inhibitor & Mutation & \multicolumn{2}{c}{$\DDG$} & \multicolumn{2}{c}{$\RF^{exp}_R$} & $\RF^{calc}_R$\\
	\midrule
    APV & 	M46I	& $-0.35$	& $0.4$	    & 3.14 & 7.33*   & 13.13$^{77.54}_{0.3}$ \\
    IDV &	M46I	& $-0.34$	& $0.72$ 	& 3.67 & 7.33& 12.69$^{75.84}_{0.27}$ \\
    APV & 	I84V	& $2.06$	& $0.47$ 	& 2.56 & 4.67*   & 33.75$^{169.66}_{1.25}$ \\
    IDV & 	I84V	& $1.13$	& $0.57$ 	& 1.4  & 2.13 & 7.23$^{42.3}_{0.11}$  \\
    LPV & 	I84V	& $1.25$	& $0.39$ 	& 3.88 & 4.81    & 6.09$^{33.62}_{0.17}$ \\
    SQV & 	I84V	& $0.56$	& $0.41$ 	& 2.08 & 2.85    & 1.33$^{6.73}_{0.05}$ \\
    APV & 	N88S	& $-0.97$	& $0.7$ 	& \multicolumn{2}{c}{0.1*}  & 0.3$^{2.14}_{0.001}$ \\
    IDV & 	N88S	& $1.41$	& $0.95$ 	& \multicolumn{2}{c}{2.27}  & 41.12$^{281.92}_{0.11}$ \\
	\bottomrule
    \end{tabular}
    \label{tab:RF_ddG}
    \end{center}
\end{table}

\begin{table}[!ht]
\newcolumntype{B}{>{\bf}c}
    \begin{center}
    \caption
    {Protease RF values for L76V mutation.
	Multiple RF measurements for the same protein separated by comma.
	For each genotype the first row represents the wildtype position as in the original sample
	and the second row represents the mutation introduced at position 76.
    }
    \begin{tabular}{BBccccc}
	\toprule
    \multicolumn{2}{c}{\multirow{2}{*}{\bf{Genotype}}}	& \multicolumn{4}{c}{\bf{RF}}\\
\cmidrule(lr){3-6}
	    &	& \bfseries{ATV}	& \bfseries{SQV}	& \bfseries{IDV}	& \bfseries{LPV}\\
	\midrule
	\multirow{2}{*}{FB15}	& L76		& 63			& 74			& --			& -- 		\\
				& V76		& 0.9, 2.3, 4.5		& 3.6, 4.6, 5.8		& --			& -- 		\\
	\midrule
	\multirow{2}{*}{GH9}	& L76		& 90			& 18.6			& --			& -- 		\\
				& V76		& 1.2, 1.9, 2.4, 3.2, 3.6	& 1.2, 1.3, 1.5	& --			& -- 		\\
	\midrule
	\multirow{2}{*}{RU1}	& V76		& 2.7			& --			& --			& 157		\\
                & L76		& 8.4, 9, 10		& --			& --			& 27, 46, 47	\\
	\midrule
	\multirow{2}{*}{iZ2}	& V76		& 4.1			& --			& 59			& 71		\\
                & L76		& 7.8, 12, 30		& --			& 2.2, 7.9, 10		& 5.5, 11, 12	\\
	\bottomrule
    \end{tabular}
    \label{tab:pr:clinical_RF}
    \end{center}
\end{table}

\begin{table}[!ht]%\scriptsize
\setlength{\tabcolsep}{9pt}
    \setlength\extrarowheight{2pt}
\caption{Change of the binding free energy ($\DDG$) of inhibitors upon mutation L76V.
	All values in kcal/mol, and $\pm$ shows bootstrap error estimate.
$\RF^{exp}_R$ indicates value ranges calculated using $\RF_R = \RF_{mutant}/\RF_{wildtype}$
from experimental RF measurements, $\RF^{calc}_R$ indicates average calculated $\RF_R$ value
from the distribution described in Equation \ref{eq:RF_R_calc}
with subscript and superscript corresponding to lower and upper bound of 95\% credible interval.}
\begin{tabular}{*{1}{c}c*{1}{r@{~$\pm$~}l}r@{--}lc}
	\toprule
    Inhibitor & Genotype & \multicolumn{2}{c}{$\DDG$}  & \multicolumn{2}{c}{$\RF^{exp}_R$} & $\RF^{calc}_R$\\
	\midrule
    ATV & 	FB15	& $0.78$	& $0.64$	& 0.01 & 0.07    & 0.23$^{1.41}_{0.004}$ \\
    SQV &	FB15	& $0.69$	& $0.62$	& 0.04 & 0.08    & 0.11$^{0.71}_{0.001}$ \\
    ATV & 	GH9	& $-0.29$	& $0.49$	& 0.01     & 0.04    & 0.23$^{1.12}_{0.01}$ \\
    SQV & 	GH9	& $-0.58$	& $0.47$	& 0.06     & 0.08    & 0.03$^{0.17}_{0.001}$ \\
    ATV & 	RU1	& $0.74$	& $0.73$	& 0.27     & 0.32    & 3.52$^{22.7}_{0.04}$ \\
    LPV & 	RU1 	& $1.52$	& $0.6$	& 3        & 5  & 3.94$^{25.08}_{0.05}$ \\
    ATV & 	iZ2	& $1.11$	& $0.65$	& 0.14     & 0.53    & 35.61$^{243.1}_{0.17}$ \\
    IDV & 	iZ2	& $1.76$	& $0.66$	& 5.9      & 26.82    & 88.25$^{649.19}_{0.06}$ \\
    LPV & 	iZ2	& $0.99$	& $0.88$	& 5.91     & 12.9    & 10.95$^{76.68}_{0.01}$ \\
	\bottomrule
    \end{tabular}
    \label{tab:ddG_L76V}
\end{table}

\begin{table}[!ht]%\scriptsize
\newcolumntype{L}{>{$}l<{$}}
\newcolumntype{R}{>{$}r<{$}}
\centering
    \caption
    {Average number of hydrogen bonds between residues D30 and K45 for protease wildtype and mutant complexes.\ $\pm$
	indicates standard error of bond frequency across independent simulations.}
    \resizebox{\textwidth}{!}{
    \begin{tabular}{ @{}c@{~} @{~}c R@{~$\pm$~}L c R@{~$\pm$~}L R@{~$\pm$~}L@{} }
	\toprule
	\multirow{2}{*}{\textbf{Drug}} & \multirow{2}{*}{\textbf{Genotype}} & \multicolumn{3}{c}{\textbf{D30--K45}} & \multicolumn{4}{c}{\textbf{D30$^\prime$--K45$^\prime$}} \\
	\cmidrule(lr){3-5} \cmidrule(lr){6-9} % with \usepackage{booktabs}
	    &	& \multicolumn{2}{l}{\textcolor{white}{0000}\textbf{L76}}	& \textbf{V76}	& \multicolumn{2}{l}{\textcolor{white}{000}\textbf{L76}} & \multicolumn{2}{l}{\textcolor{white}{000}\textbf{V76}}	\\
	\midrule
	ATV & 	FB15	& 0.068 & 0.003 & $0.52\pm0.04$	& 0.07 & 0.002	& 0.58 & 0.02	\\
	SQV &	FB15	& 0.12 & 0.01	& $0.59\pm0.02$	& 0.11 & 0.005	& 0.62 & 0.004	\\
	ATV & 	GH9	& 0.07 & 0.01	& $0.49\pm0.12$	& 0.06 & 0.002	& 0.66 & 0.12	\\
	SQV & 	GH9	& 0.07 & 0.004	& $0.46\pm0.04$	& 0.01 & 4\times10^{-5}	& 0.19 & 0.04	\\
	ATV & 	RU1	& 0.05 & 0.001	& $0.43\pm0.05$	& 0.04 & 1\times10^{-4}	& 0.61 & 0.09	\\
	LPV & 	RU1 	& 0.01 & 3\times10^{-5}	& $0.16\pm0.03$	& 0.16 & 0.006	& 0.56 & 0.12	\\
	ATV & 	iZ2	& 0.16 & 0.007	& $0.57\pm0.07$	& 0.08 & 0.002	& 0.52 & 0.04	\\
	IDV & 	iZ2	& 0.1 & 0.01	& $0.39\pm0.03$	& 0.06 & 0.002	& 0.51 & 0.07	\\
	LPV & 	iZ2	& 0.04 & 5\times10^{-4}	& $0.19\pm0.02$	& 0.24 & 0.01 & 0.63 & 0.02	\\
	\bottomrule
	\end{tabular}
  }
    \label{tab:hbond_L76V}
\end{table}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Additional Files              %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section*{Additional Files}
\subsection*{Additional file 1 --- Tab.~S1. The names of the isolates, whose RF data was used in this study, as reported in HIVdb, and the reference of the study where RF measurements were performed. Tab.~S2. Inhibitor binding free energy change upon switching the proton from the reference protonated active site residue to the active site residue on the opposite subunit for wildtype and mutant proteins. $\pm$ shows bootstrap error estimate, all values in kcal/mol. Tab.~S3. Average hydrogen bonds number between residues D30, T31, and T74 with N88 and S88 for wildtype	and mutant complexes, respectively. Columns 3 and 4 of the table corresponds to hydrogen bonds within monomer A of protease and columns 5 and 6 of the table corresponds to hydrogen bonds within monomer B of protease (residues marked with prime symbol).~$\pm$ indicates standard error	of bond frequency across independent simulations. Fig.~S3. Inhibitor binding free energy change upon switching the proton from the reference protonated active site residue to the active site residue on the opposite subunit for wildtype and mutant proteins. $\pm$ shows bootstrap error estimate, all values in kcal/mol. Fig.~S3. Interpolation between the extremes of the FMA models for the corresponding complexes. Blue-to-magenta bands correspond to the interpolation along the mode as represented as cartoon for backbone and as sticks for residues 30, 45, and 58, with blue corresponding to L76 state and magenta to V76 state. Mutated residue 76 is not part of the model and is represented here as gray dash. Tab.~S4. Inhibitor binding free energy change upon switching the proton from the reference protonated active site residue to the active site residue on the opposite subunit for wildtype and mutant proteins. $\pm$ shows bootstrap error estimate, all values in kcal/mol. Fig.~S4. Energy differences of non-bonded interactions between protein and inhibitor in wildtype and mutant complexes. Only residues, for which the difference between the wildtype and the mutant complexes is higher than the propagated error and its absolute value higher than 0.1~kcal/mol are shown.}

\end{backmatter}
\end{document}
