Aline Taoum
email: aline.taoum@utt.fr

Farah Mourad-Chehade
email: farah.chehade@utt.fr

Hassan Amoud
email: hassan.amoud@ul.edu

Evidence-based Model for Real-time Surveillance of ARDS

Keywords: Acute respiratory distress syndrome, Evidence-based theory, Linear and nonlinear parameters, Real-time surveillance

Real-time health surveillance becomes important and necessary with the increase of the elderly population to preserve their quality of life. Real-time models aim to provide alerts before the severe illness occurs. Acute respiratory distress syndrome is a crucial disease of the respiratory system that threats the health of the elderly. This paper proposes a real-time model for the surveillance of ARDS based on belief functions theory. Non-invasive physiological signals are considered such as heart rate, respiratory rate, oxygen saturation and mean airway blood pressure. Different linear and nonlinear parameters are extracted from these signals; then a parameters selection procedure is performed to reduce their dimensionality. Afterwards, classifiers are constructed using parameters distributions defined in the evidence framework. Real-time prediction is then performed by combining all classifiers decisions. As results, high performances are obtained over the testing sets with performances of 77% and 71% for sensitivity and specificity, respectively.

Introduction

The elderly population is expected to increase from 962 million in 2017 to more than 2 billion by 2050 [START_REF]World Population Prospects: The 2017 Revision, Key Findings and Advance Tables[END_REF]. Older people are prone to many pathologies that reduce their autonomy and may threaten their lives leading to worldwide economic burden. The most known elderly-related pathologies are cardiovascular diseases, chronic respiratory 5 diseases, neurological and mental disorders [START_REF] Prince | The burden of disease in older people and implications for health policy and practice[END_REF]. A real-time surveillance for elderly people helps to detect pathologies in advance. Thus it prevents its severity, shortens lengths of hospital stay and reduces mortality [START_REF] Grace | Perceptions of seniors with heart failure regarding autonomous zero-effort monitoring of physiological parameters in the smart-home environment[END_REF]. Real-time surveillance consists of continuous analysis of health-related data collected at a fixed time interval to reveal upcoming illness. It has been the subject of debate for many researchers [START_REF] Lorenz | Mobile health monitoring for the elderly: Designing for diversity[END_REF][START_REF] Sufi | Novel methods of faster cardiovascular diagnosis in wireless telecardiology[END_REF][START_REF] Forkan | A clinical decision-making mechanism for context-aware and patient-specific remote monitoring systems using the correlations of multiple vital signs[END_REF]. In this paper, we present a real-time surveillance system for Acute Respiratory Distress Syndrome (ARDS).

ARDS is a severe disability of the respiratory system characterized by expanded pulmonary inflammation [START_REF] Ware | The acute respiratory distress syndrome[END_REF][START_REF] Steinberg | National Heart, Lung, and Blood Institute, Acute Respiratory Distress Syndrome (ARDS) Clinical Trials Network, Efficacy and Safety of Corticosteroids for Persistent Acute Respiratory Distress Syndrome[END_REF]. It results in insufficiency of gas exchange with the blood, associated with high short-term mortality [START_REF] Ashbaugh | Acute Respiratory Distress in Adults[END_REF]. Over years, several definitions had been published until the Berlin Definition was proposed in 2012 [START_REF] Ashbaugh | Acute Respiratory Distress in Adults[END_REF][START_REF] Bernard | The American-European Consensus Conference on ARDS. Definitions, mechanisms, relevant outcomes, and clinical trial coordination[END_REF][START_REF] Ranieri | Acute Respiratory Distress Syndrome: The Berlin Definition[END_REF]. Accordingly, ARDS is a hypoxemia that develops within one week of new respiratory inflammation symptoms or a known clinical insult. It is diagnosed using the ratio of the partial pressure of oxygen in the arterial blood (P aO 2 ) over the fraction of oxygen in the inspired air (F iO 2 ). Moreover, ARDS was classified as mild (200 mmHg < P aO 2 /F iO 2 ≤ 300 mmHg), moderate (100 mmHg < P aO 2 /F iO 2 ≤ 200 mmHg), or severe ARDS (P aO 2 /F iO 2 ≤ 100 mmHg). All having a minimum positive endexpiratory pressure (PEEP) of 5 cm H 2 O.

The risk of ARDS increases for patients who have undergone cardiac surgery or were hospitalized for severe health problems. Depending on the case, some patients were transferred from intensive care units (ICU) to wards or were completely discharged from hospitals. Therefore, the monitoring and early-identification of those patients is crucial in order to implement preventative measures. Generally, studies conducted on ARDS have used clinical variables to characterize it or to predict its onset. Some studies have investigated the risk factors associated with ARDS [START_REF] Yu | Role of diabetes in the development of acute respiratory distress syndrome[END_REF][START_REF] Singh | Incidence and outcome of acute lung injury and acute respiratory distress syndrome in the surgical intensive care unit[END_REF], others have developed predictive algorithms for ARDS using clinical data of ventilated ICU patients, as in [START_REF] Ennett | Predicting respiratory instability in the ICU[END_REF]. Other researchers have worked on mortality prediction related to ARDS [START_REF] Navarrete-Navarro | Acute respiratory distress syndrome in trauma patients: Icu mortality and prediction factors[END_REF][START_REF] Ichikado | Prediction of prognosis for acute respiratory distress syndrome with thin-section ct: Validation in 44 cases 1[END_REF] or respiratory complications [START_REF] Ravishankar | An early respiratory distress detection method with Markov models[END_REF][START_REF] Velardo | Automatic generation of personalised alert thresholds for patients with COPD[END_REF].

Furthermore, data mining and machine learning techniques are widely adopted in biomedical applications to analyze, classify and predict medical data [START_REF] Yoo | Data mining in healthcare and biomedicine: A survey of the literature[END_REF]. Classical machine learning techniques such as k-nearest neighbors [START_REF] Ler | Combining classical hrv indices with wavelet entropy measures improves to performance in diagnosing congestive heart failure[END_REF], naive Bayes [START_REF] Seera | Classification of electrocardiogram and auscultatory blood pressure signals using machine learning models[END_REF] and support vector machine [START_REF] Clifton | Identification of patient deterioration in vital-sign data using one-class support vector machines[END_REF] are applied on medical data to diagnose or predict a clinical abnormality in physiological signals. However, there is a lack in the literature of prediction models and real-time surveillance for ARDS.

This paper proposes a novel algorithm for predicting the occurrence of ARDS.

A key point in this surveillance is considering multiple vital signs that are easy to collect [START_REF] Jassas | A smart system connecting e-health sensors and the cloud[END_REF][START_REF] Fang | Multi-parameter health monitoring watch[END_REF]. Identifying ARDS can be done using just one signal [START_REF] Klingeberg | Mobile wearable device for long term monitoring of vital signs[END_REF], but most of the serious diseases occur with concurrent irregularities in multiple vital signs [START_REF] Forkan | A clinical decision-making mechanism for context-aware and patient-specific remote monitoring systems using the correlations of multiple vital signs[END_REF].

Hence, in this study, physiological signals are analyzed continuously to identify the current health state of a patient. Multiple vital signs are used, that are the heart rate, the respiratory rate, the oxygen saturation and the mean airway blood pressure. Instead of using the whole signals, linear and non-linear parameters are extracted in real-time and the most informative ones are selected. After that, a prediction model based on evidence theory is defined in order to combine information of the extracted parameters and generate a decision regarding the health state of the patient in study. This model is constructed using the multi-parameter intelligent monitoring of intensive care II database. The outline of the paper is as follows. In the next section, the patient extraction procedure is presented, and then the evidence-based method is introduced. Thereafter, the results are illustrated and interpreted. Finally, a conclusion summarizes the paper and proposes some perspectives.

Materials

Subjects are extracted from the multi-parameter intelligent monitoring of intensive care II (MIMIC II) database. It contains two types of data: monitor waveforms and clinical data, collected over a seven years period. A detailed explanation of the database can be found in [START_REF] Goldberger | PhysioToolkit, and PhysioNet : Components of a New Research Resource for Complex Physiologic Signals[END_REF][START_REF] Saeed | Multiparameter Intelligent Monitoring in Intensive Care II (MIMIC-II): A public-access intensive care unit database[END_REF]. In this study, the clinical database was considered first to select subjects, using the values of their P aO 2 /F iO 2 ratios. Hence, an ARDS group is formed using the definition of both moderate and severe ARDS with P aO 2 /F iO 2 ≤200 mmHg in addition to PEEP levels higher than 5 cmH 2 O; while a non-ARDS group included subjects having mild and non ARDS patients with PEEP<5 cm H 2 O. After that, the waveform database is matched to the selected subjects of the two groups, and their time series are extracted; that are the heart rate (HR), the respiratory rate (RR), the oxygen saturation (SpO 2 ) and the mean airway blood pressure (MABP). These time series have a frequency of 1 sample/min. This leads to 140 ARDS subjects and 135 non-ARDS subjects.

Among the ARDS subjects, there are 50 subjects that started their records after ARDS diagnosis. These subjects are irrelevant for this study and thus were eliminated. Other ARDS subjects started their records before the onset of ARDS but stopped before ARDS diagnosis; these subjects, of number equal to 18, are also eliminated. Therefore, only the 72 subjects having started records before the ARDS onset and last at least until ARDS occurrence are considered. After that, a preprocessing step was required to smooth the signals and remove the noise. Some subjects were removed due to noise or too short signals. This leads finally to 50 ARDS subjects. The non-ARDS group is then reduced to have equal number of subjects, leading to 50 non-ARDS subjects as well. Figure 1 illustrates an example of the four extracted time series for an ARDS subjects. Within the database of 100 subjects, a k-fold cross-validation technique is considered to set training databases and testing databases [START_REF] Shao | Linear model selection by cross-validation[END_REF]. In the following, the terms training set and testing set would denote any training database obtained at any step of the cross-validation process with its corresponding testing database respectively.

Evidence-based surveillance approach

The method described in this paper is a real-time method that monitors the ongoing health state of patients using their vital signs. Let x s,i (ℓ), i = 1, . . . , 4, denote the four time series considered in this study for a subject s, recorded at time ℓ, and let t denote the real-time. Moreover, let x s,i (a : b) = (x s,i (a), ..., x s,i (b)) be a segment of signal i for subject s going from time a to time b.

Description of the approach

The aim of this work is to set a model Ψ from a training set that distinguishes in real-time between subjects who are going to develop ARDS labeled " + 1", and subjects who will not, labeled " -1". The model Ψ takes as input the collected signals x s,i (1 : t), going from the beginning of records until the real-time t, and yields as 95 output a decision whether the subject in study is going to develop ARDS. Instead of using the whole signals, the proposed approach considers a fixed window of length τ taken at the end of the signals, leading to four segments x s,i (tτ + 1 : t) to be analyzed at each time t. The reason of analyzing the data of a window of length τ is to have fixed length signals in computations. Moreover, an abnormality is better and faster detected within a short segment, compared to the high length signal. The value of τ could be taken equal to 6, 12 or 24 hours. In this paper, it is obtained by an optimization over the training sets, as shown in the following. Then, a vector of parameters is extracted from these segments at real-time t. The beginning of the signals is yet needed to normalize these parameters. This vector of parameters is then used with the belief functions theory to assign evidence to each health state according to each parameter. The parameters are considered as sources of information, providing an evidence regarding each state of the patient at the given time. Thus, masses are assigned to each state " + 1" and " -1" according to each parameter. The masses are obtained by means of mass functions defined using the parameters extracted from the training set at an earlier offline phase. The masses are then combined to assign a final mass to each state. A higher mass assigned to the state " + 1" at time t means that the subject is predicted to develop ARDS in the future, whereas a higher mass to " -1" means that the ARDS is not yet predicted at the current time. Therefore, for a new test subject given it signals in real-time, a decision will be made on its future health state at the earliest possible time.

Parameters Extraction

The first step of the prediction model is the extraction of parameters. Having a segment x s,i (tτ + 1 : t) from each signal at the real-time t, several linear and nonlinear parameters are extracted as follows. In the following, t -would denote the lower bound of the considered time interval, that is, t -= tτ + 1. For linear parameters, we consider:

• The mean defined as

µ s,i (t) = 1 τ ℓ=t ℓ=t -
x s,i (ℓ);

• The standard deviation

σ s,i (t) = 1 τ ℓ=t ℓ=t - (x s,i (ℓ) -µ s,i (t)) 2 ;
The variance is equal to the total power of spectral analysis. Thus, it reflects all the cyclic components responsible for the variability in the data [START_REF] Camm | Heart rate variability: standards of measurement, physiological interpretation and clinical use. task force of the european society of cardiology and the north american society of pacing and electrophysiology[END_REF].

• The skewness known as the third moment. It measures the symmetry of the data and it is computed using

Sk s,i (t) = ℓ=t ℓ=t - (x s,i (ℓ) -µ s,i (t)) 3 τ × σ 3 s,i (t)
;

It has been proposed that asymmetric tails to the left shown by positive skewness or to the right with negative skewness reflect, respectively, the acceleration or deceleration capacity of the time series data, e.g. for heart rate signal, it is an approximate distinction of vagal and sympathetic effects on the cardiac modulations [START_REF] Bauer | Deceleration capacity of heart rate as a predictor of mortality after myocardial infarction: cohort study[END_REF].

• The kurtosis that describes the shape of the probability distribution and is calculated using the fourth moment of the data,

Kt s,i (t) = ℓ=t ℓ=t - (x s,i (ℓ) -µ s,i (t)) 4 τ × σ 4 s,i (t)
.

The kurtosis measures the concentration of the data around the mean and reflects the rigidity of the signal [START_REF] Rivera | Heart rate and systolic blood pressure variability in the time domain in patients with recent and longstanding diabetes mellitus[END_REF].

On the other hand, non linear parameters are also considered, such as:

• The sample entropy (SampEn), that illustrates the amount of complexity in data. It is easily applied to any type of time series, including physiological data such as heart rate variability and blood pressure variability [START_REF] Moorman | Cardiovascular oscillations at the bedside: early diagnosis of neonatal sepsis using heart rate characteristics monitoring[END_REF][START_REF] Jalali | Discovering hidden relationships in physiological signals for prediction of periventricular leukomalacia[END_REF]. It is computed using the negative natural logarithm of the conditional probability that a segment of length N having repeated itself for u samples within a tolerance r will also repeat itself for u + 1 samples [START_REF] Lake | Sample entropy analysis of neonatal heart rate variability[END_REF]. Hence, for the segment x s,i of length N = τ , SampEn is computed as follows,

SampEn s,i (u, r, τ, t) = -ln h s,i (u + 1, r) h s,i (u, r) ,
where u is the length of sequences to be compared, r the tolerance for accepting and h(u, r) is the number of pairs of patterns sequences of x s,i of length u whose distance is less than r.

• The detrended fluctuation analysis (DFA), that measures the scaling behavior of a segment [START_REF] Francis | Physiological basis of fractal complexity properties of heart rate variability[END_REF]. It is defined by the short and long range scaling exponents DFA 1 and DFA 2 . To perform detrended fluctuation analysis for a segment x s,i (tτ + 1 : t), the integrated series y(u), u = 1, . . . , τ , is first computed

y(u) = ℓ=u ℓ=1 [x s,i (ℓ + t -τ ) -µ s,i (t)].
Then, y(u) is divided into time segments of size n, and a local trend y n (u) is obtained by a least-squares line fit and subtracted from y(u). The fluctuation f (n) for segments of length n is calculated as follows:

f (n) = 1 τ u=τ u=1 (y(u) -y n (u)) 2 .
The fluctuation f (n) is computed for different segment sizes n, then a graph of log f (n) against log n is constructed. Finally, the scaling exponents DFA 1 and DFA 2 define, as their names denote, the slope of the relation between log f (n) and log n for ranges of low and high values of n.

DFA has been used to assess the self-similar correlations of the time series [START_REF] Peng | Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series[END_REF].

It has also been described as a measure of "roughness" in the signals, with higher scaling exponents representing a smoother time series [START_REF] Millar | Isometric handgrip training lowers blood pressure and increases heart rate complexity in medicated hypertensive patients[END_REF].

Having the four measured signals, these computations lead to seven parameters per signal and thus twenty-eight parameters, for each subject s at each time t. Let p s,j (t), j = 1, ..., 28, denote these parameters. For many health conditions, the trends in vital signs differ among patients. Thus, normal is a relative state for each patient, that differs with age, weight, medical history, etc. Hence, instead of using the raw values of the parameters, a normalization phase is performed according to each patient. Therefore, a first segment of length τ is considered from each signal, x s,i (1 : τ ). It corresponds to the initial state of the patient that is assumed to be its stable one.

Then, the different parameters are computed for the initial segments. Let p (0) s,j denote their initial values. Then, for each real-time measurement, the obtained parameters are normalized, by computing their ratio to the parameters of the first segment p s,j (t)/p (0) s,j . In the following, p s,j would denote the normalized parameter, for simplicity.

Note that, a parameter selection would be performed in the offline phase using the training set to keep only informative one. More details are given in Section 3.4. In the following, let J be the set of indices of the selected parameters, and thus p s,j (t), j ∈ J, are the normalized parameters extracted at real-time t from a subject s signals.

ARDS prediction with Belief functions

The theory of belief functions, also called Dempster-Shafer theory, is introduced for the analysis of imperfect information [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF][START_REF] Shafer | A mathematical theory of evidence[END_REF]. This theory makes it possible to illustrate the uncertainty and imprecision of information. It also takes into account ambiguities and conflicts between sources. The belief functions theory operates on a frame of discernment Ω which consists of a finite set containing mutually exclusive propositions. In this study, there are two groups of subjects, ARDS labeled "+1" and non-ARDS labeled " -1". Then Ω = {+1, -1}. Let 2 Ω be the set of all possible subsets of Ω, it is then defined by:

2 Ω = {∅, {+1}, {-1}, {+1, -1}}.
While ∅ denotes that the state of the subject is neither ARDS, nor non-ARDS, which means impossibility, {+1, -1} denotes ignorance and thus ambiguity.

Basic concept

The main concept relies on the modeling of the evidence provided by the extracted parameters. Information given by a source, that is a parameter, can be represented by a basic belief assignment (BBA), also named mass function [START_REF] Smets | The transferable belief model[END_REF]. Let m j (•) be the mass function associated to the parameter p j . Then, m j (ω, s, t) reflects the strength of evidence supporting a subset ω ∈ 2 Ω according to the value of the parameter p j at time t. In other words, m j (ω, s, t) represents the part of evidence saying that the state of the subject s falls in ω at time t. The mass function should have the following properties ∀s, t:

       m j (∅, s, t) = 0, m j (ω, s, t) → [0, 1], for ω ∈ 2 Ω , ω∈2 Ω m j (ω, s, t) = 1. (1)
In order to construct the mass functions related to each parameter, the training sets are considered. For ARDS subjects of a given training set, the last segments of their signals are used to compute their parameters p s,j = p s,j (T s ), T s being the length of the signals of an ARDS subject s. The reason of this choice is that the end of the signals precedes directly the occurrence of ARDS and thus contains mostly its instability. For non-ARDS subjects, the more the segments approaches the beginning of the signals, the more we are sure of the stability of the subject, since it is possible that a subject labeled non-ARDS develops the ARDS after recordings. Then, the parameters of non-ARDS are computed using segments selected from the wholly beginning of the signals, following the first segments. Let S +1 and S -1 be the sets of indices of ARDS and non-ARDS subjects of the training set respectively. Then, normalized parameters databases p s,j , s ∈ S +1 , and p s ′ ,j , s ′ ∈ S -1 , are obtained. Afterwards, for a parameter type j, the pdfs are estimated using (1) only ARDS parameters values p s,j , s ∈ S +1 , (2) only non-ARDS values p s ′ ,j , s ′ ∈ S -1 and (3) parameters of both classes p s,j , s ∈

S +1 ∪ S -1 . Let Q j,{+1} (•), Q j,{-1} (•) and Q j,{+1,-1} (•) denote the estimated pdfs.
These distributions could be obtained by fitting the histograms of the parameters to all the classical distributions. Then, in the online phase, having a computed parameter p s,j (t) at the real-time t, the mass assigned to each subset of 2 Ω is given by

m j (ω, s, t) = Q j,ω (p s,j (t)) ω ′ ∈2 Ω ,ω ′ =∅ Q j,ω ′ (p s,j (t)) , ω ∈ 2 Ω , ω = ∅, (2)
whereas m j (∅, s, t) = 0.

2 . The first parameter value falls in the distribution of ARDS subjects and the second is in the middle of ARDS and non-ARDS distributions. When assigning the masses of each value from Eq. ( 2), the assigned mass of p (1) 2 from Q j,{+1} is higher than those of Q j,{-1} and Q j,{+1,-1} . Thus, the decision obtained from the parameter p

(1)
2 is the state of {+1}. On the other hand, for parameter p

(2)
2 , the distributions of {+1} and {-1} are too close, leading to similar masses from both. However, higher mass is associated in this case to the distribution of {+1, -1}, thus introducing the ambiguity of such case. The mass assigned to {+1, -1} is higher than those of {+1} and {-1} in all the cases where these distributions are too close. This figure illustrates the effectiveness of the evidence-theory concept by inserting ambiguities when dealing with close data and thus avoiding erroneous decisions.

Discounting

Since we are working on parameters extracted from human bodies, the information cannot be completely reliable. Hence, a discounting operation is needed to take into consideration the reliability of the information provided by each parameter [START_REF] Shafer | A mathematical theory of evidence[END_REF][START_REF] Smets | Belief function: the disjunctive rule of combination and the generalized bayesian theorem[END_REF]. This operation transforms each mass function to a less informative one, based on the degree of reliability of the source. Moreover, in this approach, it is possible to estimate its reliability according to each state, leading to a contextual discounting [START_REF] Mercier | Belief functions contextual discounting and canonical decompositions[END_REF][START_REF] Mercier | Corrigendum to belief functions contextual discounting and canonical decompositions[END_REF]. Let β j,{+1} = 1α j,{+1} and β j,{-1} = 1α j,{-1} be the degrees of reliability of a parameter p j knowing that the true state is " + 1" and " -1" respectively. Then the Probability distributions discounted mass function α m j assigned to p j is given by

         α m j ({+1}, s, t) = β j,{-1} m j ({+1}, s, t), α m j ({-1}, s, t) = β j,{+1} m j ({-1}, s, t),
α m j (Ω, s, t) = m j (Ω, s, t) + α j,{-1} m j ({+1}, s, t) + α j,{+1} m j ({-1}, s, t), α m j (∅, s, t) = 0.

(3) The error rates α j,{+1} and α j,{-1} could be estimated using the distribution functions Q j,ω . Knowing that a subject s state is equal to " + 1", an error is made in evidence assignment if a higher mass is given to any subset of Ω other than {+1}. It is also the same for {-1}. Then, α j,ω , ω = {+1} or {-1}, is computed for all parameter values, where Q j,ω (p) is less than any Q j,ω ′ (p), for ω ′ = ω ∈ 2 Ω . Therefore, it is defined as follows:

α j,ω = ǫ j,ω = Dj,ω Q j,ω (p)dp, (4)
where

D j,ω = {p|Q j,ω (p) < Q j,ω ′ (p), ∀ω ′ ∈ 2 Ω , ω ′ = ω}.

Combination

Once the discounted masses of the different parameters are computed, they can be combined using the normalized conjunctive rule of combination from the Dempster-Shafer theory. By doing this computation, one obtains a more informative mass function, leading to more efficient estimation. The combination of these masses is denoted by ⊕ and defined by the following equation:

m (ω, s, t) = ∩ω ′(j) =ω j∈J (α) m j (ω ′(j) , s, t) 1 - ∩ω ′(j) =∅ j∈J (α) m j (ω ′(j) , s, t)
.

(5)

Decision making

Finally, to make a decision on the health state of a subject s, the mass function m defined on Ω is transformed to a probability measure over the singletons with the pignistic transformation defined by:

BetP (ω, s, t) = ω ′ ∈2 Ω ,ω⊆ω ′ m (ω ′ , s, t) | ω ′ | , ∀ω ∈ Ω, (6)
for ω = {+1} or {-1}. In this way, the final mass of {+1, -1} subset is divided between {+1} and {-1}. The state having the highest pignistic level is selected for each subject at every time t. Thus, a decision is made over the health state of subject s whether he will develop ARDS based on its recordings, and an alert is generated for a positive decision.

To make a final decision from the real-time analysis about the health state of a patient, it is also possible to wait for successive positive decisions. Here, a threshold must be defined as being the number of successive ARDS decisions needed to decide if the patient is going to develop ARDS. This is done by applying the algorithm iteratively on the total length of the signals s of the training set, then different numbers of successive positive decisions are considered and the performance indexes are computed for each number, that are the sensitivity (Se) and specificity (Sp) indexes, also called true positive rate and true negative rate, respectively. The threshold of successive ARDS decisions is obtained thereafter by maximizing the Youden index, Youden = Se+Sp-100, over all the considered numbers.

Waiting for successive positive estimations overcomes detection of islolated abnormalities and leads to a higher specificity of the algorithm.

Selection of parameters

Instead of considering all the extracted parameters in the model, it is relevant to select the most informative ones with less error rates. This part of the algorithm consists of ranking the parameters according to their error rates and then performing a sequential forward selection (SFS) of parameters [START_REF] Devijver | Pattern recognition: A statistical approach[END_REF]. Hence, an error rate ǫ j is computed for each parameter p j , which is the average of conditional errors, ǫ j = (ǫ j,{+1} + ǫ j,{-1} )/2. The conditional errors ǫ j,ω are computed using (4). A global ranking of parameters is performed. Parameters from all signals are ordered from the lowest error obtained to the highest one.

Once the parameters are ranked, a sequential forward selection procedure is performed. One parameter is added at a time and performance indexes of the model are computed. For each added parameter, the belief model is applied in real-time for all the subjects of the training set and a decision is generated for each of them. The sensitivity and specificity indexes are computed. Thereafter, A curve is plotted and updated at each iteration illustrating the sum of the performance indexes obtained versus the number of parameters. Whenever a new value added to the curve shows a decrease in the performance higher than a value ε, the SFS procedure ends and the parameters obtained are maintained.

In addition, in order to take into consideration the dependence of the parameters, the correlation between the added parameter and the existing parameters is tested and 0.9985 ± 0.0178 1.0022 ± 0.0397 compared to the selection procedure alone. If any correlation coefficient between the parameters presents a p value ≤ 0.05 then the parameter is considered correlated with an existing one, thus it is eliminated from the selection procedure.

Results

As mentioned above, two types of parameters were extracted from each time series. Linear parameters are the mean µ, the standard deviation σ, the skewness Sk and the kurtosis Kt. Non linear parameters are the sample entropy SampEn, and both factors from the detrended fluctuation analysis DFA 1 and DFA 2 . For sample entropy, the value of u is taken equal to 5, and r equals 0.2, as mentioned in [START_REF] Fele-Žorž | A comparison of various linear and non-linear signal processing techniques to separate uterine emg records of term and preterm delivery groups[END_REF]. For the detrended fluctuation analysis, DFA 1 is calculated between the range n = 4 and n = 16 and DFA 2 is calculated between n = 16 and n = 64, as shown in [START_REF] Francis | Physiological basis of fractal complexity properties of heart rate variability[END_REF]. The results are obtained by a 5 fold cross validation repeated 10 times. For the illustration of the approach, presented in Sections 4.1, 4.2 and 4.4, only one training and one test sets are considered, then all the sets are considered for the performance evaluation in Section 4.5.

Statistical analysis

First, a statistical analysis is performed to identify parameters that show a significant difference between ARDS and non-ARDS groups. For each parameter, values were compared between both groups using the two-sample F-test. A p-value < 0.02 was considered as significant. Results are presented in Table non-linear parameters for each signal has shown significant difference between groups.

σ R R S k M A B P S k H R S a m p E n H R S a m p E n R R S k S p O 2 S k R R D F A 1 M A B P K t R R D F A 2 M A B P D F A 1 S p O 2 D F A 2 S p O 2 K t R R S a m p E n M A B P µ S p O 2 σ M A B P D F A 2 H R D F A 2 R R
It can be noted from the table that distributions of non ARDS group are more asymmetric than ARDS group for all the signals from the values of skewness. In addition to being asymmetric, kurtosis shows that distributions of non ARDS group for HR, RR and SpO 2 have heavy tails. Moreover, RR's distribution is widely spread in ARDS group according to σ. For ARDS group, sample entropy is higher in HR and RR signals; while it is lower in SpO 2 . Detrended fluctuation analysis factors are generally higher in non ARDS group for SpO 2 and MABP signal.

Parameters Selection

Since there is a high number of parameters, a ranking of these parameters based on the error rates is performed. The ranking is illustrated in Figure 3. Parameters having error rates higher than 0.4 are not presented. The results obtained from the error rate ranking in the figure are confirmed from the statistical analysis shown in Table 1. As shown, most of the parameters that have low errors in Figure 3 also showed a significant difference in the table. This figure shows that the standard deviation extracted from RR signals presented the lowest error rate among all parameters extracted from the four signals, thus ranking at the top.

From these ranked parameters, a sequential feature selection is considered to reduce the input vector of the belief computations. Figure 4 illustrates the selection procedure on the training set for both cases with and without correlations. The plot shows the accuracy as a function of the indicators of the ranked parameters. The accuracy is defined by the proportion of correctly identified subjects among the total number of subjects. A marker is placed on the curves when a parameter is selected. Starting with σ from RR, the parameters presented in Figure 3 the accuracy decreases more than ǫ > 5%. For the correlation study, if a parameter shows a correlation with any existing parameter with a p value lower than 0.05 it is excluded from the selection. The simple selection procedure finished on the insertion of the seventh parameter that is Sk from RR. Thus, six parameters are selected in the simple SFS procedure which are Sk and SampEn from HR, σ and SampEn from RR, Sk from both SpO 2 and MABP. However, there is no decrease in the accuracy of the model in the correlation study but the accuracy is lower than that obtained from the simple selection.

Selection of the alert generation threshold

In order to make a decision on the health state of a new test subject, a threshold must be defined from the training set. Therefore, we sweep the value of successive decisions on the interval [1, ..., 30] hours and the performances Se and Sp are computed. The threshold that maximizes the Youden index is then considered. Figure 5 illustrates the ROC curve of the proposed algorithm by taking the fusion of the first six parameters obtained in the previous section. The plot shows that the optimal cut-off value is obtained for threshold = 20 hours with Se = 81.4% and Sp = 83.78%.

Selection of window size

As described in the real-time model we have used a fixed-length window for extracting the parameters vector for all subjects. Thus, the selection of the size of the window is important in this study. A short window may not contain enough information and a long window may eliminate an abnormal change in the signal. Different window sizes are evaluated in this study, values are considered in the range [4, 6, . . . , 24] hours. threshold). As shown in the figure, window size = 24 hours presents the highest accuracy of 82.5% on the considered training set. These values are then considered for the following.

Performance Evaluation

After the illustration of the extraction and selection of the parameters, the model is tested on the test sets. In this section, results are obtained by a 5-fold cross validation repeated 10 times. Table 2 presents the performance of each added phase to obtain at the end the performance of the complete proposed model. The modeling of the data by associating and combining masses from all the parameters leads to a sensitivity of 70.23% and a specificity of 71.74% over the training sets and Se = 70.16% and Sp = 67.89% over the test sets. However, when the parameter selection phase was performed, an enhancement in the sensitivity was noted over the training sets (Se = 78.10%, Sp = 66.91%) and the test sets (Se = 72.59%, Sp = 61.16%). In addition, discounting the mass functions using the error rates improved the performance to 81.65% of sensitivity and 74.72% of specificity over the training sets and 77.24% and 71.25% for sensitivity and specificity respectively over the test sets. However, when the correlation was considered in the selection phase the performances decreased with Se = 80.18% and Sp = 69.29% over the training sets and 74.10%, 61.12% over the test sets.

Comparison to the state-of-the-art methods

The proposed model is compared to some classification techniques such as Bayesian hypothesis test, k-nearest neighbors (kNN) and SVM. the training set. A ten-fold cross validation was performed to obtain k that minimizes the test loss. The optimal k was found to be 8. As for SVM, it has proven to be a good model to classify data. It can be used to differentiate the data using a decision bound-365 ary that produces the optimal separation of classes. Table 3 shows the performance indexes of the proposed model compared to these described techniques. As shown, the proposed model outperforms the kNN, the SVM as well as the bayesian tests, that are a simplest form of the evidence theory where the ambiguities between the sets are not considered. These results show the importance of the proposed model.

Discussion

This study proposed a real-time model, based on the belief functions theory, for predicting ARDS prior to its onset. ARDS results from inflammatory alveolar injury occurring progressively [START_REF] Thille | Chronology of histological lesions in acute respiratory distress syndrome with diffuse alveolar damage: a prospective cohort study of clinical autopsies[END_REF]. The process begins with the destruction of the barriers of the alveolar membranes resulting by increasing the alveolar permeability. Then, injuries progress causing decreased pulmonary compliance which leads to pulmonary hypertension and refractory hypoxia [START_REF] Sweeney | Acute respiratory distress syndrome[END_REF][START_REF] Patroniti | Clinical management of severely hypoxemic patients[END_REF]. Moreover, ARDS has been linked to clinical features such as severe dyspnea, tachypnea and hypoxemia according to [START_REF] Gajic | Early identification of patients at risk of acute lung injury: evaluation of lung injury prediction score in a multicenter cohort study[END_REF]. Other studies on the clinical profile of ARDS had shown that it is associated with breathlessness, hypotension, tachypnea [START_REF] Chaudhury | A study on clinical profile and outcome of patients with acute respiratory distress syndrome in a tertiary care hospital in north east india[END_REF] and tachycardia [START_REF] Park | Incidence, risk factors, and mortality associated with acute respiratory distress syndrome in combat casualty care[END_REF]. Although, there are no clear clinical features for ARDS, it is evident that ARDS is linked to cardiorespiratory and cardiovascular mechanisms. Therefore, the four vital signs that are considered are the heart rate, the respiratory rate, the oxygen saturation and the mean airway blood pressure. However, there is a lack in the literature for the characterization of ARDS using feature extraction.

In feature extraction techniques, parameters are extracted from raw data signals and pertinent ones are selected to enhance the fusion performance. The characterization of signals using the extraction of parameters has been widely applied in biomedical applications. For instance, an analysis of heart rate variability during syncope was conducted by extracting parameters such as the mean, the standard deviation of the RR interval and non linear parameters such as the sample entropy and the detrended fluctuation analysis [START_REF] Khodor | Kernel based support vector machine for the early detection of syncope during head-up tilt test[END_REF]. Moreover, linear and non-linear parameters were extracted from heart rate variability signals for ARDS subjects to analyze their response to alveolar recruitment manoeuver [START_REF] Borghi-Silva | Analysis of heart rate variability and cardiovascular response in the alveolar recruitment manoeuver in acute respiratory distress syndrome[END_REF]. Only few studies have considered other cardiovascular signals than heart rate, like blood pressure and oxygen saturation signals [START_REF] Rivera | Heart rate and systolic blood pressure variability in the time domain in patients with recent and longstanding diabetes mellitus[END_REF][START_REF] Turianikova | The effect of orthostatic stress on multiscale entropy of heart rate and blood pressure[END_REF]. A characterization of heart beat and blood pressure for diabetes patients using time domain parameters such as the mean, the standard deviation, the skewness and the kurtosis was studied in [START_REF] Rivera | Heart rate and systolic blood pressure variability in the time domain in patients with recent and longstanding diabetes mellitus[END_REF]. In addition, detrended fluctuation analysis was extracted from blood pressure signal to study its dynamics through surgical procedure in rats [START_REF] Galhardo | Detrended fluctuation analysis of a systolic blood pressure control loop[END_REF]. In [START_REF] Jalali | Discovering hidden relationships in physiological signals for prediction of periventricular leukomalacia[END_REF], a characterization of periventricular leukomalacia, a brain injury affecting infants, is proposed by extracting minimum, maximum, mean, variance, skewness, kurtosis, energy of wavelet coefficients and sample entropy from heart rate, mean arterial blood pressure and oxygen saturation. From what preceded, we can conclude that the extraction of parameters from different physiological signals might help in the characterization of physiological aspects related to several pathologies, such as ARDS.

Hence, in this paper, linear and non-linear parameters were extracted from the time series as mean, standard deviation, skewness, kurtosis, sample entropy and detrended fluctuation analysis. Linear parameters provide information regarding the distributions of the data. Distributions were shown to be asymmetric, having wide tails and being widely spread. This verifies that biomedical data are not gaussian [START_REF] Chua | Application of higher order statistics/spectra in biomedical signalsa review[END_REF]. Conversely, it was found that heart rate and respiratory rate had higher sample entropy in ARDS group, thus increased complexity [START_REF] Richman | Physiological time-series analysis using approximate entropy and sample entropy[END_REF]. However, oxygen saturation presented lower complexity in ARDS subjects. In other terms, SpO 2 presents more similarity in the data. From DFA 1 and DFA 2 values, we can notice that SpO 2 and MABP had lower correlations in ARDS subjects [START_REF] Peng | Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series[END_REF].

Moreover, feature fusion is performed using the belief functions framework. The proposed belief functions framework defines ARDS and non-ARDS groups based on the probability distributions of the data according to each parameter. In addition, these models consider supper subsets instead of considering only singletons subsets, thus they introduce the notion of ambiguity between choosing an ARDS or non-ARDS states. The belief theory also takes into account the reliability of each source of information presented by the discounting operation. The reliability of the sources is a very important characteristic, especially when working with parameters extracted from physiological signals. On top of these comes the importance of belief functions framework, that is the modeling of the data even with missing recordings. It often happens that some signals are missing in some part of the recordings because of disconnected electrodes or a sudden movement of the patient. Unlike classical data mining techniques, the belief functions guarantee that missing parameters due to missing signals will not affect the overall performance of features fusion.

From the results, the reduction of the number of parameters, by performing a selection procedure, improved the sensitivity of the model and reduced the complexity of the computations as well. This way only parameters having low error rates are included in the model and the combination that presents the best local accuracy is considered. Thus, the modeling of the data is done using the minimum amount of parameters instead of using all the extracted parameters. Then, the inclusion of the conditional reliability of each source has led to an enhancement in the accuracy of the model over the training and the test sets. Each parameter gives information regarding the patient and the reliability of the parameter itself in distinguishing between states. Finally, 77% of ARDS patients are correctly identified in real-time before the occurrence of ARDS and only 29% of non-ARDS patients are misclassified and considered as possible developers of ARDS. The obtained results demonstrate, beside the already existing studies on ARDS, that ARDS is linked to cardiorespiratory and cardiovascular mechanisms. More importantly, ARDS can be predicted using non-invasive physiological signals instead of using clinical exams or variables that necessitate the presence of subjects in hospitals for long period [START_REF] Navarrete-Navarro | Acute respiratory distress syndrome in trauma patients: Icu mortality and prediction factors[END_REF][START_REF] Ichikado | Prediction of prognosis for acute respiratory distress syndrome with thin-section ct: Validation in 44 cases 1[END_REF][START_REF] Ennett | Predicting respiratory instability in the ICU[END_REF].

It is shown from the results that this model outperformed the classical classifiers tested similarly on the data. It is worth noting that the bayesian tests perform only on singletons and the state membership is determined by computing the likelihood ratio. However, the added ambiguity set in the proposed evidence-based model may affect the final decision since it is considered in both discounting and combination phases. Thus, the advantage of such models over the classical classifiers, especially the bayesian tests, is the inclusion of higher order sets that illustrate the ambiguity intervals.

Previous studies were conducted on the prediction of ARDS or the mortality related to ARDS using computed tomographic images [START_REF] Ichikado | Prediction of prognosis for acute respiratory distress syndrome with thin-section ct: Validation in 44 cases 1[END_REF] or clinical variables [START_REF] Navarrete-Navarro | Acute respiratory distress syndrome in trauma patients: Icu mortality and prediction factors[END_REF][START_REF] Ennett | Predicting respiratory instability in the ICU[END_REF], that necessitate clinical exams or the presence of subjects in the hospitals for long periods.

However, the approach proposed in this paper uses only vital signs that can be collected using non-invasive techniques, thus without disturbing the subjects. Moreover, these studies present general statistical models to detect ARDS assuming that all subjects will present similar changes; while this approach is a subject-based one that considers an initial stable state for each subject and computes the ongoing parameters with respect to the stable state. In addition, an information fusion was handled using the evidence theory to consider the ambiguity and unreliability of the sources. Finally, this approach achieved good performances for both sensitivity and specificity indexes for the realtime surveillance of ARDS.

Conclusion

This paper proposes an evidence-based model for the prediction of ARDS in intensive care unit patients. This work presents several contributions. Non invasive vital signs are included in the study for many reasons, the facility of their acquisition, the possible integration in home surveillance systems and the link between these signals and the risk factors associated with ARDS. Linear and non-linear parameters were computed and considered the source of information, since they can give information about the properties of a signal more than the signal itself. Imprecision and unreliability of information sources were considered by the belief functions model. It assigns masses to each subset computed from a predefined training set. Then, these masses were discounted and combined according to the reliability of each source of information. This model is then extended to reduce the dimension of the input by performing a parameter selection procedure. This model has achieved high performances in both ARDS and non-ARDS groups. In addition, it outperforms state of the art techniques. Further work must be done to extract more parameters that provide evidence on different characteristics in the signals. Moreover, it would be interesting to establish the relationships between the changes in time series and the cardiovascular mechanism. A further study on the parameters selection will also be done, to take into consideration the dependance of parameters from the beginning of the selection process.

Figure 1 :

 1 Figure 1: An example of the extracted time series for an ARDS subject

  Figure 2 illustrates an example of such a computation. The plot shows the distributions of the normalized standard deviation of respiratory rate signals parameter for ARDS subjects of a training database in straight thick line, the one for non-ARDS subjects in dashed line and the one of all subjects in straight thin line. Moreover, the plot shows the mass assignment for two new parameter values p

2 Figure 2 :

 22 Figure 2: Mass assignment of new computations in the case of the normalized standard deviation of respiratory rate signals.

Figure 3 :

 3 Figure 3: Ranking of parameters according to error rates. Only parameters having an error lower than 0.4 are included in the figure.

Figure 4 :

 4 Figure 4: The parameters selection procedure using the simple procedure and the correlation extension procedure.

Figure 6 Figure 5 :

 65 Figure 5: ROC curve for the selection of the threshold for alert generating.

  370

Table

Table 1 :

 1 The most significant parameters for each signal, parameters having a p-value < 0.02 are shown. Mean ± standard deviation are shown for each parameter.

	Signal Parameters	ARDS group	non ARDS group
		mean	1.0297 ± 0.1439	1.0218 ± 0.0935
	HR	Sk Kt	-1.875 ± 6 1.0628 ± 0.4082	-3.3629 ± 31.2769 1.2803 ± 0.6892
		SampEn	1.2277 ± 0.6079	1.0231 ± 0.1371
		σ	1.6226 ± 1.4729	0.9209 ± 0.2851
	RR	Sk	-0.8592 ± 10.6885	1.0231 ± 3.3186
		Kt	1.0172 ± 0.5052	1.3381 ± 1.2061
		SampEn	1.4088 ± 0.5191	1.0281 ± 0.2
		Sk	0.6881 ± 4.0655	3.1707 ± 14.2335
	SpO 2	Kt SampEn	1.0871 ± 0.5252 1.2836 ± 0.6092	1.0907 ± 0.8359 1.5598 ± 0.982
		DF A 1	0.8678 ± 0.1777	1.0172 ± 0.2767
		DF A 2	1.0289 ± 0.1349	1.0194 ± 0.0667
	MABP	µ Sk	1.0262 ± 0.1246 0.6396 ± 1.6238	0.9868 ± 0.072 0.5796 ± 8.2055
		DF A 1	0.9432 ± 0.1459	1.0743 ± 0.3192
		DF A 2

  . Different linear and

		0.45
		0.4
		0.35
	Error rates	0.25 0.3
		0.2
		0.15
		0.1

Table 2 :

 2 The Bayesian hypothesis test is based on Bayes theorem. It is a probability model that consists of factorizing a joint probability distribution into a set of conditional distributions for each variable. The knearest neighbors is one of the basic algorithms in machine learning. It classifies a new observation to a class by choosing the majority class among the k nearest neighbors in Changes in the accuracy over the training set with windows length. Values in the figure denote (number of parameters, threshold). Performance of the proposed model and the influence of each phase.

		85									(6,20)
	Accuracy (%)	65 70 75 80	(14,12)	(8,12)	(8,18)	(4,20)		(7,18)	(8,24)	(16,1)	(16,1)
							(10,24)
		60
		(13,24)
		4 55	6	8	10	12	14	16	18	20	22	24
					window length (in hours)
	Figure 6: Method		training set Se(%) Sp(%) Se(%) Sp(%) test set
		BF without discounting, nor selection 70.23 71.74 70.16 67.89
			BF without discounting	78.10 66.91 72.59 61.16
			complete BF model		81.65 74.72 77.24 71.25
			complete BF model+corr	80.18 69.29 74.10 61.12

  3: Comparison to state-of-the-art classification methods

	Method	Sensitivity (%) Specificity (%)
	KNN	59.20	79.75
	Bayesian test	67.74	65
	SVM	58.67	74.06
	The proposed Model	77.24	71.25

Acknowledgment

The authors would like to thank Université de Technologie de Troyes, France and the Lebanese University, Lebanon for their funding of this research.


Rhys M Adams

Justin B Kinney

Aleksandra M Walczak
email: awalczak@lpt.ens.fr

Thierry Mora
email: tmora@lps.ens.fr

Epistasis in a Fitness Landscape Defined by Antibody-Antigen Binding Free Energy

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

To ensure a reliable response and to neutralize foreign pathogens, the adaptive immune system relies on affinity maturation. In this process, antibody receptors expressed by B cells undergo an accelerated Darwinian evolution through random mutations and selection for affinity against foreign epitopes [START_REF] Cobey | The evolution within us[END_REF]. Mature antibodies can accumulate up to 20% hypermutations from their germline sequence [START_REF] Marcou | High-throughput immune repertoire analysis with IGoR[END_REF], leading to up to a 10,000 fold improvement in binding affinity [START_REF] Eisen | Variations in Affinities of Antibodies during the Immune Response[END_REF]. Affinity maturation also produces broadly neutralizing antibodies that target conserved regions of the pathogen, of particular importance for vaccine design against fast evolving viruses [START_REF] Corti | Broadly Neutralizing Antiviral Antibodies[END_REF]. Despite extensive experimental and theoretical work, the key determinants of antibody specificity and evolvability are still poorly understood, mainly because the sequence-to-affinity relationship is difficult to measure comprehensively or to predict computationally [START_REF] Esmaielbeiki | Progress and challenges in predicting protein interfaces[END_REF].

A major confounding factor in characterizing the sequence dependence of any protein function, including affinity, is the pervasiveness of epistasis, the phenomenon by which different mutations interact with each other [START_REF] Phillips | Epistasis -The essential role of gene interactions in the structure and evolution of genetic systems[END_REF]. Theory [START_REF] Carter | The role of epistatic gene interactions in the response to selection and the evolution of evolvability[END_REF][START_REF] Good | The impact of macroscopic epistasis on long-term evolutionary dynamics[END_REF][START_REF] Paixão | The effect of gene interactions on the longterm response to selection[END_REF] and genomic data [START_REF] Breen | Epistasis as the primary factor in molecular evolution[END_REF] suggest that inter-and intragenic epistasis plays a major role in molecular evolution, by constraining the set of accessible evolutionary trajectories towards adapted phenotypes [START_REF] Weinreich | Darwinian Evolution Can Follow Only Very Few Mutational Paths to Fitter Proteins[END_REF][START_REF] Poelwijk | Empirical fitness landscapes reveal accessible evolutionary paths[END_REF][START_REF] Gong | Stability-mediated epistasis constrains the evolution of an influenza protein[END_REF][START_REF] Anderson | Intermolecular epistasis shaped the function and evolution of an ancient transcription factor and its DNA binding sites[END_REF][START_REF] Podgornaia | Protein evolution. Pervasive degeneracy and epistasis in a protein-protein interface[END_REF], enhancing evolvability through stabilizing mutations [START_REF] Bloom | Protein stability promotes evolvability[END_REF][START_REF] Bloom | Permissive Secondary Mutations Enable the Evolution of Influenza Oseltamivir Resistance[END_REF], or slowing down adaptation by the law of diminishing returns [START_REF] Chou | Diminishing Returns Epistasis Among Beneficial Mutations Decelerates Adaptation[END_REF][START_REF] Kryazhimskiy | Global Epistasis Makes Adaptation Predictable Despite Sequence-Level Stochasticity[END_REF]. Evidence for epistasis in antibody affinity include direct observations of cooperativity between mutations [START_REF] Midelfort | Substantial energetic improvement with minimal structural perturbation in a high affinity mutant antibody[END_REF][START_REF] Koenig | Deep Sequencing-guided Design of a High Affinity Dual Specificity Antibody to Target Two Angiogenic Factors in Neovascular Age-related Macular Degeneration[END_REF], the dependence of mutational effects on sequence background [START_REF] Boyer | Hierarchy and extremes in selections from pools of randomized proteins[END_REF], and statistical co-variation of residues in large sequence datasets [START_REF] Mora | Maximum entropy models for antibody diversity[END_REF][START_REF] Asti | Maximum-Entropy Models of Sequenced Immune Repertoires Predict Antigen-Antibody Affinity[END_REF].

Intragenic epistasis has mostly been studied either by measuring the fitness of all possible mutational intermediates between two variants [START_REF] Weinreich | Darwinian Evolution Can Follow Only Very Few Mutational Paths to Fitter Proteins[END_REF][START_REF] Schenk | Patterns of epistasis between beneficial mutations in an antibiotic resistance gene[END_REF][START_REF] Szendro | Quantitative analyses of empirical fitness landscapes[END_REF][START_REF] De Visser | Empirical fitness landscapes and the predictability of evolution[END_REF], or by comparing the effect of mutations in different backgrounds [START_REF] Jacquier | Capturing the mutational landscape of the beta-lactamase TEM-1[END_REF][START_REF] Bank | A Systematic Survey of an Intragenic Epistatic Landscape[END_REF][START_REF] Boyer | Hierarchy and extremes in selections from pools of randomized proteins[END_REF]. Many such studies rely on a particular measure of fitness rather than a well-defined physical phenotype. Deep mutational scans (DMS) [START_REF] Fowler | Deep mutational scanning: a new style of protein science[END_REF][START_REF] Sarkisyan | Local fitness landscape of the green fluorescent protein[END_REF] can comprehensively map out the epistatic landscape of many genetic variants [START_REF] Araya | A fundamental protein property, thermodynamic stability, revealed solely from large-scale measurements of protein function[END_REF][START_REF] Olson | A Comprehensive Biophysical Description of Pairwise Epistasis throughout an Entire Protein Domain[END_REF][START_REF] Podgornaia | Protein evolution. Pervasive degeneracy and epistasis in a protein-protein interface[END_REF]. However, most DMS methods do not measure the biophysical quantity of interest directly [START_REF] Vodnik | Phage Display: Selecting Straws Instead of a Needle from a Haystack[END_REF], introducing both nonlinearities and noise that could be misinterpreted as epistasis.

Here we analyze the detailed epistatic landscape of an antibody's binding free energy, which we define as the logarithm of the Tite-Seq measurement of the dissociation constant, to its cognate antigen (the 4-4-20 antibody fragment against fluorescein), using data previously obtained by Tite-Seq, a recently introduced DMS variant that accurately measures protein binding affinity in physical units of molarity [START_REF] Adams | Measuring the sequence-affinity landscape of antibodies with massively parallel titration curves[END_REF]. By comparing to a simple additive model of mutations on the binding free energy, and carefully controlling for measurement noise and nonlinearities, we find that epistasis significantly contributes to the antibody's affinity. This epistasis is not uniformly distributed, but instead favors certain residue pairs across the protein. We use our results to analyze how epistasis both constrains and enlarges the set of possible evolutionary paths leading to high-affinity sequences.

Results

Position Weight Matrix model of affinity

We analyzed data from [START_REF] Adams | Measuring the sequence-affinity landscape of antibodies with massively parallel titration curves[END_REF]) (https://github.com/jbkinney/16_titeseq), where Tite-Seq was applied to measure the binding affinities of variants of the 4-4-20 fluorescein-binding scFv antibody, henceforth called 'wildtype'. Libraries were generated by introducing mutations to either the CDR1H or CDR3H domains restricted to 10 amino acid stretches called 1H and 3H (Fig. 1A). All single amino acid mutants, 1100 random double amino acid mutants, and 150 triple amino acid mutants were generated in multiple synonymous variants and measured, (Fig. 1B). Using a combination of yeast display and high-throughput sequencing at various antigen concentrations, Tite-Seq yielded the binding dissociation constant K d (in M or mol/L) of each variant with the fluorescein antigen.

We first tried to predict the K d of double and triple mutants from single mutant measurements. Mutations are expected to act on the binding free energy in an approximately additive way [START_REF] Wells | Additivity of mutational effects in proteins[END_REF][START_REF] Olson | A Comprehensive Biophysical Description of Pairwise Epistasis throughout an Entire Protein Domain[END_REF]. One may thus write the free energy of binding,

F = ln K d c 0 (
) (defined up to constant in units of k B T), as a sum over mutations in the mutagenized region, s = s 1 ,..., s l

( ) : F s ( ) ≈ F PWM s ( ) = F WT + h i s i ( ) i=1 l ∑ (1)
where F WT is the wildtype sequence energy, and h i s i is the effect of a mutation at position i to residue s i . The elements of the Position- Weight Matrix (PWM) h i s i ( ) are obtained from the K d of single mutants shown in Fig. 1C. Since Tite-Seq measurements are limited to values of K d ranging from 10 -9.5 to 10 S1).

The unexplained variance missed by the PWM model may have four origins: convolution with expression, nonlinear effects, measurement noise, epistasis. Tite-seq was developed specifically to separate the measurement K d and expression. From Pearson's correlation between expression and log(K d ), we find that expression explains 6 and 12% of the R 2 for CDR1H and CDR3H, respectively (Fig. S2). Furthermore, the residual from the PWM prediction, F -F PWM , had almost no correlation with expression (R 2 <0.1% for 1H, R 2 <0.2% for 3H, see Fig. S2). Log transformed expression values yielded similar but smaller contributions. The second "nonlinear effects" case corresponds to the hypothesis of additivity not being valid for F = ln K d c 0

(

) , but for some other nonlinear transformation of F.

Such a nonlinearity, also called "scale," can lead to spurious epistasis [START_REF] Fisher | The Correlation between Relatives on the Supposition of Mendelian Inheritance[END_REF][START_REF] Phillips | Epistasis -The essential role of gene interactions in the structure and evolution of genetic systems[END_REF]. We first checked that additivity did not apply to the untransformed dissociation constant, K d : a PWM model learned from K d instead of F could only explain 34% of the variance of all 1H and 3H multiple mutants, down from 62% when learning from F (Fig. S3). Refitting Tite-seq values with no boundary constraints yielded much worse PWM models, largely attributable to poor estimates of poorly binding antibodies (Fig. S4). We then looked for the non-linear transformation E(F) that would give the PWM model with the best predictive power (Methods and Fig. S5). This optimization yielded only a modest improvement to 65% of the explained variance. In addition, the optimal function E was very close to the logarithm (R 2 =97%, Fig. S6). Since nonlinear effects do not play a significant role, henceforth we only consider the PWM model defined on the free energy.

Epistasis affects affinity

To identify epistasis, we estimated the difference between the measured binding free energies of double and triple mutants, F(s), and the PWM prediction, F PWM s ( ) (Fig. 2A). However, these small differences can be confounded by measurement noise (see Fig. S7), which can be mistaken for epistasis. To control for this noise, we defined Z-scores between two estimates of the free energy, F a and F b ,

as Z = F a -F b ( ) / σ a 2 +σ b 2
, where σ a 2 and σ b 2 are their estimates of uncertainty. Uncertainty was either measured as variance from replicate measurements and synonymous mutations, or as the sum of variances from additive PWM contributions depending on context. We first computed Z-scores between independent estimates of the same free energy using synonymous variants ( Z error , Methods). Excluding mutants at the reliable readout boundary (10 -9.5 M ≤ K d ≤10 -5 M ), we found that the distribution of Z error was normal with variance ≈1 (Fig. 2B, orange line), as expected from Gaussian measurement noise.

A comparison between Tite-seq measurements and clones K d measured from flow cytometry was also approximately normal with variance ≈1 (Fig. S8), meaning that Tite-Seq introduced no systematic errors in addition to those estimated from replicates.

We then estimated the effect of epistasis by calculating Z-scores ( Z epi ) from the difference between the PWM prediction, F PWM (Eq. 1), and the measured F. The resulting distributions of Z-scores (Fig. 2A, blue and red lines) had much larger variances than expected from measurement noise (standard deviation 1.76 for 1H, and 3.18 for 3H), indicating strong epistasis. These epistatic effects were on average slightly beneficial (positive Z): 18% of double mutants inside the reliable readout boundaries (10 -9.5 M ≤ K d ≤10 -5 M ) showed significant beneficial epistasis ( Z epi >1.64 , corresponding to p<0.05 in a one-sided Z-test), and 12% significant deleterious epistasis ( Z epi < -1.64).

Comparing the variance of Z epi with that of Z error gives a large fraction of the unexplained variance that is attributable to epistasis, 1-Var Z error ( ) / Var Z epi ( ) = 60% for 1H, and 88% for 3H. While clones at the reliable readout boundaries under-estimated measurement error, their inclusion yielded more extreme results (Fig. S9). PWMs trained from optimal transformations had almost no effect on epistasis estimates (Fig. S10).

To determine whether certain positions along the sequence concentrated epistatic effects, we computed the mean squared Z-score for all double mutations at each pair of positions (excluding median boundary values), revealing a complex and heterogeneous landscape of epistasis (Fig. 2C, see Fig. S11 for the epistasis magnitude superimposed on the wildtype's crystal structure). CDR3H, which interacts directly with the antigen, is observed to have more epistatically interacting sites than CDR1H. Interestingly, the three most epistatic pairs in 3H -between positions 101, 106 and 108 -are mutated in the previously described super-optimized 4m5.3 antibody [START_REF] Boder | Directed evolution of antibody fragments with monovalent femtomolar antigen-binding affinity[END_REF] (mutations shown in green in Fig. 1B), consistent with previous suggestions that positions 101 and 106 interact together and with position 108 via hydrogen bonds [START_REF] Midelfort | Substantial energetic improvement with minimal structural perturbation in a high affinity mutant antibody[END_REF][START_REF] Adams | Measuring the sequence-affinity landscape of antibodies with massively parallel titration curves[END_REF]. Epistasis is usually expected between residues that are in contact in the protein structure [START_REF] Romero | Navigating the protein fitness landscape with Gaussian processes[END_REF][START_REF] Morcos | Direct-coupling analysis of residue coevolution captures native contacts across many protein families[END_REF][START_REF] Mclaughlin | The spatial architecture of protein function and adaptation[END_REF][START_REF] Zhang | Evolution of protein structures and interactions from the perspective of residue contact networks[END_REF][START_REF] Melamed | Deep mutational scanning of an RRM domain of the Saccharomyces cerevisiae poly(A)-binding protein[END_REF], as for instance between positions 101 and 106. However, the mean squared Z-score is only weakly correlated with residue distance (r=-0.13,p=0.22 for 1H, r=-0.27,p=0.022 for 3H, Fig. S12). Additionally, while distance to antigen have been shown to predict how strongly mutations affect binding affinity [START_REF] Brenke | Application of asymmetric statistical potentials to antibody-€"protein docking[END_REF][START_REF] Kepler | Immunoglobulin Gene Insertions and Deletions in the Affinity Maturation of HIV-1 Broadly Reactive Neutralizing Antibodies[END_REF], we did not detect a strong relationship between epistatic contributions and distance to antigen (r=-0.24,p=0.511 for 1H, r=0.19,p=0.603 for 3H, Fig. S12F). This may be due in part to this study's mutated region being too close to the antigen to detect an association: 12 out of the 20 mutated residues are within 10 Å of the antigen, and all 20 are within 16 Å.

We next looked for evidence of "sign epistasis," where one mutation reverses the sign of the effect of another mutation (Fig. 2A). Sign epistasis can constrain evolution by blocking paths to fit sequences [START_REF] Weinreich | Darwinian Evolution Can Follow Only Very Few Mutational Paths to Fitter Proteins[END_REF][START_REF] Poelwijk | Reciprocal sign epistasis is a necessary condition for multi-peaked fitness landscapes[END_REF][START_REF] Weinreich | Perspective:sign epistasis and genetic constraint on evolutionary trajectories[END_REF]. We defined a Z-score for a single mutation A quantifying the beneficial effect of that mutation relative to the noise,

Z A = F WT -F A ( ) / σ A ,
where F WT and F A are the wildtype and mutant free energies, and σ is the measurement error estimated as before. Since we are only interested in the sign of the effect, we kept single mutants at the reliable readout boundary. An equivalent Z-score was defined for a mutation A in the background of an existing mutation B:

Z A|B = F B -F AB ( ) / σ A 2 +σ AB 2
, where F AB is the free energy of the double mutant AB. Significant sign epistasis was defined by Z A|B Z A < 0 and Z A|B , Z A > 1.64, and reciprocal sign epistasis by the additional symmetric condition A ↔ B .

With a 5% false discovery rate (Benjamin Hochberg procedure) we found 52 significant sign epistasis examples. These are listed in S1_table_sign_epistasis.csv and summarized in Tables S1 andS2. Deleterious sign epistasis was exceptional, with 3 instances in 1H and 6 in 3H. These cases, as well as the four most significant cases of benefecial sign epistasis are depicted in Fig. S13. These mutants represent evolutionary trajectories blocked due to the poor binding affinity of their single mutations. Among cases where both single mutations were deleterious, we found 4% ( p < 10 -15 , Binomial test) of mutants in 1H and 0.8% ( p < 10 -7 , Binomial test) in 3H with significant beneficial epistasis, versus 0.06% expected by chance (the null expectation, which takes into account the constraint that Z A + Z B|A = Z B + Z A|B , is defined in the Methods); 1% ( p < 10 -15 , Binomial test) were reciprocal in 1H, and 0.4% ( p < 10 -10 , Binomial test) in 3H, versus 0.01% expected by chance. To evaluate how these epistatic interactions may affect affinity maturation, we estimated how often "viable" double mutants were separated from the wildtype by nonviable single mutants, where viability is defined by K d < 10 -6 M [START_REF] Batista | Affinity Dependence of the B Cell Response to Antigen: A Threshold, a Ceiling, and the Importance of Off-Rate[END_REF][START_REF] Foote | Kinetic and affinity limits on antibodies produced during immune responses[END_REF][START_REF] Roost | Early high-affinity neutralizing anti-viral IgG responses without further overall improvements of affinity[END_REF], forming possible roadblocks to affinity maturation. This strong instance of "rescue" epistasis occurred in roughly half of the mutants with beneficial sign epistasis (Table S1 andS2).

Modeling epistasis and its impact on affinity maturation

To integrate the observed epistatic interactions into a predictive model of affinity, we introduced a model of binding free energy as:

F s ( ) ≈ F pairwise s ( ) = F PWM s ( ) + J ij s i , s j ( ) i< j ∑ , ( 2
)
where J ij is the interaction strength between residues at positions i and j. To avoid overfitting and allow for independent validation (in the absence of a sufficient number of triple mutants), we grouped residues into 4 biochemical categories [START_REF] Voet | Biochemistry, 4th Edition[END_REF]) (polar, nonpolar, acidic, basic, see Methods) and let the entries of J only depend on that category.

We trained the model on the 1208 1H or 1216 3H double and triple mutants, using a Lasso penalty to control for overfitting. The optimal penalty was set by 10 fold cross-validation, i.e. by maximizing the explained variance of a subset comprising 1/10 of the mutants by using a model trained on the remaining 9/10, averaged over the 10 subsets (Fig. S14A andMethods). Interacting pairs with posterior probabilities >0.95 as determined by Bayesian Lasso [START_REF] Park | The Bayesian Lasso[END_REF] are shown in Figs. 2D andE.

Out of the 360 possible terms, 52 1H and 45 3H interaction terms were identified by this method. Although these interactions, whose number is limited by the number of measured variants, only modestly improved the explained variance relative to the PWM in all multiple mutants (from 62% to 64% for 1H and from 58% to 60% for 3H), it substantially improved the affinity prediction of the mutants with significant epistasis (R 2 from 27% to 50% in 1H, from 13% to 44% for 3H, Fig. S14B-C). In contrast, a null linear regression model based on the Miyazawa-Jernigen matrix had negligible improvement on R 2 (see Fig. S15) [START_REF] George | Effects of thymic selection on T cell recognition of foreign and tumor antigenic peptides[END_REF] Notably, two mutations of the superoptimized 4m5.3 antibody are predicted by the model to have epistatic interactions: a slightly deleterious effect between A 101 and L 108 , and a strongly beneficial one between S 102 and L 108 . While these results show some generalizability of biochemical properties, the corresponding model only accounts for a small faction of the variance explainable by epistasis. A more sophisticated approach may be warranted for fully predicting epistatic contributions.

Next we used our models to estimate the diversity, or "degeneracy", of antibodies with good binding affinity. Specifically, we evaluated the degeneracy volume V of high-affinity sequences as the number of sequences with K d <B, using either the PWM (Eq. 1) or pairwise (Eq. 2) models, using a combination of exhaustive and Monte-Carlo sampling (Methods). Compared to the coarse-grained pairwise model trained previously, the interaction strength J was learned directly for each residue pair, without grouping by biochemical category and with no Lasso penalty. The volume of 1H mutants was larger than that of 3H mutants (Fig. 2F), in agreement with the fact that CDR3H plays a more important role in binding affinity. Epistasis increased the recognition volume for both domains, consistent with the previous observation that epistatic effects are, on average, more beneficial than deleterious. To explore the diversity of evolutionary paths leading to recognition, we computed the neutral mutational flux A in and out of the high-affinity region as the probability that a random mutation in a high-affinity sequence (K d < B) causes loss of recognition (K d > B), summed over all high-affinity sequences (Methods). Again, our models show increased mutational flux due to epistasis, even after normalizing by volume, A/V (Fig. 2G). The effect is small, but only reflects the impact of epistasis from the limited, randomly chosen set of double mutations that we measured, which comprises only 7% of all possible double mutations of a 10 amino acid sequence (1058 for CDR1H, and 1066 for CDR3H, out of 16245). We speculate that differences in flux arising from all epistatic interactions may be up to 15 fold stronger. Adding explicit selection to the mutational model would also affect the results, but would require to make additional assumptions about how binding affinity and selection are linked. The neutral mutational flux analysis allows for comparisons that do not depend on such an explicit model of selection.

We checked that these differences were robust to sampling noise and overfitting by performing a jackknife analysis ( p < 10 -5 for the difference in A and V between the PWM and pairwise models, see Methods), and verified that similar conclusions were obtained based on the optimized non-linear transformation (Fig. S16).

Discussion

By analyzing massively parallel affinity measurements obtained by Tite-Seq, we obtained a detailed picture of epistasis in a well-defined physical phenotype -the binding free energy of an antibody to an antigen. Here we define the free energy as the logarithm of the dissociation constant as measured by Tite-Seq and the yeast display assay. While our analysis of epistasis strictly applies to that measured phenotype and not directly to the free energy, these two quantities are believed to be equivalent [START_REF] Van Antwerp | Fine affinity discrimination by yeast surface display and flow cytometry[END_REF]. We showed that antibody sequences contain many epistatic interactions, and that many of these interactions increase affinity. Our approach involves first training an additive (PWM) model as a baseline, and identifying departures from that model as epistasis. In this comparison, a crucial step was to correct for the two issues of scale and measurement noise.

The first issue, identified by Fisher [START_REF] Fisher | The Correlation between Relatives on the Supposition of Mendelian Inheritance[END_REF] and also called unidimensional epistasis [START_REF] Szendro | Quantitative analyses of empirical fitness landscapes[END_REF], is the idea that an epistatic trait becomes additive upon a different parametrization (Sailer and Harms, 2017a). For instance, protein stability, which often determines fitness, is a nonlinear function of the folding free energy difference, which is expected to be roughly additive [START_REF] Bloom | Thermodynamic prediction of protein neutrality[END_REF][START_REF] Bershtein | Robustness-epistasis link shapes the fitness landscape of a randomly drifting protein[END_REF][START_REF] Jacquier | Capturing the mutational landscape of the beta-lactamase TEM-1[END_REF][START_REF] Gong | Stability-mediated epistasis constrains the evolution of an influenza protein[END_REF][START_REF] Serohijos | Merging molecular mechanism and evolution: Theory and computation at the interface of biophysics and evolutionary population genetics[END_REF][START_REF] Bank | A Systematic Survey of an Intragenic Epistatic Landscape[END_REF][START_REF] Sarkisyan | Local fitness landscape of the green fluorescent protein[END_REF]. This leads to both a law of diminishing returns [START_REF] Bank | A Systematic Survey of an Intragenic Epistatic Landscape[END_REF] and robustness to mutations when the protein is very stable [START_REF] Bloom | Thermodynamic prediction of protein neutrality[END_REF]. To disentangle these potential artifacts, we defined our PWM on the binding free energy, which we expect to be additive in sequence content, and we checked that this parametrization was close to minimizing epistasis.

To tackle the second and perhaps more important issue of noise, especially in the context of deep mutational scans where many variants are tested [START_REF] Araya | A fundamental protein property, thermodynamic stability, revealed solely from large-scale measurements of protein function[END_REF], we developed a robust methodology based on Z-scores to identify epistatic interactions as significant outliers. This analysis showed that the variability in binding free energy consists of 60% of additive effects, 25-35% of epistatic effects, and the rest of experimental noise, making the epistatic contribution to the phenotype substantial compared to that of single mutations. A large fraction of that epistasis was beneficial, in contrast with previous reports of mostly negative epistasis owing to the concavity of the scale [START_REF] Bershtein | Robustness-epistasis link shapes the fitness landscape of a randomly drifting protein[END_REF][START_REF] Schenk | Patterns of epistasis between beneficial mutations in an antibiotic resistance gene[END_REF][START_REF] Bank | A Systematic Survey of an Intragenic Epistatic Landscape[END_REF], which we here circumvent by directly considering the free energy.

Epistasis is key to understanding the predictability and reproducibility of evolutionary paths [START_REF] Lässig | Predicting evolution[END_REF][START_REF] Kryazhimskiy | Global Epistasis Makes Adaptation Predictable Despite Sequence-Level Stochasticity[END_REF]. Previous studies have shown that much of the unexplained phenotypic variance could be explained by secondorder epistatic terms, although specific evolutionary trajectories may be sensitive to higher order epistatic terms [START_REF] Sailer | High-order epistasis shapes evolutionary trajectories[END_REF][START_REF] Poelwijk | Learning the pattern of epistasis linking genotype and phenotype in a protein[END_REF]. Our results show how second-order epistasis could constrain the space of possible hypermutation trajectories during affinity maturation, with important consequences for antibody and vaccine design, as the importance of eliciting responses of antibodies that are not just strongly binding but also evolvable is being increasingly recognized [START_REF] Wang | Manipulating the selection forces during affinity maturation to generate crossreactive HIV antibodies[END_REF]. Targeting epistatic interactions may provide an alternative strategy for optimizing antibody affinity: among the 2 epistatic hotspots in CDR1H and 11 in CDR3H that we identified ( Z epi 2 1 2 > 3), 4 involved positions mutated in the super-optimized 4m5.3 antibody sequence, with a higher epistatic contribution than expected by chance. This is consistent with a previous study where an antibody with multiple conformations acquired mutations that stabilized the antibody structure, resulting in a single conformation [START_REF] Wedemayer | Structural Insights into the Evolution of an Antibody Combining Site[END_REF]. We also identified 3 cases of beneficial sign epistasis, in which the double mutant was fit despite the single mutant being deleterious. For instance, the D108E mutations in 4m5.3 is deleterious by itself but is rescued beyond the wildtype value by the S101A mutation [START_REF] Midelfort | Substantial energetic improvement with minimal structural perturbation in a high affinity mutant antibody[END_REF], which occurred first in the directed evolution process [START_REF] Boder | Directed evolution of antibody fragments with monovalent femtomolar antigen-binding affinity[END_REF]. We report 15 extreme cases of viable (K d <10 -6 M) double mutants whose single-mutant intermediates are nonviable, possibly blocking affinity maturation. However, our analysis of the volume and mutational flux of the region of low binding free energies in sequence space suggests that epistasis facilitates the evolution of high-affinity antibodies (Fig. 2F,G). Therefore, we speculate that interactions with the non-mutated parts of the sequence and evolution of the antigen binding partner can either add further constraints or open up additional paths.

Antibodies pose unique questions about the evolvability of evolution [START_REF] Wagner | Perspective: Complex Adaptations and the Evolution of Evolvability[END_REF]. What sort of sequence space would favor quick accessibility while maintaining a small number of viable sequences [START_REF] Perelson | Theoretical studies of clonal selection: Minimal antibody repertoire size and reliability of self-non-self discrimination[END_REF]? A landscape could contain small, permeable, easily evolvable recognition spaces, or could contain large, hard to access, and robust sequences depending on the prevalent epistatic contributions [START_REF] Carter | The role of epistatic gene interactions in the response to selection and the evolution of evolvability[END_REF]. Such a bias could suggest that the tradeoff of an epistatic landscape would include fast initial evolution, followed by slow incremental evolution [START_REF] Good | The Impact of Macroscopic Epistasis on Long-Term Evolutionary Dynamics[END_REF]. Our observations, deviating around a viable antibody variant, support a model where epistasis plays a lesser role in determining binding affinity than PWM terms, but could still have far ranging impacts. These impacts would include fast initial/slow long term evolution, the blocking of paths to beneficial mutations, while paradoxically an increase in the accessibility of antigen recognizing sequences.

Taken together, our results show the importance of taking into account epistasis when predicting antibody evolution and guiding vaccine design. We note that, while the yeast display system we studied is highly correlated to alternative measures of affinity [START_REF] Gai | Yeast surface display for protein engineering and characterization[END_REF], antibodies could behave differently under alternative environments. Specifically, a soluble antibody may experience different interactions with the environment and be subject to different constraints than an antibody expressed on a cellular surface. Nevertheless, our systematic approach for identifying and quantifying epistasis, which controls for scale and noise, can be used by other investigators to analyze deep mutational scans of protein function in a wide variety of biological contexts.

Main figure legends

Figure 1 -Additive model of binding affinity. (A) 4-4-20 scFv antibody sequence. Six complementarity determining regions (CDR: 1L, 2L, 3L, 1H, 2H, 3H) are particularly important for antibody binding affinity. A library of antibody sequences with mutations in 10 amino-acid regions around the CDR1H and CDR3H domains were expressed using yeast display. (B) Using Tite-Seq, the binding constants K d of all 600 single codon mutants, 1100 random double codon mutants, and 150 random triple codon mutants, were measured. (C) The K d of single mutants for 1H and 3H domains were used to create position weight matrices (PWMs) to predict the affinity of double and triple mutants. K d measurements were restricted to the reliable measurements interval of 10 -9.5 -10 -5 M. WT sequences are marked with purple dots, optimized 4m5.3 mutations are marked in green dots. (D) Comparison between the PWM prediction and the measurement of K d on double and triple mutants. PWM predictions outside of our reliable readout interval (10 -9.5 -10 -5 M) were evaluated at the interval boundaries. The PWMs explained a significant portion of the variance, as quantified by the explained variance R 2 ( p < 10 -61 for CDR1H, p < 10 -48 for CDR3H, F-test for reduction in variance due to PWM). PWMs trained from the binding free energy, F = ln K d c 0

(

) , outperformed PWM trained from K d (Fig. S3) as well as models without boundaries (Fig. S4). Positive Z-scores indicate epistasis increased affinity. The Z score standard deviation was much higher than expected from measurement errors ( Z error ) for CDR1H (1.78, p <10 -16 , Levene's test) and CDR3H

(3.18, p < 10 -48 ), meaning that the discrepancy between the PWM and measurement is mainly due to true epistasis. (C) Standard Z-score deviation for each pair of positions along the sequence. This deviation is higher at pairs of positions mutated in the super-optimized 4m5.3 antibody (green dots) in 3H (p=0.005, Mann-Whitney), but not in 1H (p=0.36). (D-E) A model of biochemical epistatic interactions between polar, nonpolar, acidic, and basic residues was fitted to the data using LASSO regularization and tested by cross-validation, yielding D) 34 beneficial and E) 32 deleterious interaction terms. Line width denotes interaction strength. (F) Number V of amino-acid sequences of the 1H (blue) and 3H (red) regions with dissociation constant below K d , as estimated by the PWM model (dark color) or the epistatic model (light color). Epistasis enlarges the number of variants with good affinity for both 1H and 3H. (G) Mutational flux A (defined as the average number of random mutation events from all possible sequences to cause the dissociation constant to cross K d ), normalized by V, showing that epistasis also increases the accessibility of the region of good binders in sequence space. Differences between the PWM and epistatic models were robust to errors in the estimate of the interaction parameters ( p < 10 -5 , Jackknife analysis).

,

where f = E -1 is the inverse function of E,

, and

E(s) is evaluated similarly to Eq. 3: E s ( ) = 1 / N s ( ) ( ) E ln K d v ( ) / c 0 ( )     v∈S a s ( ) ∑ a ∑ .
To find the transformation E that gives the highest explained variance while avoiding overfitting, we aimed to minimize the following objective function:

O E [ ] = E PWM s ( ) -E s ( )     2 +α dF ′′ E F ( ) 2 ∫ S ∑ , (5)
where the sum in s runs over double and triple mutants, and α is a tunable parameter. Numerically, we parametrize the function E(F) as piecewise linear:

E F ( ) = E i × F i+1 -F ( ) / δF + E i+1 × F -F i ( ) δF for F i ≤ F ≤ F i+1
, where F i are equally spaced grid point along F, δF = F i+1 -F i , and E i the value of E at these points. The smoothing penalty is approximated by a sum over the squared discretized second derivative:

df ′′ E F ( ) 2 ≈ E i+1 + E i-1 -2E i ( ) 2 / δF 3 u ∑ ∫ . We minimize O E [ ] ≈ O E 1 , , E N [
] as a quadratic function of its arguments E i ( ) , while imposing boundary constraints on the PWM prediction and the requirement that E is a increasing function of F (i.e.

E i+1 > E i ), using the python package cvxopt (Andersen et al., 2013).

The hyper-parameter α is evaluated by maximizing the generalized cross-validation of the coefficient of determination

R 2 = 1- E PWM \ S s ( ) -E \ S s ( )     2 s∈S ∑ Var s∈S E \ S s ( )     S ( 6
)
where E \ S and E PWM \ S are learned through optimizing Eq. 5, but after removing from the dataset a subset S of the multiple mutants comprising one tenth of the total. The average is over ten nonoverlapping subsets S.

This method was first tested on simulated data. Each PWM element was drawn from a normal distribution of zero mean and variance 1, and then E PWM s ( ) was computed for each of the antibody sequences .

(10)

The Lasso penalty λ was learned by 10-fold cross-validation, and energy terms found in less than 2 sequences were excluded from the fit.

Posterior values for J terms were calculated using Bayesian Lasso [START_REF] Park | The Bayesian Lasso[END_REF].

The volume and mutational flux were defined as:

V B ( ) = Θ B -K d s ( ) ( ) s ∑ (11) A B ( ) = Θ B -K d s ( ) ( ) 1 19l Θ K d ′ s ( ) -B ( ) ′ s |d s, ′ s ( )=1 ∑ s ∑ , (12)
where τ(x) is the Heaviside function, i.e. τ(x)=1 if x≥0 and 0 otherwise; d(s,s') is the Hamming distance between two sequences; and ℓ=10 is the sequence length. The normalization 19×ℓ corresponds to the number of mutants s' at Hamming distance 1 from s. The sums over s in Eqs. 11-12 have 20 10 elements and are computationally intractable.

To overcome this, we approximated the sum using a mixture of Monte-Carlo and complete enumeration, depending on the distance of s from the wildtype. Calling C d the set of sequences s at Hamming distance d from wildtype, we used:

Figure 2 -

 2 Figure 2 -Quantification of epistasis. A) Epistasis is defined as deviation from the PWM model, which assumes an additive effect of single mutations on the binding free energy F = ln K d c 0 ( ) expressed in units of k B T. Deleterious epistasis occurs when the measured energy exceeds the PWM prediction, beneficial epistasis occurs when the energy is less than the PWM prediction. Sign and reciprocal sign epistasis examples are shown for a beneficial interaction. (B) Distribution of Z-scores, defined as the normalized deviation from the PWM prediction, Z epi = F PWM -F ( ) / σ 2 +σ PWM 2

  s) is a function of s to be summed such as in V or A in Eqs. 11-12, and is a random subset of C d of size min C d , where P is the maximum number of sequences one is willing to sample completely at each d to perform the estimation, and where C d = For small d, when C d ≤ P d , the enumeration is complete, while for large d and C d > P d , the sum is estimated from a uniformly distributed Monte Carlo sample of C d .

Acknowledgements

We would like to thank Yuanzhe Guan and Carlos Talaveira for their suggestions. R.M.A., T.M. and A.M.W. were supported by grant ERCStG n. 306312. JBK was supported by NIH Cancer Center Support Grant 5P30CA045508.

Author Contributions

Contact for reagent and resources sharing

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact Thierry Mora (tmora@lps.ens.fr).

Method details

Values of K d as measured by Tite-Seq for variants of the 4-4-20 fluorescein-binding antibody [START_REF] Adams | Measuring the sequence-affinity landscape of antibodies with massively parallel titration curves[END_REF] can be found at https://github.com/jbkinney/16_titeseq. The scripts used for the analyses presented here are available at https://github.com/rhys-m-adams/epistasis_4_4_20.

Position Weight Matrix

The amino-acid sequence of the 10 amino acid stretches of the CDR1H or CDR3H domains are denoted by s = s 1 , , s 10

(

) . The corresponding 30-long nucleotide sequences are denoted by v. The binding free energy F(s) of an amino-acid variant is obtained as the mean over 3 replicate experiments, and over all its synonymous variants:

where S a s ( ) is the set of measured nucleotide sequences that translate to s in replicate a, and N s

The elements of the PWM are defined as h i q

where s i,q

( ) is the single mutant mutated at position i to residue q, and h i q ( ) = 0 when q is the wildtype residue at position i.

Optimal nonlinear transformation of the free energy

To test whether transforming F through a nonlinear function E(F) before learning the PWM could improve its predictive power, we defined the nonlinear additive model: present in our data. Our simulated "measurement" was defined as a function of a noisy realization of E = E PWM +ε (where ε is some Gaussian noise) in three different ways: linear F=E, exponential F=exp(E), high-frequency F=2E+sin(2E), and logistic F=1/[1+exp (-E)]. ε was drawn from a centered normal distribution with 1/2 the standard deviation of E PWM . F was then truncated to the 200th lowest and 200th highest values, to mimick the boundary cutoff in our measurements. Comparing our original E PWM to our fit Ê shows that our method is able to infer the true PWM model and a smooth nonlinearity from noisy data (Fig. S5).

We then applied the method to the experimental data. The crossvalidation R 2 is represented as a function of the smoothing parameter α in Fig. S6A, and the corresponding optimal function E(F) in Fig. S6B. The comparison between measurement and the PWM model is shown in Fig. S6C.

Z-scores

We used synonymous mutants to estimate our measurement error. The mean free energy of a nucleotide sequence is defined as the mean over replicate measurements:

, and the standard error σ(v) is defined accordingly as the pooled error over replicates. Antibodies with K d or single mutant PWM contributions having median values at the boundary values of 10 -9.5 or 10 -5 were excluded from the analysis since these values artificially cluster at the boundary, leading to underestimates of error.

The error Z-score was calculated between pairs of nucleotide sequences with the same amino acid translation:

Epistatic Z-scores were estimated by calculating the measurement error over both replicates and synonymous variants, as in Eq. 3:

and the pooled standard error for a PWM prediction, calculated as the sum of measurement errors from single mutations:

where σ i q ( ) = σ s i,q ( ) ( ) , and σ i q ( ) = 0 when q is the wildtype residue at i. The epistatic Z-score is defined as:

Null model for sign epistasis

To calculate p-values for sign epistasis, we used the following null model for sets of four Z-scores satisfying

x i has zero mean and variance one, with the constraint

. The distribution with maximum entropy satisfying these requirements is a centered multi-variate Gaussian uniquely defined by its covariance matrix x i 2 = 1 and x i x j = -1 / 3 for i≠j. The p-value for sign epistasis, Z A >1.65 and Z A|B < 1. The threshold of 1.65, in contrast to the previous threshold of 1.64, was determined by applying the Benjamin-Hochberg method for limiting false discovery to 0.05.

Epistatic model

The epistatic terms of the pairwise model were made to depend on the biochemical categories of the interacting residues, , with b(s)=nonpolar for s=AFGILMPVW, b(s)=polar for s=CNQSTY, b(s)=acidic for s=DE, and b(s)=basic for s=HKR. A fifth category was added to correspond to the wildtype residue, so that . The model was trained by minimizing the mean squared error with a regularization penalty over all matrices :

Quantification and Statistical Analysis

All analyses were performed using Python 3.6. Standard statistical tests such as F-test for model comparison (Figure 1C,D), Kolmogorov-Smirnov test for normality (Figure 2B), Levene's test for equal variance (Figure 2B), Mann-Whitney test for equal epistatic contributions by residue position (Figure 2C), and Binomial test for enrichment of sign epistasis (Method Details: Null model for sign epistasis) were performed using the scipy.stats module. Bayesian Lasso (Figure 2D,E), cross-validation (Figure 2D,E), and jackknifing (Figure 2F,G) estimates of probabilities are described in method details, were implemented using Python 3.6, and can be found at https://github.com/rhys-m-adams/epistasis_4_4_20 .