Computer implemented method for predicting a monoclonal antibody specific for a certain antigen of a pathogen

Abstract

The present invention refers to a computer implemented method used to identify in an early phase the affinity of a protein-protein contact on the basis of multiple parameters such as the chemical bond affinity, of the interaction of the residues of contact antibody-protein, of the physical-chemical properties of the interactions. In particular, it is able to predict which clinically available monoclonal antibody is more capable of neutralizing a new emerging variant of pathogen such for example Sars Cov-2.

Description

[0001]SIB BW1280R COMPUTER IMPLEMENTED METHOD FOR PREDICTING A MONOCLONAL ANTIBODY SPECIFIC FOR A CERTAIN ANTIGEN OF A PATHOGEN TECHNICAL FIELD OF THE INVENTION The present invention refers to a computer implemented method used to identify in an early phase the affinity of a protein-protein contact on the basis of multiple parameters such as the chemical bond affinity, of the interaction of the residues of contact antibody- protein, of the physical-chemical properties of the interactions. In particular, it is able to predict which clinically available monoclonal antibody is more capable of neutralizing a new emerging variant of pathogen such for example Sars Cov-2. BACKGROUND OF THE INVENTION Exploiting structural bioinformatics tools to analyse both epitopes and protein-protein interactions, allows to understand the effectiveness of antibodies on new fast-evolving variants of the virus with rapidity, reliability as well as innovation and affordability. The mapping of the mAbs on Receptor Binding Domain (RBD) to identify epitopes-based classes is a first step to mitigate the harmful impact of the emergence of antibody- resistant Sars-Cov-2 (SC2) variants, allowing the right antibody to neutralize the variant, targeting the specific epitope. However to date, the only methods for evaluating the efficacy of monoclonal antibodies against Sars-Cov-2 were in vitro assays which required experimental times and expenditure of reagents, in some cases of considerable cost. It followed that the emergence of a new variant was necessary to test the antibodies experimentally having to wait several weeks. On the other hand, designing new antibodies effective against new variants required numerous experiments, lengthening the times by months before obtaining an advantageous antibody. There is therefore a need of providing a tool able to predict in an early phase which clinically available mAb is more capable of neutralizing a new emerging variants of a virus, in particular of SC2. SUMMARY OF THE INVENTION SIB BW1280R The Authors of the present invention have developed a computer implemented method able to predict in an early phase the affinity of a protein-protein contact on the basis of multiple parameters such as the chemical bond affinity, of the interaction of the residues of contact antibody-protein, of the physical-chemical properties of the interactions. In particular the method can predict which clinically available monoclonal antibody is more capable of neutralizing a new emerging variant of a certain pathogen, in particular Sars Cov-2. The proposed technical solution to the problems identified above provides a method to complement sequencing experiments that identify new variants, without the need to carry out a molecular dynamics simulation, thus employing considerably shorter times than those of the techniques known in the prior art, with superimposable results. When a new variant arises, the proposed method is able to evaluate the efficacy of the approved antibodies in therapy in a few minutes, providing an immediate response in terms of response to the pandemic. Furthermore, in case a new variant is not identified, all the contact areas of the antibodies, their physicochemical properties and contact sites are available in order to provide information regarding the design of new antibodies. With these antibodies already undergoing clinical trials, the present invention allows to expedite the selection of those most effective against potential threats, such as, for example, the Omicron variant. A first object of the invention is a computer implemented method to predict which available monoclonal antibody against a target protein of a virus binds with the highest affinity a new variant of the target protein of said virus, the method comprising the following steps: i) preparing in silico a first dataset comprising all the available structures of complexes of said antibodies with said target protein and its known variants (virus-antibody complexes); ii) performing more than one molecular dynamic simulations of the complexes of said first dataset; iii) performing a statistical analysis on the molecular dynamic simulations of step (ii) obtaining an affinity contact library; iv) refining the affinity contact library of step iii) in order to provide a training set for the prediction of the interaction affinity of said virus-antibody complexes; SIB BW1280R v) training a neural network with said training set obtained in step iv) to enhance the accuracy of the prediction of the affinity of said virus-antibody complex thereby obtaining a neural network trained on the affinity contact library; vi) running said trained neural network using as input the structure of the complex between the available antibodies and said new variant of the target protein of the virus, thereby identifying the available monoclonal antibody specific for the new variant of said target protein. A second object of the present invention is a computer program comprising instructions which, when the program is executed by a computer, cause the computer to carry out the steps of the method according to the present invention. A third object of the present invention is a computer-readable storage medium comprising instructions which, when executed by a computer, cause the computer to carry out the steps of the methods according to the present invention. Further advantages and / or embodiments of the present invention will be evident from the following detailed description. DETAILED DESCRIPTION OF THE FIGURES Fig.1: Affinity scores in relation to time (ns) for each monoclonal antibody in complex with WT, Beta, Omicron, Delta, Alpha RBDs. a) Plots mAbs belonging class 1. b) Plots mAbs belonging class 2. c) Plots mAbs belonging class 3. d) Plots mAbs belonging class 4. Fig.2: Binding region of mAb most and less effective, mapped on RBDWT b) Binding region of mAb most and less effective, mapped on RBDALPHA. c) Binding region of mAb most and less effective, mapped on RBDBETA. d) Binding region of mAb most and less effective, mapped on RBDDELTA. e) Binding region of mAb most and less effective, mapped on RBDOMICRON. Images were obtained with the Molecular Surface panel implemented in Maestro (Schrodinger Suits), a solid surface style and color scheme based on Atom Color were chosen as viewing settings. Fig.3: Dendrogram of the time series clustering on the molecular dynamics simulations of the Fab-RBD complex. Fig.4: Heatmaps of the median of affinity score by residues for each LyCov555-RBD complex. SIB BW1280R Fig.5: Heatmaps of the median of affinity score by residues for each LyCov016-RBD complex. Fig. 6: Heatmaps of the median of affinity score by residues for each CTP59-RBD complex. Fig.7: The affinity score threshold (dashed vertical line) is obtained by computing the optimal (equal sensitivity and specificity) cut-off between the affinity distribution of stable complex clusters (1 and 3, blue curve) against unstable ones (2 and 4, red curve). Fig. 8: Structural bioinformatic pipeline. From PDB download of selected molecular complex to protein preparation leading to Molecular Dynamic Simulation. In the end, data generated were processed to obtain Affinity Score value as a measure of binding stability between monoclonal antibodies and Sars-cov2 variants. GLOSSARY In the present invention, the term “in silico” refers to computer simulations or modelling of biological processes, performed on a computer or computer program. In the present invention, the term “structure of a virus-antibody complex” refers to the three-dimensional arrangement of atoms and molecules that make up the complex formed between a virus and an antibody. In the present invention, the expression “enhance the accuracy of the prediction” refers to increasing precision at the cost of recall. In the present invention, the term “affinity contact library” refers to a collection of molecular interactions between an antibody and a target protein, e.g binding energy. In a preferred embodiment, it is a collection of protein-protein contact affinity scores, that take into account the binding distance between residues and the temporal evolution of these interactions. In the present invention, the term “statistical analysis” refers to the process of using statistical methods by means of computer software to analyse and interpret data. In the present invention, the term “molecular dynamic simulation” refers to a computer- based simulation of the physical and chemical interactions between molecules over time. The approach of the present invention minimizes computational burden by utilizing, in a preferred embodiment, short 100 ns molecular dynamics (MD) simulations, leveraging clinically validated proteins. SIB BW1280R In the present invention, the term “monoclonal antibody” refers to a type of antibody produced by a single type of immune cell, which recognizes a specific target protein. In the present invention, the term “affinity between an antibody and a target protein” refers to the strength of binding of the antibody to the target proteins well as the stability of their interaction, which can contribute to increased therapeutic efficacy. In the present invention, the term “variant of a target protein” refers to a variant or mutated form of the target protein. In the present invention, variants of a target protein are variants of the Spike protein of Sars-Cov2, i.e. Alpha, Beta, Gamma, Delta, Omicron variants. In the present invention, the term “training set” refers to a dataset used to train a machine learning model. In the present invention, the term “neural network” refers to a type of artificial intelligence algorithm modelled after the structure and function of the human brain. It consists of interconnected nodes, called artificial neurons, that are organized into layers. These nodes process and transmit information by passing signals through the network, which are modified by weights assigned to each connection. The aim of a neural network is to learn patterns in input data and use that knowledge to make predictions or decisions based on new, unknown data. In the present invention, the term “Fab-RBD complexes” refers to complexes formed between a fragment of an antibody (Fab) and the receptor-binding domain (RBD) of a virus. In the present invention, the expression “have a marketing authorization” refers to the process of obtaining approval from a regulatory body to sell a product. Regarding monoclonal antibodies binding the spike protein of the variants of the Spike protein of Sars-Cov2, the following antibodies have a marketing authorization: Bamlanivimab, Etesevimab,Regdanvimab, Tixagevimab, Cligaviman, Sotrovimab In the present invention, the term “PDB (Protein Data Bank)” refers to a large and widely used database that collects and stores information about the three-dimensional structures of proteins, available on the internet site http: / / www.wwpdb.org / . In the present invention, the expression “chain in the crystal structure” refers to a single chain of amino acids in a protein crystal structure. SIB BW1280R In the present invention, the term “HETATM” refers to a type of molecule or atom that is present in the crystal structure of a protein but is not part of the amino acid sequence of the protein itself (HETeroAToMs). In the present invention, the term “RBD mutation” refers to a specific alteration or modification of the genetic sequence in the receptor-binding domain (RBD) of a virus that can affect the virus's ability to bind to and enter host cells, thereby altering its infectivity and pathogenesis. In the present invention, the term “targeted mutation” refers to a specific alteration or modification of the genetic sequence made to a target protein specifically to study its effects. In the present invention, the term "binding scenario" refers to the specific and unique way in which an antibody and its target protein interact and associate with each other. This includes the precise molecular interactions, such as hydrogen bonding and van der Waals forces, as well as the conformations and orientations of the antibody and target protein, which together determine the strength and specificity of the binding interaction. The binding scenario can have a significant impact on the efficacy and performance of the antibody-protein complex, making it a crucial aspect to consider in research and development. In the present invention, the expression “adding missing hydrogen bonds” refers to the process of supplementing the hydrogen bond interactions in a molecular simulation that are not already present in the initial model. This process aims to increase the accuracy of the simulation by taking into account the missing hydrogen bond interactions that occur in reality, thereby improving the simulation's representation of the molecular system being studied. In the present invention, the term “stability of complex” refers to the ability of the complex formed between an antibody and its target protein to persist and maintain its structural integrity over time under different environmental conditions. The stability of the complex is a crucial factor in determining the effectiveness of the antibody in binding to its target protein. In the present invention, the expression “uniforming the nomenclature of the chains” refers to the process of standardizing and consistently naming the chains within a protein crystal structure. This involves assigning a consistent and standardized label to each SIB BW1280R chain in the structure, allowing for easier and more precise communication and analysis of the protein crystal structure. In the present invention, the term “Change Atom Property tool” refers to a software tool used to alter the characteristics or attributes of atoms in a molecular simulation or modelling environment. In the present invention, the term “Maestro Schrodinger Suite” refers to a collection of software programs designed to support molecular modelling and simulation tasks, typically used in biochemistry and computational chemistry. In the present invention, the term “minimizing and equilibrating the system” refers to the process of finding the lowest possible energy state for a molecular simulation and reaching a state of thermal balance in a molecular simulation or model. In the present invention, the term “NPT conditions” refers to a molecular dynamics simulation ensemble in which the number of particles, pressure, and temperature of the system are held constant. In the present invention, the term “adjusting the position of the atoms” refers to a process of altering the spatial arrangements of atoms in a molecular simulation or model. In the present invention, the term “energy local minimum” refers to a state of minimum energy within a molecular simulation or model, which reflects the stability of the system. When the energy of the system reaches a minimum, the system is said to be in an "energy local minimum". This state represents a stable configuration of the system, as any small perturbations to the system will cause an increase in energy and a return to the minimum energy state. In other words, the stability of the system can be inferred from the depth and shape of the energy local minimum. A shallow minimum indicates that the system is only weakly stable and can easily be perturbed to another state, whereas a deep and narrow minimum indicates that the system is very stable and difficult to perturb. In the present invention, the term “steepest descent method” refers to an algorithm used to minimize the potential energy of a system in molecular dynamics simulations. It is an iterative method that starts with an initial configuration of the system and moves it towards a local minimum of the potential energy surface. In the present invention, the term “GROMACS” refers toaopen-source software package for performing molecular dynamics simulations, used to study the dynamics and interactions of large biological systems. SIB BW1280R In the present invention, the term “canonical ensemble (NVT) conditions” refers to a molecular dynamics simulation ensemble in which the number of particles, volume, and temperature of the system are held constant. In the present invention, the term “Noose-Hoover thermostat” refers toa temperature control algorithm commonly used in molecular dynamics simulations to maintain a constant temperature. The thermostat works by adding or removing energy from the system in order to regulate the temperature. It does this by acting as a heat sink or heat source as needed, to ensure that the temperature of the system remains at a desired set point. The Noose-Hoover thermostat is implemented as a feedback loop, continuously monitoring the temperature of the system and adjusting the amount of energy added or removed as necessary. This ensures that the temperature of the system remains constant and stable, which is important for accurately representing the thermodynamics of the molecular system being simulated. In the present invention, the term “cluster analysis” refers toa statistical method used to group together similar data points based on their characteristics or features. In the context of molecular simulations, it is used to identify patterns and relationships within the data. The method uses algorithms to sort the data into different clusters, each cluster containing data points that are similar to one another. In the present invention, the term “common classes of complex” refers to groups of molecular complexes that share similar characteristics or properties, used to categorize and study similar systems in molecular simulations. In the present invention, the term “time series” refers to a set of data points recorded at regular intervals over time, used to study the dynamics of a system in molecular simulations. In the present invention, the term “linear regression modelling” refers to a type of statistical analysis that is used to determine the linear relationship between two or more variables. The method involves fitting a line, called the regression line, to the data points in a scatter plot and using that line to make predictions about the value of one variable based on the value of another. In the context of molecular simulation and experimental data, linear regression model is used to study the relationship between simulation data and experimental data by analyzing the correlation between the two sets of data. The regression line can then be used to predict the values of the experimental data based on the values of the simulation data, in order to validate the accuracy of the simulation data. SIB BW1280R In the present invention, the term “dose response trend analysis” refers to a systematic examination of the relationship between the quantity of a substance (referred to as the "dose") and the resulting changes in a specific response or effect. In molecular simulations, this type of analysis is used to understand how changes in the dose of a compound affects the behavior of the system being modeled. In the present invention, the term “analysis of residues causing low performance” refers to the process of identifying which residues within a RBD are responsible for lower binding stability and / or affinity. In the present invention, the term “R software” refers to open-source software environment for statistical computing and graphics, commonly used for data analysis in molecular simulations, which allows to search contacts with known affinity profile (the contact library) that are similar in certain parameters to those constituting the query Fab- RBD complex. In the present invention, the term “R module” refers to the implementation of the search and the filtering module. In the present invention, the term “side chain biochemical similarity” refers to the comparison of the chemical and biological properties of the side chains of different molecules in order to determine their level of similarity. This comparison includes aspects such as the size, shape, charge, and hydrophobicity of the side chains, as well as any other relevant biochemical properties. These parameters include amino acid side chain polarity, electric charge, and steric hindrance; contact interaction class (I, II, III, IV), defining the interaction location over the RBD surface; known contact mutation, as retrieved from literature. In the present invention, the expression “refining the results by interaction class” refers toa process of improving the accuracy and precision of the results of a molecular simulation or analysis by dividing the interactions into different groups based on their type or class. This process involves analyzing the results for each group of interactions separately and then combining the results to form a more complete and refined picture of the interactions in the complex. In the present invention, the term “contribution to the overall binding stability and / or affinity of the complex” refers to the capability of a single contact to increase the Fab- RBD complex affinity score. SIB BW1280R In the present invention, the term “affinity score” refers to a multifaceted metric which combines binding distance with chemical features and residue properties.In the present invention, the term “do not contribute significantly” refers to the inability of a contact to increase the Fab-RBD complex affinity score. In the present invention, the term “expected spatial localization” refers to the predicted location of an interaction within a molecular simulation or model. This prediction is based on the molecular properties of the elements involved in the interaction, as well as their relative arrangement in the simulation or model. In the present invention, the term “matching with mutations in known variants” refers to the exact correspondence between an amino acid mutation found in an input contact and one already present in literature. In the present invention, the term “prioritizing specific amino acid residues” refers tothe process of assigning a priority to specific amino acids within a protein based on their contribution to the stability and / or affinity of the Fab-RBD interaction. In the present invention, the term “refining by mutation hotspots” refers to the process of refining the results of an analysis by focusing on regions or areas within the molecular structure that have a high concentration of mutations. By focusing on these hotspots, the analysis aims to identify if the mutations are having a significant impact on the stability and / or affinity of the protein-protein interaction. In the present invention, the term “odds ratio of belonging to an interaction localization class” refers tothe ratio that is calculated to determine the odds that an interaction between two molecules belongs to a particular localization class, compared to the odds of the interaction not belonging to that class. This ratio helps to identify which localization class an interaction is most likely to belong to. In the present invention, the term “most likely to affect the binding stability” refers tothe molecules or factors that are believed to have the greatest impact on the stability of the complex. These molecules or factors can interact with the residues in the complex and cause changes in the energy state of the system, which can lead to changes in the binding stability of the complex. These changes can be due to changes in the conformation of the residues, changes in the interactions between the residues and other molecules, or changes in the electrostatic potential of the residues. In the present invention, the term “training-prediction iteration” refers to the process of repeatedly training a machine learning model and making predictions based on that SIB BW1280R training, in order to improve the accuracy of the predictions. The process typically involves adjusting the parameters of the model, such as the weights of the neurons, based on the results of the predictions, and then using the adjusted model to make new predictions. This process is repeated multiple times, until the model reaches a state where its predictions are accurate enough for the desired application. In the present invention, the term “stratifying” refers to the process of dividing data into smaller groups based on certain characteristics, in order to analyze the data in a deeper way. This can involve grouping the data based on attributes such as the type of molecule, the location of the residues in the complex, or the binding affinity of the residues. In the present invention, the term “validation set” refers to a set of data that is used to validate a machine learning model or method, by comparing its results with the expected outcomes. This set of data is typically independent from the data used to train the model and is used to estimate the performance of the model when applied to new, unseen data. By comparing the results of the model with the expected outcomes for the validation set, it is possible to determine whether the model is overfitting to the training data, or whether it is generalizing well to new data. In the present invention, the term “Python3 module” refers to a module written in the Python3 programming language. Python3 is a high-level, interpreted language that is widely used for scientific computing, data analysis, and machine learning. The module can be a collection of functions or classes that perform specific tasks, and can be imported and used in other Python programs. In the present invention, the term “LSTM module” refers toa module based on a type of artificial neural network called a Long Short-Term Memory (LSTM) network. LSTMs are a type of recurrent neural network that are designed to capture long-term dependencies in sequential data. They are particularly well-suited for applications in natural language processing, speech recognition, and time series prediction. In the present invention, the term “leave-one-out cross-validation strategy” refers toa technique for validating a machine learning model, where the model is trained on all but one data point, and the prediction for that one data point is used as the validation. This process is repeated for each data point, such that the model is trained on all but one data point at a time. In the present invention, the term “LSTM recurrent neural network” refers toa type of artificial neural network that is based on the Long Short-Term Memory (LSTM) SIB BW1280R architecture and is used for modelling time-series data, such as speech and text, where there are dependencies between inputs and outputs over time. The LSTM recurrent neural network uses feedback connections and memory cells to store and access information over long sequences of data, allowing it to make predictions based on a long history of inputs. In the present invention, this technology is used to analyze and predict the binding stability and / or affinity of complexes in molecular simulations. DETAILED DESCRIPTION The present invention relates to a computer implemented method to predict which available monoclonal antibody against a target protein of a virus binds with the highest affinity a new variant of the target protein of said virus, the method comprising the following steps: i) preparing in silico a first dataset comprising all the available structures of the complexes of said antibodies with said target protein and its known variants; ii) performing more than one molecular dynamic simulations of the complexes of said first dataset; iii) performing a statistical analysis on the molecular dynamic simulations of step ii) obtaining an affinity contact library; iv) refining the affinity contact library of step iii) in order to provide a training set for the prediction of the interaction affinity of said virus-antibody complex; v) training a neural network with said training set obtained in step iv) to enhance the accuracy of the prediction of the affinity of said antigen-antibody complex thereby obtaining a neural network trained on the affinity contact library; vi) running said trained neural network using as input the complex structure between the known antibodies and said new variant of the target protein of the virus, thereby identifying the available monoclonal antibody specific for the new variant of said target protein. The expression “specific for the new variant” means that the identified monoclonal antibody is the monoclonal antibody with the highest affinity value to the new variant of the target protein of the virus. In an embodiment, said virus is a Coronavirus. In a preferred embodiment, said virus is a Sars-Cov2. In another embodiment, said target protein is the Receptor Binding Domain (RBD) of the Sars-Cov2. SIB BW1280R In an embodiment of the present invention, said complexes are Fab-RBD complexes (Fragment antigen-binding – Receptor Binding Domain). In a preferred embodiment, said antibodies have a marketing authorization for the treatment of a disease caused by said virus. In particular, said antibodies are monoclonal antibodies approved for the therapy against Sars-Cov-2. In another embodiment, said first dataset comprises the complex of each available monoclonal antibody with each available variant of the target protein, in particular the Fab-RBD complex. In a more preferred embodiment, said Fab-RBD complexes are complexes between the Fabs of all available monoclonal antibodies against the RBD of all the variants of the Sars-Cov2. In an embodiment, the affinity contact library is a collection of protein-protein contact affinity scores. In a preferred embodiment, said contact affinity scores take into account the binding distance between residues and the temporal evolution of these interactions. In a preferred embodiment, the switch function is applied to the binding distance (r) to calculate the "Affinity Score." This score reflects the strength of interaction, with lower distances (smaller r values) corresponding to stronger affinity or interaction, using different variations of the switch function, the affinity score id adjusted for different type of interactions. The switching function helps smooth the transition of the affinity score between 0 and 1, providing a continuous measure to characterize weighted interactions in the molecular dynamic simulations. 1. Binding Distance: In an embodiment, the Affinity Score relies on the binding distance (denoted as r), which represents the spatial proximity between interacting molecular entities. This distance is a critical parameter, reflecting the physical closeness of atoms or residues involved in the interaction. 2. Chemical Features: In an embodiment, the Affinity Score is designed to capture the nuanced nature of molecular interactions by incorporating various chemical features: a. Type of Binding: SIB BW1280R In an embodiment, the nature of the interaction, such as hydrophobic, pi-stacking, salt bridges, etc., is considered. This allows for differentiation between diverse interaction types with distinct chemical characteristics. b. Residue Properties: In an embodiment, chemical properties of interacting residues play a crucial role. Residues involved in interactions are characterized by their chemical composition, polarity, and other relevant properties as charge. This information enriches the Affinity Score by accounting for the specificities of residue-residue interactions. 3. Weighting Scheme: In an embodiment, the present invention employs a weighting scheme to emphasize the significance of certain interactions or chemical features in the Affinity Score calculation. This allows for a more nuanced interpretation of the molecular interactions, giving greater importance to specific binding events or residue characteristics. 4. Switching Function: In an embodiment, the Affinity Score is computed using a switching function, which ensures a smooth transition between 0 and 1 based on the binding distance. This function, customized with parameters like r_0, a, and b, imparts flexibility to the Affinity Score, allowing it to adapt to the unique characteristics of different molecular interactions. 5. Continuous Measure of Affinity: In an embodiment, the resulting Affinity Score is a continuous measure, providing a quantitative representation of the strength of interactions. Lower Affinity Scores correspond to stronger interactions, while higher scores indicate weaker or negligible interactions. 6. Customization for Interaction Types: In an embodiment, recognizing the diversity of molecular interactions, the present invention tailors the Affinity Score calculation for specific interaction types. In a preferred embodiment, pi-stacking interactions may have different parameter settings (r_0, a, and b) compared to hydrophobic interactions, ensuring adaptability and accuracy. SIB BW1280R In a particular embodiment of the present invention, said step (i) of preparation complexes of said antibodies with said target protein and its known variants comprises the steps of: (i-a) downloading the crystal structures of complexes of said Fab-RBD complexes from the database such as PDB (Protein Data Bank); (i-b) reassigning chain names by identifying and renaming the chains in the crystal structure to correspond to the Fab and RBD components of the complex; (i-c) removing HETATM residues by identifying and removing any non-standard amino acids or other molecules present in the crystal structure that do not correspond to the Fab or RBD; (i-d) implementing RBD mutations by introducing targeted mutations to the RBD component of the complex to simulate various binding scenarios using the software known to the skilled person; (i-e) adding and optimizing hydrogen bonds between residues by identifying and adding missing hydrogen bonds between the Fab and RBD residues to optimize the stability of the complex. In the method of the present invention, said step (i-b) consists in uniforming the nomenclature of the chains in each Fab-RBD complex, and is a crucial step in ensuring the proper management of the simulation code. The Change Atom Property tool, available in the Maestro Schrodinger Suite, can be utilized to standardize and homogenize the data. Furthermore, said step (i-c) of removing HETATM residues is useful for the simplification of the system, to improve the accuracy (the interactions between these atoms and the main protein structures may introduce artifacts into the simulation) and to focus on specific interaction of interest, such as those between Fab and RBD. According to the present invention, said step ©was performed by in-house script in Python3. In one embodiment of the method of the present invention, said step (ii) of performing more than one molecular dynamic simulations of said complexes of said first dataset comprises the steps of: (ii-a) setting up simulation parameters; SIB BW1280R (ii-b) minimizing and equilibrating the systems by adjusting the position of the atoms reducing the initial energy of the system and bringing it to a stable state before running the simulation; (ii-c) running simulations for 100 ns under NPT conditions by running the simulation for a specific duration and under a pressure of 1 atm, and a temperature from 20°C to 40°C, preferably 30°C, obtain a representative sample of the system's dynamics as close as possible to the real system, assuming that the experimental teston the real system was done at room temperature. In an embodiment of the present invention, said simulation parameters for the simulation of step (ii-a) are the following: dt, nsteps, cutoff-scheme, nstlist, rlist, coulombtype, rcoulomb, vdwtype, vdw-modifier, rvdw_switch, rvdw, tcoupl, tc_grps, tau_t, ref_t, pcoupl, pcoupletype, tau_p, compressibility, ref_p, constraints, constraint_algorithm, continuation, nstcomm, comm_mode, comm_grps, refcoord_scaling. In a preferred embodiment of the present invention, said simulation parameters for the simulation of step (ii-a) are: “dt = 0.002 ps nsteps = 50000000 cutoff-scheme = Verlet nstlist = 20 rlist = 1.2 nm coulombtype = pme rcoulomb = 1.2 nm vdwtype = Cut-off vdw-modifier = Force-switch rvdw_switch = 1.0 nm rvdw = 1.2 nm tcoupl = Nose-Hoover tc_grps = SYSTEM tau_t = 1.0 ps ref_t = 303.15 K pcoupl = Parrinello-Rahman pcoupltype = isotropic tau_p = 5.0 ps SIB BW1280R compressibility = 4.5e-5 bar-1ref_p = 1.0 bar constraints = h-bonds constraint_algorithm = LINCS continuation = yes nstcomm = 100 comm_mode = linear comm_grps = SYSTEM refcoord_scaling = com” In an embodiment of the present invention, said step (ii-b) is carried out by adjusting the position of the atoms until the energy reaches a local minimum, through the steepest descent method, available in GROMACS. Equilibration is performed to bring the system to thermal equilibrium. The equilibration procedure in GROMACS can be executed in canonical ensemble (NVT) conditions for 125 ps, using Nose-Hoover thermostat. Furthermore, in one embodiment of the present invention, said step (iii) of performing statistical analysis on the molecular dynamic simulations of step ii) comprises the steps of: (iii-a) cluster analysis by grouping together similar conformations of the complex according to their stability and / or binding affinityin order to detect the presence of common classes of complex, based on their time series stability and / or binding affinity; (iii-b) linear regression modelling and dose response trend analysis by identifying the relationship between changes in the RBD and corresponding changes in binding stability and / or affinity in order to characterize and better interpret the common classes of complex detected in previous analysis step. Concerning that, a number of parameters were estimated by statistical modelling to describe all of time series generated by in silico simulations. Particularly, by linear regression modelling were estimated the intercept (i.e., the affinity score value at baseline) and the slope (i.e., the affinity score expected variation per an additional nanosecond) parameters. Regarding the dose response trend analysis, the estimated parameters was four, labelled as “b”, “c”, “d” and “e”: “b” was the slope (of the dose response modelling), “c” and “d” the lower and upper limits of the affinity score expected from the model and the “e” was the nanosecond corresponding to the expected half-way affinity score between the upper and lower limit (also called ED50). SIB BW1280R (iii-c) analysis of residues causing low performance clusters by identifying which residues in the RBD are responsible for lower binding stability and / or affinity in order to identify the residues causing low performance classes of complex (clusters),by measuring descriptive statistics as the median value were computed on the residues affinity returned by complex. In an embodiment of the present invention, said steps (iii-a), (iii-b) and / or (iii-c)are carried out by a R software. In another embodiment of the invention, said step v) of refining said affinity contact library based on the statistical analysis results comprises the steps of: (iv-a) refining the affinity contact library by interaction class by grouping together similar interactions between the Fab and RBD residues according to their contribution to the overall binding stability and / or affinity of the complex, summing up the contribution of side chain biochemical similarity, expected spatial localization, and matching with mutations in known variants by means of a similarity search module, thereby obtaining a collection of contacts (referred to as “filtered contact library”) that are exactly matching or similar to the input ones. This allows the neural network to learn from known time series with similar features, increasing prediction specificity for a given Fab-RBD pair, reducing prediction error, and eliminating those interactions that do not contribute significantly to the stability and / or affinity of the complex; (iv-b) refining said first refined affinity contact library by mutation hotspots by identifying and prioritizing specific amino acid residues in the RBDthat are most likely to affect the binding stability and / or affinity when mutated, by evaluating the contact priority on the base of the odds ratio of belonging to an interaction localization class by means of a filtering module which provides a quality value of the similarity search on a scale from 0 (totally unreliable) to 15 (exact match in the contact library). This is done to further refine the training set, since partial similarity, shift in localization, or matches with known mutants, may reveal RBD residues that are most likely to affect the binding stability and / or affinity when mutated; (iv-c) stratifying the affinity contact library into a training set and a validation set by dividing the library into two distinct sets for the purpose of training and testing a neural network for predicting the stability of the antigen-antibody complex. This is achieved using a leave-one-out cross-validation strategy where, for each training-prediction iteration, a set of contacts corresponding to a single Fab-RBD complex is excluded from the training set and used as validation set. SIB BW1280R The aim of the validation set is to evaluate the accuracy of the neural network using a set of known contacts, biochemically similar to those constituting the query Fab-RBD complex, in predicting the true Fab-RBD affinity values over a time series of 100 nanoseconds. Model training, prediction of the affinity time series, and cross-validation steps are performed by a LSTM module (Python3 module). In a preferred embodiment the method of the present invention uses a two-step approach based on a long short-term memory (LSTM) recurrent neural network. In the first step, affinity scores from single RBD- receptor interactions are filtered on the base of side chain biochemical similarity, spatial localization of the contact, and presence of possible mutation hotspots. This will select known affinity profiles that are predicted to be similar to the unknown interaction. In the second step, the LSTM predicts the affinity profile of the unknown contact using the training set. If the median predicted affinity score is above or equal to 0.88, the interaction is predicted to be stable within 100 ns, it is unstable otherwise. The affinity threshold is derived from cluster analysis of the affinity profiles from 4 variants (alpha, beta, delta, omicron) and the wild type. In an embodiment, said steps (iv-a) and (iv-b) are carried out by a R module. In an embodiment of the present invention, said neural network is a long short-term memory (LSTM) recurrent neural network. In one embodiment of the present invention, said step vi) of identification of the available monoclonal antibody specific for a new variant of a target protein is carried out with the calculation of the median affinity value of each time series (100 nanoseconds) once every Fab-RBD complex has been predicted. Antibodies can then be ranked by decreasing median affinity. In any case, if the median affinity value drops below 0.88, the Fab-RBD complex is labelled as unstable. In a preferred embodiment of the present invention, the method of the present invention comprises the following steps: i) preparing in silico a first dataset comprising all the available structures of complexes of said antibodies with said target protein and its known variants (virus-antibody complexes), in particular Fab-RBD complexes; (i-a) downloading crystal structures of complexes of said Fab-RBD complexes from a database, such as PDB (Protein Data Bank); SIB BW1280R (i-b) reassigning chain names by identifying and renaming the chains in the crystal structure to correspond to the Fab and RBD components of the complex; (i-c) removing HETATM residues by identifying and removing any non-standard amino acids or other molecules present in the crystal structure that do not correspond to the Fab or RBD; (i-d) implementing RBD mutations by introducing targeted mutations to the RBD component of the complex to simulate various binding scenarios; (i-e) adding and optimizing hydrogen bonds between residues by identifying and adding missing hydrogen bonds between the Fab and RBD residues to optimize the stability of the complex. ii) performing more than one molecular dynamic simulations of theFab-RBD complexes of said first dataset; (ii-a) setting up simulation parameters; (ii-b) minimizing the systems preferably through a steepest descent algorithm and equilibrating the systems by adjusting the position of the atoms thereby reducing the initial energy of the system and bringing it to a stable state in canonical ensemble (NVT) conditions for 125 ps, preferably using Nose-Hoover thermostat, before running the simulation; (ii-c) running simulations for 100 ns under NPT conditions30°C, preferably using Nose-Hoover thermostat, with a T-coupling constant of 1 ps, and a Parrinello-Raman barostat at 1 atm, to obtain a representative sample of the system's dynamics. iii) performing a statistical analysis on the molecular dynamic simulations of step (ii) obtaining an affinity contact library; (iii-a) cluster analysis by grouping together similar conformations of the complex according to their stability and / or binding affinity in order to detect the presence of common classes of complex, based on their time series stability and / or binding affinity; (iii-b) linear regression modelling and dose response trend analysis by identifying the relationship between changes in the RBD and corresponding changes in binding stability and / or affinity in order to characterize and better interpret the common classes of complex detected in previous analysis step by measuring, for the linear regression modelling the affinity score value at baseline, the affinity score expected variation per an additional nanosecond, and for the dose response trend analysis the slope (of the dose SIB BW1280R response modelling), the lower and upper limits of the affinity score expected from the model and the nanosecond corresponding to the expected half-way affinity score between the upper and lower limit (also called ED50). (iii-c) analysis of residues causing low performance clusters by identifying which residues in the RBD are responsible for lower binding stability and / or affinity in order to identify the residues causing low performance classes of complex (clusters). Concerning that, the median of the affinity score values of the mean trajectory (across the three independent replicates) were computed. The median was used because of its robustness to the outliers and considering the decreasing trend (approximately) of the trajectory, whose starting point at baseline was affinity score=1. iv) refining the affinity contact library of step iii) in order to provide a training set for the prediction of the interaction affinity of said Fab-RBD complexes; (iv-a) refining the affinity contact library by interaction class by grouping together similar interactions between the Fab and RBD residues according to their contribution to the overall binding stability and / or affinity of the complex, considering at least one of side chain biochemical similarity, expected spatial localization, and matching with mutations in known variants by means of a similarity search module and eliminating those interactions that do not contribute significantly to the stability and / or affinity of the complex (similarity search module)thereby obtaining a first refined training set; (iv-b) refining said first refined affinity contact library obtained in (v-a) by mutation hotspots by identifying and prioritizing specific amino acid residues in the RBD that are most likely to affect the binding stability and / or affinity when mutated, by evaluating the contact priority on the base of the odds ratio of belonging to an interaction localization class by means of a filtering module which provides a penalty value of the similarity search on a scale from 0 (exact match in the contact library) to 15 (totally unreliable) (filtering module), in order to further refine the affinity contact library; (iv-c) stratifying the further refined affinity contact library into a training set and a validation set by dividing the library into two distinct sets for the purpose in order to of training and testing a neural network for predicting the stability of the antigen Fab– antibody RBD complex by means of a leave-one-out-cross-validation strategy where for each training-prediction iteration, a set of contacts corresponding to a single Fab-RBD complex is excluded from the training set and used as validation set. SIB BW1280R v) training a neural network with said training set obtained in step iv) to enhance the accuracy of the prediction of the affinity of said Fab-RBD complex thereby obtaining a neural network trained on the affinity contact library; vi) running said trained neural network using as input the structure of the complex between the available antibodies and said new variant of the target protein of the virus, thereby identifying the available monoclonal antibody specific for the new variant of said target protein. In one embodiment of the present invention, said step of refining said affinity contact library is performed using a R module. The aim of the R module is to search contacts with known affinity profile (the contact library) that are biochemically similar to those constituting the query Fab-RBD complex. Once similar contacts have been identified, the contact library is filtered to include only these contacts, constituting the training set. This filtering step is done to improve affinity prediction specificity. In another embodiment, said training of said neural network is performed using a Python3 module. This module learns the affinity trend from known contacts to predict contact affinity in complexes with unknown affinity profile. Model training and affinity time series prediction are done with long short-term recurrent neural networks, to achieve accurate in silico predictive performances, in absence of wet-lab and / or computational resources, as a tool for fast evaluation of virus-antibody complex stability. Another object of the present invention is a computer program comprising instructions which, when the program is executed by a computer, cause the computer to carry out the steps of the method according to any of the embodiments of the present invention. A further object of the present invention is a computer-readable storage medium comprising instructions which, when executed by a computer, cause the computer to carry out the steps to any of the embodiments of the present invention. Throughout this specification and the claims, the term “comprising” may be replaced by the term “consisting of”. Some examples are given below which have the purpose of better illustrating the methodologies disclosed in the present description, such examples are in no way to be considered a limitation of the previous description and of the subsequent claims. EXAMPLES RESULTS SIB BW1280R 1 - Evaulation of Molecular Dynamics Simulation The setting of the computational model provides for the selection of starting point on the full coverage of the RBD surface and the mAbs approved for clinical use. Therefore, 7 crystallographic structures of the RBD-Fab complexes, obtained with X-Ray Diffraction method, were selected from the Protein Data Bank (PDB) (Table 1). The selection was based on the resolution of the crystal (≤ 3 Å) and on the experimental evidence of affinity with the RBD portion of the non-mutated Spike protein (wild type - WT), in order to evaluate how the interaction of antibodies changes with the different variants of the protein (Alpha, Beta, Delta, Omicron) (Table 2). Furthermore, to have a more explanatory vision of the reality, the 7 selected monoclonal antibodies fall into differentiated epitopes that homogeneously cover the RBD domain, providing at least one representative for each of the four identified class on the basis of RBD binding. Table 1: Main features of monoclonal antibodies selected for analysis, starting from corresponding PDB code, through crystal’s characteristics, up to clinical details. Note. G.U.: official gazette. AIFA: Agenzia Italiana del Farmaco. EMA: European Medicines Agency. Table 2: Details of Sars-Cov-2 variants selected for analysis Thirty-five complexes of mAb and RBD variants were obtained and processed as described in the Methods section. MD simulations’ results showed a different approach SIB BW1280R in terms of RBD neutralizations depending on the mAb binding class, that could provide a general trend for the MDs: the figure 1 shows the plots of the MD stratified by antibody (figure 1). Affinity scores time series achieved from MD were similar to neutralization titers of in vitro experiments emerged in literature, only AZD8895-RBDDELTA complex showed an opposite trend. Our results pointed out the interaction between AZD8895 (Tixagevimab) and RDBWT as the most effective. By the contrast, Ly-Cov-555 (Bamlanivimab) complex with RBDBETA, accounted for the worst affinity. On the basis of the RBD variants, epitopes that were advantageous for effective mAb binding, as well as epitopes that were linked to detrimental outcomes are shown in figure 2. 2 – Statistical analysis 2.1 – Descriptive statistics Notably, all the complex time series were normalized (range affinity score 0-1) and comparable because starting from baseline affinity score equal to 1 (max value) and having a not increasing trajectory. Preliminary, table 1S in supplementary shows the descriptive statistics for each trajectory in terms of quartiles and range. 2.2 – Clustering for time series Mean 0-100 time series of the molecular dynamics for each Fab-RBD complex are shown in supplementary (figure 8S). Time series clustering provided the clusters of Fab- RBD complexes in relation to the trends of affinity score across 100 nanoseconds. Dendrogram of the clustering is shown in figure 3 and presented four distinct clusters cutting at distance approximatively equal to 10-14: the first (green cluster from here on) contained 15 Fab-RBD complexes and it was distinguished by a trend approximately constant with values of affinity very high, near to one. In the second one (blue cluster, 11 complexes) the mean trajectories had high values of affinity and denoting a slight decreasing till values around 0.8. The third one (red cluster, 6 complexes) encompassed time series more decreasing than those contained in green and blue clusters. These decreasing trajectories stabilized to affinity scores near to 0.7. The last cluster (i.e., black cluster) contained the worst mean trajectories, i.e., the trends most decreasing that presented the worst performances of affinity between variant and antibody. Of note, this cluster involved the following 3 Fab-RBD complexes: Beta-Lycov555, Omicron- Lycov016 and Omicron-CTP59. SIB BW1280R Notably, the 4-cluster pattern provided the solution more consistent and interpretable. 2.3 – Linear and dose response trend analysis In order to quantitatively characterize the time series clusters, a number of statistical methods were applied. At first, linear analysis showed that the lower (or bigger in absolute value) estimated slopes (b1) belonged to the Fab-RBD complexes in the red and black clusters. On the other hand, the green and blue clusters show the bigger slopes (b1) values. Thus, 4-parameter dose response modelling were applied and estimates of the four parameters (b, c, d, e) (see methods) were obtained as reported in table 1 in relation to the clusters. Of note, the parameters c, d and e were the most important in order to compare the estimated trend. Concerning that, the Fab-RBD complexes belonging to the black cluster showed the estimated values of c parameters lower (and significant) than others. Next, change point analysis and broken line regression were performed to extract informative hot spots, in terms of nanosecond and affinity score, where the trajectory changed slope in relevant way. Regarding that, optimal change points are shown in supplementary table 1S, and among these, the strongest change points was also provided. Notably, the number of optimal change points were a proxy of the regularity of the time series. Regarding that, the Fab-RBD complexes belonging to the black cluster had the strongest change points and the x-parameters of the broken line regression lower (and significant) than others. Notably it is worth to point out that the dose response model of BETA-LyCov016 and WT-CPT59 complexes fits in not optimal manner as reported from the estimates on c parameter. Moreover, table 2 shows the results of the joint 4-parameter dose response modelling of the antibodies stratified by variant, vice versa the table 3 shows the results of the modelling of the variants stratified by antibody. The supplemental figures 2 and 3 report the corresponding time series plots and 4-parameter dose response curves of the stratified analysis by both variant and antibody. (BETA-LyCov016, WT-CPT59) 2.4 – Analysis of the residues Concerning the Fab-RBD complexes, focusing on the bad (black) cluster, all residues of the Beta-Lycov555 (Figure 4) have median values under 0.88; particularly, the residues heavy chain R50 (0.2366) and light chain Y92 (0.4633) and R96 (0.2533) reported values under 0.5. Similarly, the Omicron-Lycov016 (Figure 5) showed low values on the R97 SIB BW1280R (0.37) of the heavy chain and Y32 (0.44) and Y92 (0.3466) of the light chain. Finally, Omicron-CTP59 (figure 6) showed low values on heavy chain S32 (0.4267), D57 (0.4), and R109 (0.47), and light chain Y50 (0.4366). Notably, stratifying for antibody, there were some residues that provided low median affinity scores across variant, both overall and relatively to specific Fab-RBD complexes. For example, LyCov016 residues presented low median values on R97 of the heavy chain and Y92 of the light chain. Similarly, S32 residue of the light chain showed low values relatively to the AZD1061 antibody. Moreover, the residue Y54 of the heavy chain of SOTROVIMAB and Q27 of the light one of EY64 reported median affinity lower. 3 – Complex stability prediction Considering time series belonging to clusters 1 and 3 as stable dynamics, and clusters 2 and 4 as unstable, we obtained an affinity score threshold of 0.88 (Figure 7), with an AUC of 94% (95% CI: 90-98%). If the median affinity score of a time series is above or equal 0.88, the complex is stable; it is unstable otherwise. According to this criterion, LSTM predictions are always in accordance with molecular dynamics simulations, except for the delta-AZD8895, with an overall accuracy of 93% of correct predictions. In the case of delta-AZD8895, the molecular dynamics simulation defines an unstable time series, while the LSTM predicted time series is stable, in accordance with literature. METHODS For each of Fab-RBD complexes (table 1, Results) the following workflow was performed to generate data for the construction of computational model (Fig. 8). All of the steps were performed from the command line, wrapped in a sh file to automatically carry out the analysis (Protocol). 1 – In silico preparation of Fab-RBD complexes Firstly, the selected crystal structures of Fab-RBDWT complex were downloaded from PDB. Change Atom Property tool implemented in Maestro, Schrodinger Suite, was used to reassign chain name in order to uniform the data (Supplementary Information). Protein preparation protein was performed starting from PDB cleaning, through an HETATM removal in-house script (NOHET.py). Subsequently, on the RBD portion of each complex, the RBD mutations [https: / / www.who.int / activities / tracking-SARS-CoV-2-variants], reported in Table 2.2 (Results section), were implemented through Pymol's python Mutagenesis Wizard library. The molecular complexes were prepared by Protein Preparation Wizard tool [Schrodinger suite], through which the atoms of hydrogen were SIB BW1280R added and their positions were modified by optimizing the network of hydrogen bonds between residues based on their calculated pKa. All parameters were set by default. 2 – Molecular dynamics simulation 2.1 – Parameters set-up Molecular dynamics simulations trajectories were carried out using GROMACS 2020.4 package with OPLS-AA force-field and water solvent, choosing TIP3 model for water molecules. The protein–protein systems were solvated in a rhombic dodecahedron (xy- square) water box under periodic boundary conditions, setting the image distance to 1 nm. The total charge of the system was neutralized by randomly substituting water molecules with Na+ and Cl– ions, to obtain neutrality with 0.15 M salt concentration. System was minimized through a steepest descent algorithm and equilibrated in canonical ensemble (NVT) conditions for 125 ps, using Nose-Hoover thermostat. Then, the molecular dynamics simulation was run for 100 ns under NPT (Number of particles, Pressure and Temperature) conditions at 303.15 K, using Nose-Hoover thermostat, with a T-coupling constant of 1 ps, and a Parrinello-Raman barostat at 1 atm. Particle-Mesh Ewald method was used for the calculation of electrostatic interactions of the systems. A normal cut-off of 1.2 nm was implied for Van der Waals interaction. All bonds with H- atoms were constrained using the LINCS algorithm. The time step employed was 2 fs, and the coordinates of trajectory were saved every 10 ps for analysis. Three independent molecular dynamic simulations were performed for each Fab-RBD complex examined. 3 – Statistical analysis The statistical analysis explored the pattern of the molecular dynamics simulations in relation to the affinity score across nanosecond for each RBD-Fab complex. The simulations generated three independent replicates of Fab-RBD simulation trajectories that were averaged by nanosecond (0-100) and returned the mean trajectories, one by Fab-RBD complex, on which the statistical analysis was performed. The analysis encompassed 2 steps: i i) a cluster analysis for time series; ii ii) a linear regression modelling and dose response trend analysis for characterizing the achieved clusters. SIB BW1280R In addition, a third step (iii) regarding to the analysis on the residues was carried out in order to investigate the potential RBDs causing low performances clusters of Fab-RBD complexes. All the analysis was performed by R software. 3.1 – Cluster analysis Cluster analysis for time series was performed to detect the complexes (Fab-RBD complexes trends) that could be clustered. Preliminarily, all the 0-100 time series (one by Fab-RBD complex) of affinity score were include in a data matrix. Hence, a dissimilarity matrix was computed by the Dynamic Time Warping (DTWARP) method, in order to perform the hierarchical cluster analysis by using the complete-linkage agglomeration method. In order to select the optimal number of clusters (k), the dendrogram was evaluated to define the suitable distance where to cut it to obtain a clear clustering partition. This step was performed by R packages Tsclust and ggdendro. ii) Secondly, a linear regression modelling and dose-response trend analysis was applied to characterize the Fab-RBD complex clusters previously obtained. Preliminary, the linear regression modelling was fitted to get the intercept (b0) and slope (b1) parameters. The intercept was interpreted as expected affinity score at baseline, while slope was expected affinity score variation per nanosecond. Regarding the dose-response trend analysis, the nanoseconds (i.e., the doses) were evaluated on the x-axis, whereas the affinity score (trend) was plot on the y-axis. Specifically, we fitted dose response models by a 4-parameter log-logistic function, where the parameters are the slope (b), the lower (c) and upper (d) limits of the affinity score expected from the model and the ED50 (e), i.e., the nanosecond corresponding to the expected half-way affinity score between the upper and lower limit, Notably, all the parameters of the statistical modelling were tested by two tailed T-tests, so P value less than 0.05 revealed values different from zero in significant way. This part of analysis was performed by R packages drc and devtools. In this step, a change point analysis was also performed to detect the structural changes of affinity score across nanoseconds for each Fab-RBD complex. In detail, the change point analysis detected the structural changes of affinity score across nanoseconds for each Fab-RBD complex. In this phase, two sub-analyses were performed, by considering frameworks of breakpoint estimation and F (Chow) statistical testing: in the first one, multiple optimal breakpoints (on the x-axis of the nanoseconds) were computed by ordinary least squares method, and Bayesian Information Criterion (BIC) was applied in order to get the best number of optimal breakpoints (i.e., optimal segment partition); in SIB BW1280R the second one a F statistics were computed in order to get and test the strongest optimal breakpoint. Of note, the null hypothesis of the F statistics was the absence of a structural change in the Fab-RBD trajectory, so P-value less than 0.05 revealed the presence in significant way. For both sub-analyses, the corresponding affinity score values were also returned. This part of analysis was performed by R package strucchange. Next, a broken-line regression analysis was performed in order to further characterize the time series clustering, by discovering the nanosecond breakpoint beyond which the time series stabilized the affinity score trend, potentially. In detail, this method was applied to the Fab-RBD complex trajectories by using a line-threshold model: the expected trend of the model was linear until to an estimated threshold on the nanosecond, then it stabilized on a constant value of affinity after such threshold. Briefly, the underlying model consisted of two straight lines joined at a breakpoint. This part of analysis was performed by R package lm.br. Finally, two stratified analyses of the dose response trend modelling were performed by using joint dose-response models, the first one was stratified by variant, the second one by antibody. In the first analysis, for each variant we statistically compared the different antibodies by the model parameters (b,c,d,e,). In the second analysis, for each antibody we statistically compared the different variants by the same parameters. In this application, the model returned the parameters b, c, d and e, whose interpretations corresponded to the four pairwise comparisons (differences) among the parameters b, c, d and e of the stratified dose response models. This part of analysis was performed by R packages drc, devtools and multcomp. 3.2 – Analysis of residues The third step of the statistical analysis was the analysis on the residues. In this phase, the residues for each Fab-RBD complex were investigated in relation to their affinity score values. Residues data were provided in a data matrix containing three independent replicates of time series on 0-100 nanoseconds. In detail, in a preliminary step the mean for every nanosecond on the three independent replicates was computed in order to obtain a mean trajectory. Hence, the median of the affinity score values of the mean trajectory were computed. The median was used because of its robustness to the outliers and considering the decreasing trend (approximately) of the trajectory, whose starting point at baseline was affinity score=1. Finally, the medians of the affinity score were also shown by heatmaps (red: low affinity score; green: high affinity score) stratified by SIB BW1280R antibody. Of note, the rows were the RBD-Fab complex, and the columns were the residues. 4 - Antigen-Antibody complex stability prediction modules: lstmContacts The lstmContactslibrary (https: / / github.com / fernandoPalluzzi / lstmContacts) has two main modules: the contact similarity search module (R code) and the prediction module, based on long short-term memory (LSTM) recurrent neural networks (Python code). The R module takes the input antigen-antibody complex and generates an affinity contact library that can be used by the LSTM module for the antigen-antibody affinity time series prediction. The input complex is represented by a list C = [v1, v2, …, vn] of n contacts. Each j-th contact is a vector vj= [x0j, x1j, …, xmj] of mjelements, where x0j is the antibody residue interacting with the x1j, …, xmjresidues on the surface of the RBD of the Spike protein variant. Every contact corresponds univocally to a vector ai = affinity(vj) of 101 affinity score values, ranging from 0 to 1, and corresponding to the 101 nanoseconds of the molecular dynamics simulation. Since the dynamics are run in triplicates, we have three affinity vector replicates per contact: ai1, ai2, and ai3, where ai is computed as the element-wise mean of the three vectors. The preprocess(vj) function of the search module assigns the best possible ai vector to the vjcontact. If the dynamics for vjare in our contact library, the assignment is referred to as an exact match. Otherwise, the [x0j, x1j, …, xmj] residues of the vjcontact are searched by similarity. Each residue xkj(with k = 0, 1, …, m) is firstly searched in a nearby position of the polypeptide chain. If this search fails, the algorithm seeks for a residue with similar chemical properties (referred to as “group”), according to table 3 (lstmContact object contact.groups). Table 3. Table reporting the characteristic group of each aminoacid, based on the chemical properties of their side chains. Residues belonging to the same group are matched in a similarity-by-group search. SIB BW1280R For the antibody residue x0j of vj, the similarity search is further constrained within the original Fab chain (either heavy or light). A warning level based on powers of 2 is used to determine the search results quality: 23 = 8, exact match failed; 22 = 4, residue match (same Fab chain as the input) failed; 21 = 2, group search failed; 20 = 1, residue match (Fab chain mismatch) failed. A warning level below 8 means that an exact match was found, whereas a warning level between 8 and 14 indicates a non-exact match. If the warning level reaches 15, the search failed at every level and the algorithm cannot go further. This warning system enables fast user monitoring of the prediction quality. By default, the search module shows the search status, commenting on results quality in human readable language. Two criteria are used to further refine similarity search results. Firstly, each contact is classified on the base of the available interaction class. An interaction class specifies which part of the surface of the RBD is involved in the antigen-antibody complex formation. By default, the lstmContact object contact.groups specifies the four available classes (I-to-IV) for the Spike Fab-RBD interaction. Each residue is assigned (either exactly or by similarity) to one or more classes. Then the odds ratio OR = odds(i-th class) / odds(not i-th class) is calculated. The contact is assigned to the class with the highest OR, referred to as the majority class. A confirmation of the goodness of the search results is when all the contacts of a complex belong to the same class. SIB BW1280R Secondly, to find the antigen variants that best fit the input complex, raw results are further refined by checking for possible mutation hotspots. The lstmContacts object contact.mutations contains mutational information for the Spike variants alpha, beta, delta, and omicron, with respect to the wtprotein. The presence of one or more residues from either the wtor the mutant Spike variant may reveal similarities with the input contact(s). Contact search can be applied iteratively, for each contact of the antigen-antibody complex C, through the contacts(C) function. The output of this function is a list of attributes that collectively describes the antigen-antibody complex, providing information to build the input dataset for the LSTM module from the contact library (by default, the lstmContacts object contact.data). The output contains: (i) the list of known antibodies that might interact with the given antigen residues, (ii) the list of known variants that might interact with the given antibody residues, (iii) the known antibody residues that might have similar affinity to those of the input complex, (iv) the interaction class(es) of the contacts in the input complex, (v) the warning level raised by each of the contacts in the input complex, and (vi) the quality level of the contact(s) search (3: exact, 2: good similarity search, 1: suboptimal similarity search). 5 - LSTM module data preparation Before launching the LSTM module, an optional step involves the preparation of the training data. This module uses an encoder-decoder strategy, in which the time series is divided into intervals of equal size (by default, 5 nanoseconds). Each interval is used to predict the next one, up to the end of the time series. If the time series length is not a multiple of the interval size, the remaining time points are removed from the end of the series. The training set consists in entries composed by a source interval (i.e., the input sequence) and a target interval (i.e., the sequence that should be predicted from the source one). Here, both the target and the source sequences are vectors of affinity scores (one value for each time point). The training data can be prepared from the internal contact library by using the prepareLibrary function: data <- prepareLibrary(data = contact.data, chunk = 5) The chunk argument defines the size of the time interval in nanoseconds. The object data is a data.frame with the following attributes: source sequences (data$a), target sequences (data$y), involved antibody residue (data$res), antibody (data$antibody), variant (data$variant). These attributes can be used to filter subsets of the training data. If an affinity score threshold is given (argument a0), an optinonaldata$group attribute will SIB BW1280R be added (this will be 0 for a stable contact and 1 for an unstable one). The LSTM input file can be then generated from this data.frame as a tab-separated text file. The lstmContacts software already comes with two learning sets derived from the internal contact library, with interval size 5 and 10 nanoseconds (contactLibrary_t5.txt and contactLibrary_t10.txt, respectively). If one of these datasets are used, the library preparation step is not required. 5.1 - Running the LSTM module The LSTM module is written to run under UNIX-based environments. This module requires the installation and activation of the Conda(https: / / docs.conda.io / en / latest / ) environment tfgpu_env, based on TensorFlow(https: / / www.tensorflow.org) and available at the lstmContacts repository. The prediction can be now performed in four steps: # 1. Training set preparation tset = contactLibrary(filename = "~ / contactsCore / contactLibrary_t5.txt", variant = variants_list, antibody = antibody_list) # 2. Model definition and training model, encoder, decoder = lstmTraining(tset, n = 101, units = 128, epochs = 100) # 3. Prediction set pset = contactLibrary(filename = "~ / contactsCore / contactLibrary_t5.txt", variant = validation_variant, antibody = validation_antibody) # 4. Prediction running profile = lstmProfile(pset, encoder, decoder, t0 = 5, t1 = 5, n = 101, method = "median") The lstmProfile function returns the predicted affinity score time series and prediction accuracy value. The variants_list and antibody_list objects are given by the output of the SIB BW1280R contacts function, from the search module. At step 2, the model is defined through the training set generated at step 1, the number of temporal features (n), the number of units of the LSTM (units), the number of training epochs (epochs; i.e., the number of forward and backward propagation cycles). The prediction set can be defined either to predict a known antibody-variant interaction (e.g., for validation purposes) or an unknown one. In case of an unknown variant, the variant and antibody arguments can be omitted and the input file may not have the y column (i.e., the expected outcome field). Finally, the predicted affinity time series can be generated using the lstmProfile function, taking the prediction (validation) set, the encoder and decoder generated during the training (step 2), the number of source (t0) and target (t1) intervals, the total number of features (n), and the method used to combine the predicted contacts into a single affinity time series (median or mean are available). 5.2 - LSTM validation The LSTM validation is based on a leave-one-out cross-validation (LOO-CV) scheme. At each iteration of the LOO-CV, the training set T is built by leaving out a given complex C. The model is then defined through the lstmTraining function, using T and default arguments (n = 101, units = 128, epochs = 100). The validation set is then composed by the source contact intervals belonging to complex C. Each k-th source interval Ck0 is predicted and compared against the given target Ck1 by computing the element wise absolute difference between them: ^= Σ(^^1^−^^0^), with L being the interval length (in the current study, L = 5 nanoseconds). If any of the ^^1^−^^0^≤^ and ^≤^^^^ (by default, c = 3 and dmax= c*L), the prediction is correct; it is wrong otherwise. The prediction accuracy is defined as the ratio n.correct / n.predictions. 6 - Affinity threshold calculation The global affinity score threshold is calculated using the R package OptimalCutpoints (version 1.1-5), as follows: library(OptimalCutpoints) optimal.cutpoints(X = "affinity", status = "y", tag.healthy = 0, methods = "SpEqualSe", data = MD) SIB BW1280R The MD object is a data.frame reporting the affinity values (attribute affinity) of each available molecular dynamics simulation. The attribute y is a binary vector equal to 0 if a given affinity value comes from a stable molecular dynamics simulation, and 1 if the value comes from an unstable simulation. The stability values are derived from the affinity time series cluster analysis (clusters 1 and 3 are stable, while clusters 2 and 4 are unstable). The criterion used to define the optimal cutpoint is the affinity value at which the equality between sensitivity and specificity is reached. This package allows also to compute point area under the ROC curve (AUC) values and related 95% confidence intervals. Protein Interactions Analysis To explain monoclonal antibodies behavior against different Sars-Cov-2 variants, we investigate the interface interactions of the macromolecular complexes. Here, the procedure for each Fab-RBD complex: (i) To have free access to files generated with Molecular Dynamics Simulations, enter the following command: $ sshfs 14802292@10.160.65.20: / home / 14802292 ~ / 14802292 (ii) Upload on VMD “MD.pdb” file and “trjadj.xtc” for trajectory coordinates. (iii) Save coordinates of frames corresponding to nanosecond where the trend of curve dose-response (refers to Statistic Methods) undergoes a significant change, in a pdb file named “RBDexpPDBcodeNS.pdb”. 1000 frames for 100 nanoseconds, first frame to select: ns x 10 – 5, last frame: ns x 10 + 5. (iv) Import and process pdb file generated on Maestro (Schrödinger) with Protein Preparation Wizard. Refine the protein structure through H-bond assignment and restrained minimization. All parameters were left by default but force field OPLS_2005. (v) Perform the analysis with Protein Interaction Analysis Panel: Chain A (RBD) was assigned to Group 1, both chain B and chain C (heavy and light chain of Fab respectively) to Group 2. The interactions between two groups are calculated listing closest interaction neighbors within 4 Å, ignoring interactions between backbone atoms. The criteria for detecting bonds were established in advanced options panel, setting a maximum distance of 3,0 Å for Hydrogen Bonds (HB), 5 Å for Salt Bridges, 5.5 Å for pi-pi stacking interactions. The parameters not specified were left by default. SIB BW1280R (vi) Export table in a csv file - “RBDexpPDBcodens.csv”: residue | closest residues | distance | Interactions | # HB | # Salt Bridges | # Pi stacking | Disulfides | # vdW clash | Surface Complementary | Buried SASA Statistical Analysis protocol Statistical analysis explored the pattern of the molecular dynamics simulations in relation to the affinity score across nanosecond for each RBD-Fab complex. The RBD-Fab simulation trajectories to be analyzed were got by averaging by nanosecond the three independent replicates, as reported in the previous steps of the protocol. The analysis encompassed 3 steps: i) a cluster analysis for time series; ii) a dose response trend analysis for characterizing the achieved clusters; iii) an analysis of the residues in order to investigate the potential RBDs causing low performances clusters of RBD-Fab complexes. All the analysis were performed by R software. i) Cluster analysis for time series was performed to detect the complexes (RBD-Fab complexes trends) that could be clustered. This step was performed by R packages TSclust and ggdendro by applying the following key commands (in courier-new font). Library(Tsclust) #load the package Tsclust Include all the 0-100 time series (one by RBD-Fab complex) of affinity score in a data matrix (class data.frame), whose rows were the nanoseconds (0 to 100), whereas the columns were the RBD-Fab complexes: Time_series_dataframe<-data.frame() Compute the dissimilarity matrix by the Dynamic Time Warping (DTWARP) method: dist_ts<- Tsclust::diss(SERIES = t(Time_series_dataframe), METHOD = "DTWARP") Perform the hierarchical cluster analysis on the dissimilarity matrix by using a complete- linkage agglomeration method: hc<- stats::hclust(dist_ts, method="complete") Evaluate the dendrogram in order to define the cut generating the number of clusters plot(hc) SIB BW1280R Cuts the dendrogram resulting from hclustinto desired number(s) of clusters (k): hclus<- stats::cutree(hc, k=…) Load the package ggdendro: library(ggdendro) Extract cluster info from dendrogram: hcdata<- ggdendro::dendro_data(hc) ii) Secondly, a dose-response trend analysis was applied to characterize the RBD-Fab complex clusters. Concerning that, the nanoseconds (i.e., doses) were evaluated on the x-axis, whereas the affinity score (trend) was plotted on the y-axis. In this step, a change point analysis was also performed to detect the structural changes of affinity score across nanoseconds for each RBD-Fab complex. Finally, a broken-line analysis was performed in order to further characterize the time series clustering, by evaluating the nanosecond breakpoint beyond which the time series stabilized the affinity score trend, potentially. Briefly, the underlying model consisted of two straight lines joined at a breakpoint. In this step, the following R packages were applied: drc, devtools, and multcomp for the dose-response trend analysis, whereas strucchange and lm.br were applied for the change point analysis and broken-line analysis, respectively. The key commands (in courier-new font) for carrying out this step are reported below: Dose-response trend analysis: Import the file.csv of the data-matrix in R environment data<-read.csv2(file.choose(),header=T,sep=",",dec=".") Tab.3. Format of the data matrix for each RBD-Fab complex Load the packages drc, devtoolsandmultcomp library(drc) SIB BW1280R library(devtools) library(multcomp) We fitted a dose response model by a 4-parameter log-logistic function, where the parameters were the slope (b), the lower (c) and upper (d) limits of the affinity score expected from the model, and the ED5 (e), i.e., the nanosecond corresponding to the expected half-way affinity score between the upper and lower limit. Beta_Ly_Cov555.LL.4 <- drm(affinity_score ~ nanosecond, data = subset(data, Variant_Antibody == "Beta-Ly-Cov555"), fct = LL.4(names = c("Slope", "Lower Limit", "Upper Limit", "ED50")) Next, the change point analysis detected the structural changes of affinity score across nanoseconds for each RBD-Fab complex. In this phase, two sub-analyses were performed, by considering frameworks of breakpoint estimation and F (Chow) statistical testing: in the first one, multiple optimal breakpoints (on the x-axis of the nanoseconds) were computed by the ordinary least squares (OLS) method, and Bayesian Information Criterion (BIC) was applied in order to get the best number of optimal breakpoints (i.e., optimal segment partition); in the second one, a F statistics were computed in order to get the strongest optimal breakpoint. Get the multiple optimal breakpoints by the OLS method. The h parameter states the minimal segment size. This command returns also the corresponding affinity scores (on y-axis). breakpoints(Affinity_score ~ 1, h = 2) Then, compute the F statistics on the time series of the RBD-Fab trajectory fs.score<- Fstats(as.ts(Affinity_score) ~ 1) Extract the strongest optimal breakpoint, i.e., the nanosecond corresponding to the structural change point of the affinity score time series nanosecond_change_point<-breakpoints(fs.score)$breakpoints Get the value of affinity score corresponding to the strongest optimal breakpoint AS_nanosecond_change_point<- Affinity_score[breakpoints(fs.score)$breakpoints] Finally, in order to provide a further characterization, a broken line regression analysis was applied to the RBD-Fab complex trajectories by using a line-threshold model: the SIB BW1280R expected trend of the model was linear until to an estimated threshold on the nanosecond, then it stabilized on a constant value of affinity after such threshold. At first, we load the package lm.br library(lm.br) Hence, we fit the model m0.br<-lm.br(Affinity_score ~ nanosecond, type ="LT") "LT" which stand for line-threshold defined the trend just described. The model fitting returned the breakpoint on nanosecond scale m0.br$coef[1] and the corresponding affinity score value m0.br$coef[2] In addition, two stratified analysis of the dose response trend modelling were performed by using joint dose-response models, the first one was stratified by variant, the second one by antibody. Concerning the stratification by variant, such an example, we report the key commands of the analysis on the BETA variant. Firstly, we set the reference category of an antibody, e.g., AZD1061. data$Variant_Antibody<- relevel(factor(data$Variant_Antibody), ref = "BETA- AZD1061") Hence, we fit the dose response model by comparing the different antibodies with the reference, related to the variant Beta: BETA.LL.4 <- drm(Affinity_score ~ nanosecond, curveid = Variant_Antibody, data = subset(data, Variant == "Beta"), fct = LL.4(), pmodels = list(~Variant_Antibody, ~ Variant_Antibody, SIB BW1280R ~ Variant_Antibody, ~ Variant_Antibody)) Secondly, we statistically test the other pairwise comparisons of the model parameters (b,c,d,e,) among the different antibodies, related to the variant Beta: ") ") ") ") In the same way, we fit the stratified dose response model by antibody. Such an example, we report the key commands of the analysis on the antibody Ly-Cov555. Firstly, we set the reference category of a variant, e.g., Wild type. data$Variant_Antibody<- relevel(factor(data$Variant_Antibody), ref = "WT-Ly- Cov555") Hence, we fit the dose response model by comparing the different variants with the reference, related to the Ly-Cov555 antibody: Ly_Cov555.LL.4 <- drm(Affinity_score ~ nanosecond, curveid = Variant_Antibody, data = subset(data, Antibody == "Ly-Cov555"), fct = LL.4(), pmodels = list(~Variant_Antibody, ~ Variant_Antibody, ~ Variant_Antibody, ~ Variant_Antibody)) Then, we statistically test the other pairwise comparisons of the model parameters (b,c,d,e,) among the different variants, related to the antibody Ly-Cov555: ") ") SIB BW1280R compParm(Ly_Cov555.LL.4, "d", "-") ") iii) The third step of the statistical analysis was the analysis of the residues. In this step, the residues for each RBD-Fab complex were investigated in relation to their affinity score values. Residues data were provided in a data matrix (.csv), that was processed as data.frameobject in R programming, containing three independent replicates of time series on 0-100 nanoseconds. Here, we show the key commands of the residue R50 of the heavy chain of the RBD- Fab complex WT_LyCov555. Lycov555_R50hWT<-data.frame (Rep1=datiLycov555$R50hWT[1:101], Rep2=datiLycov555$R50hWT[102:202], Rep3=datiLycov555$R50hWT[203:303]) In the first step, the mean for every nanosecond on the three independent replicates was computed in order to obtain a mean trajectory. apply(Lycov555_R50h1WT,1,mean) Hence, the median of the affinity score values of the mean trajectory were computed. The median was used because of its robustness to the outliers and considering the decreasing trend (approximately) of the trajectory, whose starting point at t0 was affinity score=1. median(apply(Lycov555_R50d1WT,1,mean)) Finally, the medians of the affinity score were also shown by heatmaps (red: low affinity score; green: high affinity score) stratified by antibodies. Of note, the rows were the RBD- Fab complex, and the columns were the residues. Long short-term memory recurrent neural networks protocol 1. Requirements The lstmContacts (https: / / github.com / fernandoPalluzzi / lstmContacts) has been used to evaluate the predictive power of our molecular dynamics simulation. The software has two main modules: an R-based search module for contact similarity search and contact data preprocessing, and a Python module based on long-short term memory (LSTM) SIB BW1280R recurrent neural networks, used to predict the affinity time series of an antigen-antibody complex. The software is designed for Unix-based systems. The search module only requires the base R environment (>= 4.0), that can be installed following the inst https: / / www.r-project.org. The lstmContacts library already comes with the required contact data library and the supplementary annotations. The LSTM-based library requires the installation of Python3.8 and the tensorflow-gpu library vershttps: / / pypi.org / project / tensorflow-gpu), together with core computational Python libraries, including numpy, scipy, pandas, and matplotlib. To facilitate the correct execution of the software, the installation of Conda is rhttps: / / docs.conda.io / projects / conda / en / latest / user-guide / install). The lstmContacts repository provides a tfenv.yml file that can be used to install all the required dependencies: # Activate the base Conda environment condaactivate # Install the TensorFlow environment condaenv create -f tfenv.yml # Deactivate the base environment condadeactivate # Activate the TensorFlow-GPU environment conda activate tfgpu_env # Start a Python session to verify that python 3.8 is in use python # Type CTRL-D to exit and deactivate the environment condadeactivate For completeness, the list of Python imports is shown below: import re import pandas from random import randint, uniform SIB BW1280R from numpy import median from numpy import argmax from numpy import array from numpy import array_equal from keras.layers import Input from keras.layers import LSTM from keras.layers import Dense from keras.models import Model from keras.utils import to_categorical 2. Installation To install lstmContacts, it is sufficient to clone or copy its repository within a given directory (e.g., the home directory) and add its full path to the PYTHONPATH environment variable. This can be done permanently by modifying the .bashrc file, by adding the following line: export PYTHONPATH=$PYTHONPATH:~ / lstmContacts Please, be sure to have execution rights for the ~ / lstmContacts / lstmContacts.py file. 3. Running the R module for contact library search The search module can be used by starting an R environment and loading the following code: source("~ / lstmContacts / contacts.R") The input antigen-antibody complex can be specified as a list of contacts, as in the example below: x <- list(x1 = c("h.R50", "V483", "E484"), x2 = c("h.L55", "L452", "T470", "F490"), x3 = c("h.Y101", "E484", "F490"), x4 = c("h.R104", "Q493", "S494"), x5 = c("l.Y32", "F486", "Y489"), SIB BW1280R x6 = c("l.Y92", "F486", "Y489"), x7 = c("l.R96", "V483", "E484")) Each contact is a vector. The first element of the vector is the antibody residue, defined as a lower-case letter determining the FAB chain (“h” for heavy and “l” for light), followed by a dot, the single-letter amino acid code, and the position in the polypeptide chain. The other elements are the Spike receptor binding domain (RBD) residues interacting with the antibody one. The antigen-antibody complex search can be easily done with the contacts() function: R <- contacts(x) The variable R contains: the list of antibodies matching the input complex (R$antibody), the list of variants matching the input complex (R$variant), the antibody residues that will be extracted from the contact library to generate the predicted profile (R$ab.residues), the interaction class (R$class), the search results warning codes (R$warning.level), and the search quality level (R$contact.level). The interaction class specifies which part of the surface of the Spike RBD is involved in the complex formation. By default, the object contact.class specifies the four available classes (I-to-IV) for the Spike RBD-FAB interaction. Each residue is assigned (also by similarity) to one or more classes. Then the odds ratio p = odds(i-th class) / odds(not i-th class) is calculated. The contact is assigned tothe class with the highest p, referred to as the "majority class". This cleans up the extraction from the contact library, increasing the prediction specificity. The warning level ranges from 0 to 15 and determines the quality of the results: exact match found (0 to 7), similarity search enabled (8 to 14), search failed (= 15). Analogously, the contact level provides information about the search results type: exact match (= 3), high-quality similarity search (= 2), suboptimal similarity search (= 1). On the base of these indices, the software will raise human-readable warnings. The selected contacts can be now manually extracted from the contact library with the extractProfiles() function: profile <- extractProfiles(data = contact.data, antibody = R$antibody, variants = R$variant, residues = R$ab.residues) SIB BW1280R head(profile) where contact.data is the internal contact library. The profile variable includes all the profiles that best match the input contact. Although this procedure can be used to further inspect search results, it does not represent an actual prediction, which is instead done using the LSTM module. The quality of the LSTM prediction can be improved using the antibody and variant attributes provided by the search module (see section 5). 4. LSTM input data preparation Before launching the LSTM module, an optional step involves the preparation of the training data. This module uses an encoder-decoder strategy, in which the time series is divided into intervals of equal size (by default, 5 nanoseconds). Each interval is used to predict the next one, up to the end of the time series. If the time series length is not a multiple of the interval size, the remaining time points are removed from the end of the series. The training set consists in entries composed by a source interval (i.e., the input sequence) and a target interval (i.e., the sequence that should be predicted from the source one). Here, both the target and the source sequences are vectors of affinity scores (one value for each time point). The training data can be prepared from the internal contact library by using the prepareLibrary() function: data <- prepareLibrary(data = contact.data, chunk = 5) The chunk argument defines the size of the time interval in nanoseconds. The object data is a data.frame with the following attributes: source sequences (data$a), target sequences (data$y), antibody residue (data$res), antibody (data$antibody), variant (data$variant). These attributes can be used to filter subsets of the training data. The lstmContacts library already comes with two learning sets derived from the internal contact library, with interval size 5 and 10 nanoseconds. If one of these datasets are used, the prepareLibrary() step is not required. 5. LSTM running Firstly, we need to activate the TensorFlow-GPU environment and start a Python session: conda activate tfgpu_env python SIB BW1280R Once inside the python console, we should load the LSTM module: from lstmContacts import * The analysis can be now performed in four simple steps. The training set can be prepared with the contactLibrary() function: tset = contactLibrary(filename = "~ / contactsCore / contactLibrary_t5.txt", variant = variants_list, antibody = antibody_list) Arguments variant and antibody are optional and can be used to increase the specificity of the training set according to the search module indications (see section 3). If these arguments are missing, the entire training set will be used. Model definition and training can be done simultaneously with the lstmTraining() function: model, encoder, decoder = lstmTraining(tset, n = 101, units = 128, epochs = 100) Besides the training set, this function takes the total number of features, the number of units of the LSTM (by default, 128 units are used), and the number of epochs (i.e., the number of forward-backward propagation cycles that are used to learn model parameters). By default, the epochs are set to 100, although the user may tune this argument, depending on the training set size and the available computational resources. The prediction set can be then defined similarly to the training set: pset = contactLibrary(filename = "~ / contactsCore / contactLibrary_t5.txt", variant = "beta", antibody = "7l7d") Finally, the affinity profile can be drawn using the lstmProfile() function: profile = lstmProfile(pset, encoder, decoder, t0 = 5, t1 = 5, n = 101, method = "median") This function requires: the set of source sequences to be predicted, the encoder and decoder generated during the training step, the size of the source (t0) and target (t1) intervals, the total number of features (n), and a method to combine single contacts SIB BW1280R predictions to obtain a single predicted antigen-antibody affinity time series (by default, the median is used). 6. Affinity score threshold calculation The global affinity score threshold is calculated using the R package OptimalCutpoints (version 1.1-5), as follows: library(OptimalCutpoints) optimal.cutpoints(X = "affinity", status = "y", tag.healthy = 0, methods = "SpEqualSe", data = MD) The MD object is a data.frame reporting the affinity values (attribute affinity) of each available molecular dynamics simulation. The attribute y is a binary vector equal to 0 if a given affinity value comes from a stable molecular dynamics simulation, and 1 if the value comes from an unstable simulation. The stability values are derived from the affinity time series cluster analysis (clusters 1 and 3 are stable, while clusters 2 and 4 are unstable). The criterion used to define the optimal cutpoint is the affinity value at which the equality between sensitivity and specificity is reached. This package allows also to compute point area under the ROC curve (AUC) values and related 95% confidence intervals.

Claims

SIB BW1280R CLAIMS 1. A computer implemented method to predict which available monoclonal antibody against a target protein of a virus binds with the highest affinity a new variant of the target protein of said virus, the method comprising the following steps: i) preparing in silico a first dataset comprising all the available structures of complexes of said antibodies with said target protein and its known variants (virus-antibody complexes); ii) performing more than one molecular dynamic simulations of the complexes of said first dataset; iii) performing a statistical analysis on the molecular dynamic simulations of step (ii) obtaining an affinity contact library; iv) refining the affinity contact library of step iii) in order to provide a training set for the prediction of the interaction affinity of said virus-antibody complexes; v) training a neural network with said training set obtained in step iv) to enhance the accuracy of the prediction of the affinity of said virus-antibody complex thereby obtaining a neural network trained on the affinity contact library; vi) running said trained neural network using as input the structure of the complex between the available antibodies and said new variant of the target protein of the virus, thereby identifying the available monoclonal antibody specific for the new variant of said target protein.

2. The method of claim 1, wherein said virus is a Coronavirus.

3. The method of claim 1 or 2, wherein said virus is a Sars-Cov2.

4. The method according to any one of the preceding claims, wherein said target protein is the Receptor Binding Domain (RBD) of the Sars-Cov2.

5. The method according to any one of the preceding claims, wherein said virus-antibody complexes are Fab-RBD complexes (Fragment antigen-binding – Receptor Binding Domain).

6. The method according to any one of the preceding claims, wherein said first dataset comprises the complex of each available monoclonal antibody with each available variant of the target protein.SIB BW1280R 7. The method according to any one of the preceding claims, wherein said antibodies have a marketing authorization for the treatment of a disease caused by said virus.

8. The method according to any one of the preceding claims, wherein said Fab-RBD complexes are complexes between the Fabs of all available monoclonal antibodies against the RBD of all the variants of the Sars-Cov2.

9. The method according to any one of the preceding claims, wherein said step (i) of preparation of complexes of said antibodies with said target protein and its known variants comprises the steps of: (i-a) downloading crystal structures of complexes of said Fab-RBD complexes from a database, such as PDB (Protein Data Bank); (i-b) reassigning chain names by identifying and renaming the chains in the crystal structure to correspond to the Fab and RBD components of the complex; (i-c) removing HETATM residues by identifying and removing any non-standard amino acids or other molecules present in the crystal structure that do not correspond to the Fab or RBD; (i-d) implementing RBD mutations by introducing targeted mutations to the RBD component of the complex to simulate various binding scenarios; (i-e) adding and optimizing hydrogen bonds between residues by identifying and adding missing hydrogen bonds between the Fab and RBD residues to optimize the stability of the complex.

10. The method according to any one of the preceding claims, wherein said step (ii) of performing more than one molecular dynamics simulations of said Fab-RBD complexes of said first dataset comprises the steps of: (ii-a) setting up simulation parameters; (ii-b) minimizing and equilibrating the systems by adjusting the position of the atoms thereby reducing the initial energy of the system and bringing it to a stable state before running the simulation; (ii-c) running simulations for 100 ns under NPT conditions by running the simulation for a specific duration and under pressure of 1 atm, and at a temperature from 20°C to 40°C, preferably 30°C, to obtain a representative sample of the system's dynamics.SIB BW1280R 11. The method according to any one of the preceding claims, wherein said step (iii) of performing statistical analysis on the molecular dynamic simulations of step (ii) comprises the steps of: (iii-a) cluster analysis by grouping together similar conformations of the complex according to their stability and / or binding affinity in order to detect the presence of common classes of complex, based on their time series stability and / or binding affinity; (iii-b) linear regression modelling and dose response trend analysis by identifying the relationship between changes in the RBD and corresponding changes in binding stability and / or affinity in order to characterize and better interpret the common classes of complex detected in previous analysis step by measuring linear regression modelling and dose response trend analysis parameters; (iii-c) analysis of residues causing low performance clusters by identifying which residues in the RBD are responsible for lower binding stability and / or affinity in order to identify the residues causing low performance classes of complex (clusters).

12. The method according to claim 11, wherein said linear regression modelling parameters is at least one selected from the affinity score value at baseline, and the affinity score expected variation per an additional nanosecond, and said dose response trend analysis parameter is at least one selected from the slope (of the dose response modelling), the lower and upper limits of the affinity score expected from the model and the nanosecond corresponding to the expected half-way affinity score between the upper and lower limit (ED50).

13. The method according to any one of the preceding claims, wherein said step iv) of refining said affinity contact library based on the statistical analysis results comprises the steps of: (iv-a) refining the affinity contact library by interaction class by grouping together similar interactions between the Fab and RBD residues according to their contribution to the overall binding stability and / or affinity of the complex, considering at least one of side chain biochemical similarity, expected spatial localization, and matching with mutations in known variants by means of a similarity search module thereby obtaining a first refined affinity contact library; (iv-b) refining said first refined affinity contact library obtained in (v-a) by mutation hotspots by identifying and prioritizing specific amino acid residues in the RBD that are most likely to affect the binding stability and / or affinity when mutated, by evaluating theSIB BW1280R contact priority on the base of the odds ratio of belonging to an interaction localization class by means of a filtering module which provides a quality value of the similarity search on a scale from 0 (totally unreliable) to 15 (exact match in the contact library), in order to further refine the affinity contact library; (iv-c) stratifying the further refined affinity contact library into a training set and a validation set by dividing the library into two distinct sets in order to train and test a neural network for predicting the stability of the Fab-RBD complex by means of a leave-one- out-cross-validation strategy where for each training-prediction iteration, a set of contacts corresponding to a single Fab-RBD complex is excluded from the training set and used as validation set.

14. The method according to any one of the preceding claims, wherein said neural network is a long short-term memory (LSTM) recurrent neural network.

15. A computer program comprising instructions which, when the program is executed by a computer, cause the computer to carry out the steps of the method according to claims 1-14.

16. A computer-readable storage medium comprising instructions which, when executed by a computer, cause the computer to carry out the steps of the methods according to claims 1-14.

Images