BACKGROUND
structure-based site-directed mutagenesis is a widely used approach to elucidate and modify specific aspects of protein function and to investigate the properties of protein-ligand interactions.
prediction of mutants with desirable properties is often obtained by rational design of a few specific position in the region of interest.
this approach requires an in-depth knowledge of physicochemical, structural and functional properties of the protein, which may not be exhaustively provided by x-ray crystallography data.
a different strategy is to apply the so-called scanning mutagenesis based on the systematic replacement of all residues involved in a specific function.
several examples of alanine scanning mutagenesis experiments  <cit>  are reported in literature that illustrate its relevance for the study of molecular determinants of ligand binding  <cit>  and protein function  <cit> .
even more informative is the use of the saturation mutagenesis, where all  <dig> natural amino acids are tested in place of the wild type residue, in or near the functional site.
'in vitro' scanning saturation mutagenesis has been applied to several case studies and has proved to be useful to modulate protein properties such as substrate enzyme specificity  <cit>  or to identify key residues for catalytic mechanisms  <cit> .
a particular aspect that is amenable to investigation with random mutagenesis is the study of protein-ligand interaction.
several examples are reported in literature that show how this approach has contributed to the clarification of ligand interaction in g protein coupled receptors  <cit> , in nuclear receptors  <cit>  and in epha3/ephrin-a <dig>  <cit> , to cite a few examples.
however, the large number of mutants that are obtained from the systematic replacement of all the residues of a functional site makes the experimental procedure laborious and time consuming.
this is especially true when dealing with large functional sites, because it requires the implementation of specific high-throughput methods for mutant screening or selection  <cit> .
'in silico' simulation of mutagenesis and screening may represent a possible way around this problem.
by exploiting the capabilities of high performance computers, it is possible to design and test a large numbers of variants with altered properties in order to restrict the number of models that are worth testing experimentally.
computational protein design methods have been employed in protein engineering to increase enzymes efficiency  <cit>  or antibiotic resistance: for example hayes et al.
<cit>  applied a protein optimization strategy to increase the resistance of bacteria toward the antibiotic cefotaxime by optimizing tem- <dig> β-lactamase.
moreover, computational design tools have been developed and have proved to be successful in the design of proteins with enhanced stability and specificity, and in engineering novel protein functions  <cit> .
we present here a high performance semi-automated computational procedure suited to the study of key residues of protein-ligand interactions, which is able to dissect the contribution of different residues to ligand binding.
the procedure performs a saturation mutagenesis of all residues involved in ligand binding, and the subsequent evaluation of the effect of amino acid substitutions on ligand affinity by docking simulation and on protein structure stability.
an example of application to endothelial protein c receptor is reported.
methods
pipeline implementation
a semi automated pipeline was designed to manage mutant modelling, docking simulations and data analysis procedures by integrating public software through a number of in-house developed perl scripts.
mutant modelling
'in silico' side chain replacement and modelling were carried out with the routine 'mutate_model' of modeller 9v <dig>  <cit> .
'mutate_model' introduces a single point mutation in a user-specified residue and optimizes the mutant side chain conformation by conjugated gradient and by molecular dynamics simulation.
docking
autodock <dig>   <cit>  was used for docking simulation handling both ligand and mutated side chain as flexible.
autodocktools  () facility supported the protein and ligand set up for docking.
the automatic procedure, based on a perl script, employs the following adt routines:
- 'prepare_receptor4' to add polar hydrogens, assign kollman charges and convert the protein pdb file in pdbqt ;
- 'prepare_flexdocking' and 'prepare_flexreceptor' are applied to protein and ligand to obtain the input files for flexible docking simulation;
- 'prepare_gpf' and 'prepare_dpf' generate the parameters file for autogrid  and for autodock , respectively.
docking simulations were performed with lamarkian genetic algorithm which is reported as the most efficient and reliable method of autodock  <cit> .
stability evaluation
the 'stability' routine of foldx algorithm  <cit>  was applied to evaluate the total energy of wild type structure and mutated models by the estimation of the interactions that contribute to the protein stability.
computational resource
docking simulations were performed on a shared linux cluster of  <dig> opteron amd cores  <dig> at  <dig>  ghz, provided by the eurotech group , dedicated to bioinformatics applications.
the system is composed of  <dig> chassis of  <dig> diskless blades, each equipped with an infiniband 4× network card and  <dig> gb of ram.
the cluster nodes are accessed through the openpbs queuing system.
case study
protein structure and preparation
the crystal structure of the endothelial protein c receptor , in a complex with phosphatidylethanolamine , and the gla domain of protein c solved by x-ray crystallography at  <dig>  Å  <cit>  was retrieved from the protein data bank  <cit> .
chain a with the corresponding ligand  molecule was used for this study.
the gla domain, crystallization water, ca ions, n-acetyl-d-glucosamine  and 2--2-deoxi-a-d-glucopyranose  were removed from pdb file.
polar hydrogens were added at  <dig>  ph, using insightii .
docking simulation
the adt graphics interface was employed for manual preparation of the ligand to assign the flexibility, add polar hydrogens and load gasteiger charges.
<dig> rotatable bonds were assigned to the pty ligand.
a grid box of dimensions  <dig> Å ×  <dig> Å ×  <dig> Å was constructed around the binding site, based on the co-crystallized ligand.
ligand flexible docking simulations were performed with  <dig> runs and  <dig> energy evaluations per run, while the other parameters were set to the default.
RESULTS
pipeline description
the complete scheme of the computational procedure used in this work is depicted in figure  <dig>  a preliminary step is performed using ligplot  <cit>  to identify the residues involved in ligand interactions.
the pipeline has a modular structure to improve its adaptability to different datasets.
the backbone of the pipeline is a mysql database where the results of each computational step are stored .
the pipeline is structured in two main stages devoted to the modelling and analysis of single and double mutants respectively .
the first stage begins by single mutant modelling starting from the receptor atomic coordinates and the list of residues identified in the preliminary step.
this first step makes use of the modeller routine 'mutate_model' that performs the side chain substitution and refinement.
a perl script is used to iterate the substitution of each residue with all the amino acids to achieve the saturation mutagenesis.
for each mutant, a pdb file is obtained and stored in the database.
in the following step, the pipeline prepares the input files for autogrid and autodock <dig>  according to user specified parameters: pdb files are converted to pdbqt format, adding polar hydrogens and assigning kollman charges.
substituted residues are set for flexible docking.
the ligand is manually prepared as illustrated in the methods section and given as input to the pipeline.
in the next step, autodock simulations are submitted to be run in parallel on the linux cluster.
docking results are parsed to extract, for each simulation, the most representative conformation.
this corresponds to the best energy conformation included in an autodock cluster of at least  <dig> elements and showing a binding energy within one kcal/mol compared to the first ranked solution.
this criterion is arbitrarily set in order to obtain a solution from a well represented group of conformations with the best scored results.
however, when these conditions are not satisfied , the first ranked solution is accepted.
for each selected conformation, the corresponding values of binding energy, cluster population, inhibition constant , root mean square deviation  and atomic coordinates are stored in the database.
in order to select mutants with a modified affinity for the ligand compared to the wild type, the output data are filtered taking the wild type ki value as threshold.
structural stability of selected mutants is evaluated using the foldx program that computes, for each conformation, an estimated free energy values.
the free energy difference between the mutant and the wild type  is automatically calculated and stored in the database.
pair combinations of single substitutions are analyzed in the second stage of the pipeline  to investigate the occurrence of synergistic effects in double mutants.
residues chosen for double mutants modelling are obtained from the results of ki filtering in first stage screening.
a perl script manages the residues combination, removing the redundant pairs.
the following pipeline steps, simulations, data analysis, filtering and storage, are performed as in the first stage.
the database is designed to collect and organize all the output data which remain available for further study.
specific queries can be applied in order to compare binding capabilities and stability of all the mutants, as illustrated in the example below.
case study
we used the described strategy to study the ligand binding of the endothelial protein c receptor .
epcr is a transmembrane glycoprotein homologous to the major histocompatibility complex  class 1/cd <dig> family  <cit> .
epcr plays a well-characterized role in the coagulation cascade being the receptor of the anticoagulant factor protein c  and participating to its activation mechanism.
being involved in protein c pathway, epcr also takes part in the regulation of inflammation and apoptosis in endothelial cells  <cit> .
similarly to cd <dig> molecules, it binds a phospholipid ligand in a deep hydrophobic groove composed of an eight-stranded β-sheet floor and three almost parallel α-helices that lie above the floor and form a crevice where the lipid is located   <cit> .
the specific role of the lipid for the protein structure and function has not yet been clarified.
amino acid variants located in the lipid binding region of epcr may provide useful information about the role of this interaction in protein function.
however, mutations that impair or suppress the lipid binding have not been reported so far.
in this work we predicted mutants that have reduced binding affinity for the ligand, while maintaining structural stability.
the structure of human epcr, co-crystallized with pty, exhibits a wide pattern of interactions that extend the length of binding groove.
pty makes contacts with  <dig> residues: four are located on α-helix h <dig>  two on h <dig> and nine on h <dig>  one on the loop between h <dig> and β-strand  <dig>  and eleven residues on the strands of β-sheet floor .
the binding is mostly dominated by hydrophobic interactions with the hydrocarbon chain of the ligand; a single h-bond is established between the polar group of phospholipid and arg  <dig>
employing the pipeline described above, we have carried out an 'in silico' saturation mutagenesis of the epcr residues involved in lipid binding.
each modelled mutant was subjected to molecular docking simulation and stability evaluation.
the protein-ligand complexes obtained by docking simulations were screened by their inhibition constant  values in order to estimate the impairment of ligand binding ability induced by the mutation.
the overall procedure and a summary of the results are depicted in figure  <dig>  the  <dig> interacting residues, identified using ligplot, were saturated with the  <dig> remaining amino acids, producing  <dig> mutant models.
then, wild type and mutants were analysed by docking simulations.
the wild type re-docking correctly reproduced the x-ray ligand conformation with a rmsd of  <dig>  Å, a binding energy of - <dig>  kcal/mol and a ki value of  <dig> nm.
positions and substitutions critical for binding affinity were identified with a ki-based filter.
referring to the wild type ki value, we used a threshold of  <dig> μm  to discriminate between mutants with probable impaired binding ability.
the pipeline selected  <dig> mutants with ki value higher than the threshold, substituted in  <dig> out of the  <dig> starting positions.
pair combinations of the  <dig> substitutions were generated and  <dig> double mutants were modelled.
they were analyzed by the same docking protocol of the single mutants and filtered with a ki threshold of  <dig> μm.
this value was set an order of magnitude higher than for single mutants, with the purpose of selecting those combinations with a binding ability worse than the single ones.
with this filtering criterion,  <dig> double mutants were rescued.
the large reduction in the number of double mutants , suggested that most residue combinations did not cooperate to increase the ki value compared to the single substitutions.
the structural stability of single and double variants that passed the corresponding ki thresholds was evaluated using the foldx tool, directly integrated in the pipeline.
mutants ensuing a ΔΔg value higher than  <dig> kcal/mol were considered structurally unstable and removed from the selection  <cit> .
after the screening, the pipeline retrieved  <dig> single and  <dig> double mutants with a ΔΔg <  <dig> kcal/mol.
binding energy and ki values intervals are reported in figure  <dig>  the stability filter was passed only by mutants with a ki value within  <dig> and  <dig> μm for single and double, respectively.
this value was an order of magnitude higher than the respective ki thresholds.
mutants with high ki values, although potential candidates for binding affinity studies, were rejected by the ΔΔg filter and thus considered unstable.
at the end of this analysis,  <dig> key positions for the ligand binding in epcr were identified : position  <dig> and  <dig> are localized in the h <dig> α-helix, positions  <dig>   <dig> and  <dig> are localized in h <dig> α-helix and one, the position  <dig> on the eight-stranded β-sheet floor.
the most represented substitutions of these key positions involved hydrophobic residues, such as w, in positions  <dig>   <dig> and  <dig>  and l, in positions  <dig> and  <dig>
pair combinations in double mutants were evaluated in terms of synergistic effect where the ki value of double mutant was larger than the sum of the corresponding single ones.
to estimate the synergistic effect, we applied the following equation:
where ki <dig> and ki <dig> represent the ki value of the first and the second single mutant, and ki1+ <dig> the ki value of double mutant.
kid is the difference between the effect of the double mutation  and the sum of the two corresponding single mutations .
results are shown in figure  <dig>  we found that  <dig> of the  <dig> double mutants displayed a partial additive effect as the ki value was included between the more damaging mutation and the sum of the two corresponding single mutants    <cit> .
the remaining  <dig> double mutants had ki values higher than the sum of the two corresponding single mutants , to which we assigned a synergistic effect compared to the single mutation .
we observed that the highest ki values were associated with double mutants resulting from replacement of a <dig> with either i or l and f <dig> with a small polar/apolar side chain.
between them, the highest synergistic effect was shown by a31i/f164s mutant , while the largest effect in terms of binding energy  and ki  value was associated to a31i/f164n mutant.
in conclusion, using the described approach, we have selected a limited number of single and double mutants of epcr for which we predict a reduced binding ability for the pty ligand and a structural stability comparable to the wild type.
the set of  <dig> selected single mutants and the  <dig> double mutants presenting synergistic properties can be proposed as rationally constructed candidates for experimental testing in cellular models.
pipeline performance
due to the computational load needed to accomplish all the simulations, a parallel infrastructure was employed .
in this test case, the simulation of the single and double mutant took on average  <dig> and  <dig> hours, respectively.
the whole challenge consisted of  <dig> docking experiments for the single mutants plus the wild type simulation and  <dig> for the double mutants.
depending on the resource availability, it was possible to employ  <dig> cpus on average.
therefore, in a period of about ten days we collected all the results covering about  <dig> years of computer activity.
CONCLUSIONS
we present a high performance semi-automated computational pipeline that integrates modelling and docking simulation algorithms for large scale mutant design.
additional analysis tools such as structure stability evaluation are also included.
the pipeline is organized into two main blocks for dealing with single and double mutants that are interconnected so that selected results of first analysis are used for the second part of the study.
the procedure is presently devoted to protein-ligand interaction studies, but can be adapted to investigations of protein-protein interaction.
the mutagenesis method is based on a saturation approach and does not require previous knowledge of functional or structural role of involved residues: it can therefore be applied to explore new binding features.
the implementation of the pipeline on a parallel computer infrastructure permits a greatly reduced computational time and makes it possible to test a very large number of mutants.
for example, in the reported text case we have modelled and tested  <dig> variants of epcr with the purpose of identifying single or coupled amino acid mutations that significantly impair the binding of the lipid ligand.
the whole procedure reduced the number of selected models to  <dig>  which corresponds to about 2% of the initial dataset.
selected variants provide information about the probable functional and structural effect of the mutated residues and can be proposed as subjects for experimental tests.
in addition, the large amount of data stored in the database after each step of computational analysis can be further retrieved for more detailed examination and study.
this pipeline can be useful for protein-ligand binding studies to identify key residues for the ligand interaction, to evaluate the effect of the substitutions on different ligands and to design selective mutants.
competing interests
the authors declare that they have no competing interests.
authors' contributions
fc developed the core of the pipeline, tested all the software for modelling and docking simulations, prepared the biological dataset, analyzed the output results and drafted the manuscript.
pd conceived and designed the study and critically revised the results.
im contributed to developing the pipeline, designing the database and enabling the parallel computation of the docking experiments on the cluster infrastructure.
lm supervised the project and made available the computational resources.
er coordinated the study and critically revised the results and manuscript.
all authors read and approved the final manuscript.

 structure-based site-directed mutagenesis is a widely used approach to elucidate and modify specific aspects of protein function and to investigate the properties of protein-ligand interactions. prediction of mutants with desirable properties is often obtained by rational design of a few specific position in the region of interest. this approach requires an in-depth knowledge of physicochemical, structural and functional properties of the protein, which may not be exhaustively provided by x-ray crystallography data. a different strategy is to apply the so-called scanning mutagenesis based on the systematic replacement of all residues involved in a specific function. several examples of alanine scanning mutagenesis experiments are reported in literature that illustrate its relevance for the study of molecular determinants of ligand binding and protein function. even more informative is the use of the saturation mutagenesis, where all natural amino acids are tested in place of the wild type residue, in or near the functional site.'in vitro' scanning saturation mutagenesis has been applied to several case studies and has proved to be useful to modulate protein properties such as substrate enzyme specificity or to identify key residues for catalytic mechanisms. a particular aspect that is amenable to investigation with random mutagenesis is the study of protein-ligand interaction. several examples are reported in literature that show how this approach has contributed to the clarification of ligand interaction in g protein coupled receptors, in nuclear receptors and in epha3/ephrin-a, to cite a few examples. however, the large number of mutants that are obtained from the systematic replacement of all the residues of a functional site makes the experimental procedure laborious and time consuming. this is especially true when dealing with large functional sites, because it requires the implementation of specific high-throughput methods for mutant screening or selection.'in silico' simulation of mutagenesis and screening may represent a possible way around this problem. by exploiting the capabilities of high performance computers, it is possible to design and test a large numbers of variants with altered properties in order to restrict the number of models that are worth testing experimentally. computational protein design methods have been employed in protein engineering to increase enzymes efficiency or antibiotic resistance: for example hayes et al. applied a protein optimization strategy to increase the resistance of bacteria toward the antibiotic cefotaxime by optimizing tem- β-lactamase. moreover, computational design tools have been developed and have proved to be successful in the design of proteins with enhanced stability and specificity, and in engineering novel protein functions. we present here a high performance semi-automated computational procedure suited to the study of key residues of protein-ligand interactions, which is able to dissect the contribution of different residues to ligand binding. the procedure performs a saturation mutagenesis of all residues involved in ligand binding, and the subsequent evaluation of the effect of amino acid substitutions on ligand affinity by docking simulation and on protein structure stability. an example of application to endothelial protein c receptor is reported. methods pipeline implementation a semi automated pipeline was designed to manage mutant modelling, docking simulations and data analysis procedures by integrating public software through a number of in-house developed perl scripts. mutant modelling'in silico' side chain replacement and modelling were carried out with the routine'mutate_model' of modeller 9v.'mutate_model' introduces a single point mutation in a user-specified residue and optimizes the mutant side chain conformation by conjugated gradient and by molecular dynamics simulation. docking autodock was used for docking simulation handling both ligand and mutated side chain as flexible. autodocktools () facility supported the protein and ligand set up for docking. the automatic procedure, based on a perl script, employs the following adt routines: -'prepare_receptor4' to add polar hydrogens, assign kollman charges and convert the protein pdb file in pdbqt; -'prepare_flexdocking' and'prepare_flexreceptor' are applied to protein and ligand to obtain the input files for flexible docking simulation; -'prepare_gpf' and'prepare_dpf' generate the parameters file for autogrid and for autodock, respectively. docking simulations were performed with lamarkian genetic algorithm which is reported as the most efficient and reliable method of autodock. stability evaluation the'stability' routine of foldx algorithm was applied to evaluate the total energy of wild type structure and mutated models by the estimation of the interactions that contribute to the protein stability. computational resource docking simulations were performed on a shared linux cluster of opteron amd cores at ghz, provided by the eurotech group, dedicated to bioinformatics applications. the system is composed of chassis of diskless blades, each equipped with an infiniband 4× network card and gb of ram. the cluster nodes are accessed through the openpbs queuing system. case study protein structure and preparation the crystal structure of the endothelial protein c receptor, in a complex with phosphatidylethanolamine, and the gla domain of protein c solved by x-ray crystallography at Å was retrieved from the protein data bank. chain a with the corresponding ligand molecule was used for this study. the gla domain, crystallization water, ca ions, n-acetyl-d-glucosamine and 2--2-deoxi-a-d-glucopyranose were removed from pdb file. polar hydrogens were added at ph, using insightii. docking simulation the adt graphics interface was employed for manual preparation of the ligand to assign the flexibility, add polar hydrogens and load gasteiger charges. rotatable bonds were assigned to the pty ligand. a grid box of dimensions Å × Å × Å was constructed around the binding site, based on the co-crystallized ligand. ligand flexible docking simulations were performed with runs and energy evaluations per run, while the other parameters were set to the default.  pipeline description the complete scheme of the computational procedure used in this work is depicted in figure a preliminary step is performed using ligplot to identify the residues involved in ligand interactions. the pipeline has a modular structure to improve its adaptability to different datasets. the backbone of the pipeline is a mysql database where the results of each computational step are stored. the pipeline is structured in two main stages devoted to the modelling and analysis of single and double mutants respectively. the first stage begins by single mutant modelling starting from the receptor atomic coordinates and the list of residues identified in the preliminary step. this first step makes use of the modeller routine'mutate_model' that performs the side chain substitution and refinement. a perl script is used to iterate the substitution of each residue with all the amino acids to achieve the saturation mutagenesis. for each mutant, a pdb file is obtained and stored in the database. in the following step, the pipeline prepares the input files for autogrid and autodock according to user specified parameters: pdb files are converted to pdbqt format, adding polar hydrogens and assigning kollman charges. substituted residues are set for flexible docking. the ligand is manually prepared as illustrated in the methods section and given as input to the pipeline. in the next step, autodock simulations are submitted to be run in parallel on the linux cluster. docking results are parsed to extract, for each simulation, the most representative conformation. this corresponds to the best energy conformation included in an autodock cluster of at least elements and showing a binding energy within one kcal/mol compared to the first ranked solution. this criterion is arbitrarily set in order to obtain a solution from a well represented group of conformations with the best scored results. however, when these conditions are not satisfied, the first ranked solution is accepted. for each selected conformation, the corresponding values of binding energy, cluster population, inhibition constant, root mean square deviation and atomic coordinates are stored in the database. in order to select mutants with a modified affinity for the ligand compared to the wild type, the output data are filtered taking the wild type ki value as threshold. structural stability of selected mutants is evaluated using the foldx program that computes, for each conformation, an estimated free energy values. the free energy difference between the mutant and the wild type is automatically calculated and stored in the database. pair combinations of single substitutions are analyzed in the second stage of the pipeline to investigate the occurrence of synergistic effects in double mutants. residues chosen for double mutants modelling are obtained from the results of ki filtering in first stage screening. a perl script manages the residues combination, removing the redundant pairs. the following pipeline steps, simulations, data analysis, filtering and storage, are performed as in the first stage. the database is designed to collect and organize all the output data which remain available for further study. specific queries can be applied in order to compare binding capabilities and stability of all the mutants, as illustrated in the example below. case study we used the described strategy to study the ligand binding of the endothelial protein c receptor. epcr is a transmembrane glycoprotein homologous to the major histocompatibility complex class 1/cd family. epcr plays a well-characterized role in the coagulation cascade being the receptor of the anticoagulant factor protein c and participating to its activation mechanism. being involved in protein c pathway, epcr also takes part in the regulation of inflammation and apoptosis in endothelial cells. similarly to cd molecules, it binds a phospholipid ligand in a deep hydrophobic groove composed of an eight-stranded β-sheet floor and three almost parallel α-helices that lie above the floor and form a crevice where the lipid is located. the specific role of the lipid for the protein structure and function has not yet been clarified. amino acid variants located in the lipid binding region of epcr may provide useful information about the role of this interaction in protein function. however, mutations that impair or suppress the lipid binding have not been reported so far. in this work we predicted mutants that have reduced binding affinity for the ligand, while maintaining structural stability. the structure of human epcr, co-crystallized with pty, exhibits a wide pattern of interactions that extend the length of binding groove. pty makes contacts with residues: four are located on α-helix h two on h and nine on h one on the loop between h and β-strand and eleven residues on the strands of β-sheet floor. the binding is mostly dominated by hydrophobic interactions with the hydrocarbon chain of the ligand; a single h-bond is established between the polar group of phospholipid and arg employing the pipeline described above, we have carried out an'in silico' saturation mutagenesis of the epcr residues involved in lipid binding. each modelled mutant was subjected to molecular docking simulation and stability evaluation. the protein-ligand complexes obtained by docking simulations were screened by their inhibition constant values in order to estimate the impairment of ligand binding ability induced by the mutation. the overall procedure and a summary of the results are depicted in figure the interacting residues, identified using ligplot, were saturated with the remaining amino acids, producing mutant models. then, wild type and mutants were analysed by docking simulations. the wild type re-docking correctly reproduced the x-ray ligand conformation with a rmsd of Å, a binding energy of - kcal/mol and a ki value of nm. positions and substitutions critical for binding affinity were identified with a ki-based filter. referring to the wild type ki value, we used a threshold of μm to discriminate between mutants with probable impaired binding ability. the pipeline selected mutants with ki value higher than the threshold, substituted in out of the starting positions. pair combinations of the substitutions were generated and double mutants were modelled. they were analyzed by the same docking protocol of the single mutants and filtered with a ki threshold of μm. this value was set an order of magnitude higher than for single mutants, with the purpose of selecting those combinations with a binding ability worse than the single ones. with this filtering criterion, double mutants were rescued. the large reduction in the number of double mutants, suggested that most residue combinations did not cooperate to increase the ki value compared to the single substitutions. the structural stability of single and double variants that passed the corresponding ki thresholds was evaluated using the foldx tool, directly integrated in the pipeline. mutants ensuing a ΔΔg value higher than kcal/mol were considered structurally unstable and removed from the selection. after the screening, the pipeline retrieved single and double mutants with a ΔΔg < kcal/mol. binding energy and ki values intervals are reported in figure the stability filter was passed only by mutants with a ki value within and μm for single and double, respectively. this value was an order of magnitude higher than the respective ki thresholds. mutants with high ki values, although potential candidates for binding affinity studies, were rejected by the ΔΔg filter and thus considered unstable. at the end of this analysis, key positions for the ligand binding in epcr were identified: position and are localized in the h α-helix, positions and are localized in h α-helix and one, the position on the eight-stranded β-sheet floor. the most represented substitutions of these key positions involved hydrophobic residues, such as w, in positions and and l, in positions and pair combinations in double mutants were evaluated in terms of synergistic effect where the ki value of double mutant was larger than the sum of the corresponding single ones. to estimate the synergistic effect, we applied the following equation: where ki and ki represent the ki value of the first and the second single mutant, and ki1+ the ki value of double mutant. kid is the difference between the effect of the double mutation and the sum of the two corresponding single mutations. results are shown in figure we found that of the double mutants displayed a partial additive effect as the ki value was included between the more damaging mutation and the sum of the two corresponding single mutants. the remaining double mutants had ki values higher than the sum of the two corresponding single mutants, to which we assigned a synergistic effect compared to the single mutation. we observed that the highest ki values were associated with double mutants resulting from replacement of a with either i or l and f with a small polar/apolar side chain. between them, the highest synergistic effect was shown by a31i/f164s mutant, while the largest effect in terms of binding energy and ki value was associated to a31i/f164n mutant. in conclusion, using the described approach, we have selected a limited number of single and double mutants of epcr for which we predict a reduced binding ability for the pty ligand and a structural stability comparable to the wild type. the set of selected single mutants and the double mutants presenting synergistic properties can be proposed as rationally constructed candidates for experimental testing in cellular models. pipeline performance due to the computational load needed to accomplish all the simulations, a parallel infrastructure was employed. in this test case, the simulation of the single and double mutant took on average and hours, respectively. the whole challenge consisted of docking experiments for the single mutants plus the wild type simulation and for the double mutants. depending on the resource availability, it was possible to employ cpus on average. therefore, in a period of about ten days we collected all the results covering about years of computer activity.  we present a high performance semi-automated computational pipeline that integrates modelling and docking simulation algorithms for large scale mutant design. additional analysis tools such as structure stability evaluation are also included. the pipeline is organized into two main blocks for dealing with single and double mutants that are interconnected so that selected results of first analysis are used for the second part of the study. the procedure is presently devoted to protein-ligand interaction studies, but can be adapted to investigations of protein-protein interaction. the mutagenesis method is based on a saturation approach and does not require previous knowledge of functional or structural role of involved residues: it can therefore be applied to explore new binding features. the implementation of the pipeline on a parallel computer infrastructure permits a greatly reduced computational time and makes it possible to test a very large number of mutants. for example, in the reported text case we have modelled and tested variants of epcr with the purpose of identifying single or coupled amino acid mutations that significantly impair the binding of the lipid ligand. the whole procedure reduced the number of selected models to which corresponds to about 2% of the initial dataset. selected variants provide information about the probable functional and structural effect of the mutated residues and can be proposed as subjects for experimental tests. in addition, the large amount of data stored in the database after each step of computational analysis can be further retrieved for more detailed examination and study. this pipeline can be useful for protein-ligand binding studies to identify key residues for the ligand interaction, to evaluate the effect of the substitutions on different ligands and to design selective mutants. competing interests the authors declare that they have no competing interests. authors' contributions fc developed the core of the pipeline, tested all the software for modelling and docking simulations, prepared the biological dataset, analyzed the output results and drafted the manuscript. pd conceived and designed the study and critically revised the results. im contributed to developing the pipeline, designing the database and enabling the parallel computation of the docking experiments on the cluster infrastructure. lm supervised the project and made available the computational resources. er coordinated the study and critically revised the results and manuscript. all authors read and approved the final manuscript. 
for an exhaustive analysis of selected positions of protein structure large scale mutagenesis techniques are often employed, with laborious and time consuming experimental set-up. each modelled mutant was subjected to molecular docking simulation and stability evaluation. mutants are screened by docking simulation and stability evaluation followed by a rationally driven selection of those presenting the required characteristics. the computational pipeline presented in this work is a useful tool for the design of structurally stable mutants with altered affinity for ligand binding, considerably reducing the number of mutants to be experimentally tested. the saturation mutagenesis procedure does not require previous knowledge of functional role of the residues involved and allows extensive exploration of all possible substitutions and their pairwise combinations. the mutagenesis was performed with a'saturation' approach, where all natural amino acids were tested in positions involved in ligand binding sites. the simulated protein-ligand complexes were screened for their impairment of binding ability based on change of calculated ki compared to the wild-type.'in silico' mutagenesis and screening simulation represents a valid alternative to laboratory methods to drive the'in vivo' testing toward more focused objectives. we present here a high performance computational procedure for large-scale mutant modelling and subsequent evaluation of the effect on ligand binding affinity.
 the design of mutants in protein functional regions, such as the ligand binding sites, is a powerful approach to recognize the determinants of specific protein activities in cellular pathways. for an exhaustive analysis of selected positions of protein structure large scale mutagenesis techniques are often employed, with laborious and time consuming experimental set-up.'in silico' mutagenesis and screening simulation represents a valid alternative to laboratory methods to drive the'in vivo' testing toward more focused objectives.  we present here a high performance computational procedure for large-scale mutant modelling and subsequent evaluation of the effect on ligand binding affinity. the mutagenesis was performed with a'saturation' approach, where all natural amino acids were tested in positions involved in ligand binding sites. each modelled mutant was subjected to molecular docking simulation and stability evaluation. the simulated protein-ligand complexes were screened for their impairment of binding ability based on change of calculated ki compared to the wild-type. an example of application to the endothelial protein c receptor residues involved in lipid binding is reported.  the computational pipeline presented in this work is a useful tool for the design of structurally stable mutants with altered affinity for ligand binding, considerably reducing the number of mutants to be experimentally tested. the saturation mutagenesis procedure does not require previous knowledge of functional role of the residues involved and allows extensive exploration of all possible substitutions and their pairwise combinations. mutants are screened by docking simulation and stability evaluation followed by a rationally driven selection of those presenting the required characteristics. the method can be employed in molecular recognition studies and as a preliminary approach to select models for experimental testing. 19– may bioinformatics methods for biomedical complex systems applications varenna, italy 