A proximity-based in silico approach to identify redox-labile disulfide bonds: The example of FVIII

Allosteric disulfide bonds permit highly responsive, transient ‘switch-like’ properties that are ideal for processes like coagulation and inflammation that require rapid and localised responses to damage or injury. Haemophilia A (HA) is a rare bleeding disorder managed with exogenous coagulation factor(F) VIII products. FVIII has eight disulfide bonds and is known to be redox labile, but it is not known how reduction/oxidation affects the structure-function relationship, or its immunogenicity—a serious complication for 30% severe HA patients. Understanding how redox-mediated changes influence FVIII can inform molecular engineering strategies aimed at improving activity and stability, and reducing immunogenicity. FVIII is a challenging molecule to work with owing to its poor expression and instability so, in a proof-of-concept study, we used molecular dynamics (MD) to identify which disulfide bonds were most likely to be reduced and how this would affect structure/function; results were then experimentally verified. MD identified Cys1899-Cys1903 disulfide as the most likely to undergo reduction based on energy and proximity criteria. Further MD suggested this reduction led to a more open conformation. Here we present our findings and highlight the value of MD approaches.


Introduction
Disulfide bonds are covalent bonds formed between the sulfur atoms of two cysteine residues in a protein. These bonds contribute to the mechanical stability of the protein tertiary structure and this stability can allow certain proteins to survive and function in harsh extracellular environments that can be rich in proteases and subject to large changes in salt and pH that would otherwise degrade a non-covalently stabilised protein structure [1]. However, it is becoming increasingly apparent that a subset of disulfide bonds are more susceptible to reduction than others and that these labile disulfide bonds can, when reduced, impart structural changes that translate to modulated function. These functional disulfide bonds are termed allosteric disulfide bonds and they permit rapid and reversible changes in protein structure that negates the need for energetically expensive and time consuming de novo synthesis [2] and are often referred to as "redox switches". There are many examples of functional allosteric disulfide successful in identifying labile disulfide bonds with allosteric properties in many proteins, however these can be complex and time consuming processes. In this study we describe a new predictive modelling method for identifying labile disulfide bonds within protein structures and validate it using known labile disulfide bonds reported in the literature. This approach is based on all-atom molecular dynamics (MD) simulations of proteins in explicit solvent and in the presence of tris(2-carboxyethyl)phosphine (TCEP), a reducing agent frequently used to break disulfide bonds. Our approach was then applied to predict which disulfide bonds are labile in coagulation FVIII. These were then experimentally validated and their effect on FVIII activity quantified.
All the simulations described in the following were performed using Gromacs (5.1.4 or 2018.6) [37] and Plumed (2.4.1 or 2.5.1) [38], employing the GROMOS54A7 force field [39] for the protein and the SPC/E force field [40] for water. The topology file for TCEP was obtained from the automated topology builder (ATB) server [41], and a 100:1 TCEP to protein ratio was used in all simulation boxes.
For all the simulations described in the following we used periodic boundary conditions, and the cut-off radius for both Coulombic (calculated using the PME method [42]) and Lennard-Jones interactions was set to 1.2 nm. In all simulations, 1 native protein molecule was introduced into the simulation box, and the overall charge was neutralized using Na + or Clions. Each box was then energy minimized using the steepest descent algorithm, and equilibrated for 1 ns in the NPT ensemble, using Berendsen pressure and temperature coupling [43]. The production runs were also performed in the NPT ensemble, controlling temperature with the Nose-Hoover thermostat (0.5 ps relaxation time) [44,45], and pressure with the Parrinello-Rahman barostat (3 ps relaxation time) [46]. The LINCS algorithm was employed for constraining all bonds [47].
First, 3-replica multiple walker parallel bias metadynamics (PBMetaD [48]) was used to select a good proximity-based criterion for the identification of labile disulfides. Afterwards, multiple unbiased MD simulations exploiting this criterion were performed.

Parallel bias metadynamics
Metadynamics [49] works by introducing a history-dependent bias potential V(s i ,t) that acts on selected degrees of freedom of the system s i , generally referred to as collective variables (CVs), where ω i is a deposition hill height and σ i a Gaussian width. Such bias potential pushes the simulated system out of local minima, promoting the exploration of a considerably larger fraction of the phase space compared to conventional MD simulations. The computational cost associated with the deposition of the bias, however, increases exponentially with the number of CVs selected. Parallel bias metadynamics alleviates this issue by constructing multiple one-dimensional biases, each acting on a single CV in parallel. This makes it possible to include as many CVs as needed at a reasonable computational cost.
3-replica multiple walker PBMetaD simulations were performed at 310 K and 1 bar, for CD44 in presence of 100 TCEP molecules. The simulation box was cubic with 7.9 nm side length. Seventy nanoseconds per replica (210 ns total simulation time) were performed, and a 1 fs time step was used.
The protein radius of gyration (R g ), and the coordination number of the TCEP molecules (CN) around the protein were used as CVs. The Gaussian height was set to 2 kJ/mol, the bias factor to 15, and the Gaussian deposition rate to 1 hill/ps. The σ (Gaussian width) values used were 0.02 nm and 2 for R g and CN, respectively.
The radius of gyration for a molecular structure containing N atoms is defined as, ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where R com is the position of the centre of mass (COM), and m i and R i are the mass and position of atom i, respectively. This CV was used to improve the sampling of more expanded, but not unfolded, configurations of CD44. The bias factor and Gaussian height were in fact defined so as to impede complete unfolding. The coordination number CN was defined as, where r ij is the distance between species i in group A (containing N A entities) and j in group B (including N B entities), and r 0 is a reference distance. Group A was the protein COM, while group B comprised the centre of masses of all TCEP molecules. The value of r 0 was set to 3 nm, while the exponentials n and m were set to 6 and 12, respectively. This preliminary 3-replica multiple walker PBMetaD simulation allowed us to select a proximity-based criterion for the identification of labile disulfide bonds. The selected parameter was the coordination number of TCEP molecules at 0.8 nm from the disulfide bond. This criterion was subsequently used in unbiased MD simulations.

Unbiased MD simulations to identify labile disulfide bonds
The labile disulfides in some of the model proteins used have been experimentally validated (Cys77-Cys97 for CD44 [10], Cys160-Cys209 for CD132 [11], Cys3-Cys127 for IL-4 [11], and Cys186-Cys209 for tissue factor [7]). These proteins were preliminarily simulated in order to access and refine the accuracy of the simulation approach used. The same simulation approach was then extended to the light chains of FVIII. Snapshots of the model proteins used in this work, with highlighted disulfides, are shown in Fig 1. For each of these proteins, we started at least 20 independent trajectories. Each protein was simulated in presence of 100 TCEP molecules, at 310 K and 1 bar, and using a 1 fs timestep. Simulation boxes were cubic, with side length equal to 10.7 nm for tissue factor, 10.2 nm for CD132, 7.9 nm for IL-4 and CD44, 10.6 nm for GP130 and 14.4 nm for FVIII. CD44 and CD132 were also simulated in presence of dithiothreitol (DTT), again at a 100:1 DTT:protein ratio. DTT is another common reducing agent, often used to break disulfide bonds in proteins. The topology file for DTT was also downloaded from the ATB server. The aggregated simulation time for these unbiased simulations was about 2.29 μs.
During these unbiased trajectories, we tracked the position of the TCEP (or DTT) molecules. We then applied two different criteria to identify a labile DSB [50], (1) distance-based or (2) distance+energy-based. According to the first criterion, a disulfide bond was considered to be reduced as soon as a TCEP molecule got closer than 0.8 nm. This same proximity-based condition had to be met for criterion number 2, but in this case a further requirement had to be satisfied. Specifically, the simulation was stopped, the disulfide bond was reduced and the resulting system energy minimized using the steepest descent algorithm. The initial structure, before disulfide bond reduction, was also energy minimized. The potential energies of the two systems, before (E in ) and after (E fin ) disulfide reduction, were compared and the disulfide deemed to be broken only if E fin <E in .

Analysis of FVIII conformational changes upon reduction of the labile disulfide bond
After the identification of the FVIII labile disulfide bond in our trajectories, two additional simulations were performed for B-domain deleted FVIII (PDB ID: 3CDZ [36]) in water at 300 K and 1 bar, to investigate how reduction of the labile disulfide bond would affect protein conformation. In this case, both the heavy and the light chains of the protein were simulated.
The force fields used were the same previously described, but in this case no TCEP molecules were inserted into the simulation box. For one simulation, the labile disulfide bond was reduced before starting the trajectory, while native FVIII was used in the other case. The protein molecules were inserted into two separate cubic boxes with 14.9 nm side length. The overall charge of the systems was neutralized by the addition of Clions. The protein structures were simulated for 200 ns at 300 K and 1 bar in the NPT ensemble, using a 2 fs time step.
The protein radius of gyration (R g ), number of intra-molecular hydrogen bonds (HB pp ), RMSD (root mean square deviation) with respect to the initial structure and solvent accessible surface area (SASA) as function of simulation time were computed. The conformations assumed by FVIII were grouped together by performing a cluster analysis based on the Daura algorithm [51], after discarding the first 100 ns of the trajectory. More specifically, the conformations were grouped together if the root mean square deviations (RMSDs) of the protein backbone (N-Cα-C atoms) were less than 0.2 nm compared to each other. The most sampled protein conformation during the equilibrated trajectories (last 100 ns) was selected based on this analysis for the case of both the reduced and oxidized forms of FVIII. The change in secondary structure as a result of disulfide reduction was also investigated, by comparing the equilibrated structures obtained from the two simulations.

Kinetic trapping of labile disulfide bonds in recombinant FVIII
FVIII products (diluted to 200 international units [IU] per ml) were incubated either with 1mM TCEP at room temperature (RT) for 15 minutes, or with Trx (1μM), NADPH (0.2 mM), and thioredoxin reductase (TrxR1, 1 nM) at 37˚C for 30 minutes. Alexa Fluor 488-maleimide was then added to label free thiols (30 min, room temperature, in the dark). Samples were then boiled in Lamelli buffer and separated on a 4-20% gel by SDS-PAGE and visualised using a GeneSys Imager. Total protein was evaluated using Coommassie blue.

Quantitative mass spectrometry analysis of disulfide bonds in FVIII
50μg of recombinant FVIII was added to a number of 10kDa 500 μl centrifugal concentrators (Vivacon 500, Sartorious) and reduced with one of 1mM TCEP in PBS, 1mM DTT in PBS and Trx (1μM), NADPH (0.2 mM), TrxR1 (1 nM) for 90 min at 37˚C. After washing with 500 μl PBS, free cysteines arising from the reduction of disulfide bonds were alkylated with iodoacetamide (IAA) (100 μl of 50mM in 25mM ammonium bicarbonate for 1 h at room temperature in the dark). A non-reduced control sample was treated with IAA only to obtain the base level of free cysteines in FVIII. The samples were denatured and any remaining disulphide bonds TCEP-reduced (200 μl of 8 M urea solution with 25 mM TCEP for 1 h at room temperature) after which the tubes were centrifuged to remove liquid from filters. Cysteines arising from this were alkylated with NEM (200 μl of 50mM in PBS for 60 minutes at room temperature) to discriminate them from cysteines arising from the initial 1μm TCEP reduction of the intact FVIII. A 100% IAA labelled sample of FVIII was also prepared by denaturing and reducing 50μg in a centrifugal concentrator and alkylating with IAA as above. Three biological replicates were prepared and the samples were trypsin-digested as reported previously [52]. Briefly, samples were incubated with 2 μg trypsin (Sigma) in 100 μl Ambic at 37˚C overnight. Any remaining liquid was then spun into a clean collection tube, and 100 μl 0.1% formic acid was added to each sample for 10 min before centrifuging the liquid into the collection tube. Solutions of 50% and 100% acetonitrile with 0.1% formic acid were then consecutively added to the samples, which were centrifuged to remove liquid at each stage. The filters were then removed, and samples placed in a speedvac (Thermo) at 45˚C to remove all liquid. For mass spectrometry, samples were reconstituted in 0.1% formic acid, followed by sonication.

Mass spectrometry
The tryptic peptide samples were injected as technical triplicates onto an Ultimate 3000 nano HPLC system coupled to an Orbitrap XL Discovery mass spectrometer (Thermo Scientific).
Samples were online desalted on a μ-Precolumn (C18 PepMap100, 300 μm id × 5 mm; 5 μm, 100 Å) at a flow rate of 25 μL/min, which was followed by separation on a nano analytical column (Acclaim PepMap100 C18, 75 μm id × 50 cm, 3 μm, 100 Å) (Thermo Scientific) using a 90-minute linear gradient from 5 to 40% solvent B (98% CH 3 CN/2% H 2 O/0.1% formic acid, v/v/v) versus solvent A (98% H 2 O/2% CH 3 CN/0.1% formic acid, v/v/v) at a flow rate of 300 nl/min. The mass spectrometer was operated in a data-dependent acquisition mode. The full survey scan (m/ z 400-2000) was acquired in the Orbitrap with a resolution of 30,000 at m/z 400, which was followed by five MS/MS scans in which the most abundant peptide precursor ions detected in the preceding survey scan were dynamically selected and subjected for collision-induced dissociation (CID) in the LTQ (linear ion trap) to generate MS/MS spectra ('top-5 method').

Quantitative analysis of disulfide bond REDOX state after reduction with TCEP
Data files of the mass spectrometry runs were combined and searched against the human Swiss-Prot database using Peaks 8 proteomics studio (Bioinformatics Solutions Inc. On, Canada). Precursor mass tolerance was 10 ppm and fragment ion tolerance was 0.6 Da with up to two missed trypsin cleavage sites per peptide allowed. Variable modifications were defined as deamidation on asparagine and glutamine, oxidation on methionine and alkylation with IAA or NEM (and hydrolysed variants) on cysteines. de-novo, peaks-db, SPIDER and peaks PTM algorithms were sequentially used to search against a concatenated target/decoy database, providing an empirical false discovery rate (FDR) and results are reported at a 1% target/decoy FDR for both peptides and proteins.
Peptides containing IAA labelled cysteines were selected for further analysis ensuring at least one cysteine from each disulfide bond was covered. Four control peptides from FVII, all of which contained no modifiable residues were selected and used for normalisation. All of the peptides used in the analysis are summarised in S1 Table. Precursor ion areas for the selected peptides were extracted using MS1 filtering in Skyline [53]. To normalise the data the ratio of precursor ion area of the peptide under analysis to the precursor ion areas of each of the control peptides is found for the control, Trx and 100% reduced samples. The % reduction for a given cysteine is calculated using the following equation for each peptide in each of the biological and technical replicates: n From this, the mean % reduction and standard deviation for each cysteine before and after reduction with DTT, TCEP and Trx was determined.

FVIII activity assays
Chromogenic assay. FVIII activity was quantified using the Chromogenix Coatest 1 SP4 FVIII assay kit according to the manufacturers' instructions on an ACL Top 550 blood analyser (Werfen). In brief, FVIII was diluted to~1 IU/ml in FVIII-deficient plasma and to this, activated FIXa, FX, Ca 2+ , phospholipid, and the FXa substrate S-2765 were added to the appropriate concentrations and conversion of the substrate to a coloured product was measured at 405 nm.

Enzyme-linked immunosorbent assay (ELISA)
Nunc Maxisorb plates were coated with 500 ng recombinant B domain deleted (BDD) FVIII for 2 hours at room temperature. Wells were then incubated with either PBS or 1 mM TCEP for 30 minutes at room temperature before being blocked with PBS/1% bovine serum albumin (BSA)/0.05% Tween 20 for 30 minutes (at room temperature). Recombinant anti-FVIII antibodies were added to the wells (100 μl, diluted in blocking buffer) and incubated for 1 hour at room temperature with gentle agitation. Unbound antibody was removed by washing (3 x 200 μl PBS/0.2% Tween 20) and anti-FVIII antibodies detected using an anti-human k chain antibody (Sigma Aldrich #A7164) conjugated to horse radish peroxidase (1:20,000) and 3,3 0 ,5,5 0 -Tetramethylbenzidine (TMB) substrate. Plates were read on a microplate reader (SpectraMax M2e, Molecular Devices) as an end point assay.

Statistics
Chromogenic assay and one-stage clotting assay (OSCA) data were analysed using a parallel line model (with a log transformation for the OSCA) and evaluated for linearity and parallelism using Combistats Version60. Data represent three experimental replicates.

Haemostasis-related plasma proteins are substrates of Trx and contain labile disulfide bonds
A major driving force for this work is to understand the redox lability of coagulation factor (F) VIII. The redox lability of several haemostasis-related plasma proteins has been experimentally validated and we began by using these to validate our reduction protocols. Samples were incubated under chemical or enzymic reducing conditions before free cysteines were labelled with a thiol-reactive maleimide probe conjugated to a fluorescent dye. Samples were then separated by SDS-PAGE and visualised. As shown in the upper panel of Fig 2A, both enzymic and chemical reduction results in the appearance of bands corresponding to the proteins of interest, or an increase in their intensity, indicative of an increase in free cysteines; this confirmed we were able to reduce labile disulfides using both TCEP and thioredoxin and detect cysteine labelling in proteins previously reported to contain labile disulfide bonds, including factor (F) XI, von Willebrand factor (VWF), and beta glycoprotein I [13,21,[54][55][56][57]. Furthermore, in agreement with previous reports [26], we can also detect an increase in free cysteine labelling when recombinant B domain deleted (BDD) FVIII (BDD rFVIII) is incubated under chemical or enzymic reducing conditions (Fig 2B). The light chain exhibits more intense staining than the heavy chain, which is noteworthy because the light chain conveys the majority of FVIII functionality (FIXa binding, VWF binding, and the interaction with phospholipids). We next set out to develop an in silico method to identify which of the 4 disulfide bonds in the FVIII light chain are the most likely to be redox labile using a molecular dynamics (MD) approach.

Identification of suitable proximity criteria to predict disulfide bond redox susceptibility
Using a simple proximity criterion previously used to investigate thiol-disulfide isomerization in a mutated immunoglobulin [50], we looked at the interaction of TCEP (tris(2-carboxyethyl) phosphine) with proteins with experimentally validated redox-labile disulfide bonds (Fig 1) before applying the method to the FVIII light chain. We first applied the parallel bias metadynamics method, which allows an extensive sampling of configuration space, to select a good proximity-based criterion for the identification of labile disulfides in proteins. We started with the Cys77-Cys97 disulfide bond in CD44 looking at three parameters-1) the average distance between TCEP and the disulfide bonds, 2) the minimum distance between TCEP and the disulfide bonds, and 3) the average number of TCEP molecules whose centre of mass (COM) was closer than 0.8 nm to the COM of the disulfides. This 0.8 nm cut-off was selected as it is both small enough to capture local phenomena and large enough to accommodate at least a few TCEP molecules.
The average distance between TCEP molecules and disulfide bonds (parameter 1) is lowest for Cys53-Cys118 (Fig 3A), but this disulfide bond is not reduced under experimental conditions indicating that this parameter is not suitable for predicting labile disulfide bonds. The reason for this becomes evident when considering the data in Fig 3B (parameter 2), which reveals that although TCEP can get close to the disulfide bond, TCEP molecules cannot get sufficiently close to react with it. The average distance parameter alone (Fig 3A) is not a reliable predictor for disulfide bond lability, and although the minimum distance parameter (Fig 3B) can identify labile disulfides in the framework of a proximity criterion, it is not, at least in this case, sensitive enough to distinguish between Cys77-Cys97, and Cys28-Cys129. The third criterion-the number of TCEP molecules that, on average, were closer than 0.8 nm to the COM of the disulfide bonds (Fig 3C)-allowed clearer separation of disulfide bonds with increased sensitivity compared with the minimum distance parameter. In our MD simulations, the CD44 disulfide Cys77-Cys97 is surrounded by the highest number of TCEP molecules (a proxy for concentration) when a 0.8 nm sphere around the COM is defined. This result is in line with experimental data for CD44 that identified Cys77-Cys97 as redox-labile [10].
However, the reduction of a disulfide bond is the result of both accessibility (i.e., TCEP molecules need to get close enough to react) and net energy gain (i.e., the reduced system should be energetically favoured). A set of unbiased simulations was therefore started for each protein shown in Fig 1, using either a (1) distance-or (2) a distance+energy-based criterion, as detailed in the Materials and Methods section. According to the distance-based criterion, the disulfide bond was considered to be reduced as soon as a TCEP molecule got closer than 0.8 nm. In contrast, a second requirement had to be satisfied according to criterion (2)-the potential energy of the reduced (E fin ) system had to be smaller than that of the oxidized (E in ) one (E fin <E in ). The results of these unbiased simulations are shown in Fig 4. The average solvent accessible surface area (SASA) and root mean square deviation (RMSD) of each disulfideforming cysteine during the unbiased simulations are displayed in S1 and S2 Figs.
The two criteria (distance or distance+energy) were applied to 5 experimentally validated redox-labile disulfide bonds including CD44 (Fig 4A), TF (Fig 4B), CD132 (Fig 4C), IL-4 ( Fig  4D), and the D1 and D2 domains of GP130 (Fig 4E). The distance-only criterion (plain black bars) was enough to identify unambiguously the labile disulfides Cys77-Cys97 in CD44 ( Fig  4A), Cys186-Cys209 in TF (Fig 4B), Cys160-Cys209 in CD132 (Fig 4C), Cys3-Cys127 in IL4 (Fig 4D), and Cys6-Cys32 in GP130 (Fig 4E). In this case, adding the energy-based requirement (dashed bars) did not modify the main conclusions. For instance, the distance-only criterion predicted that Cys77-Cys97 in CD44 would have 96.7% probability to be reduced, and this probability further increased to 100% when the distance+energy-based criterion was applied. The average SASA of the different disulfide bonds correlates fairly well with their lability (see S1 Fig) and is therefore a good initial approximation. However, it is important to note that the distance-based criterion is not only a measure of solvent exposure, but also takes into account the specific interactions between the probe (TCEP) and the protein surface. For instance, in our simulations of CD44 with TCEP, the solvent accessibility of Cys77-Cys97 was similar to that of Cys28-Cys129, or even lower, as shown in S1A Fig. However, TCEP molecules mostly remained closer to Cys77-Cys97, presumably because of specific interactions with this patch on the protein surface. Moreover, it is important to note that the solvent accessibility, as measured in MD simulations, is representative of protein dynamics, as such being more reliable than the solvent accessibility measured from static protein structures.
We also studied the interactions between dithiothreitol (DTT) and CD44 or CD132 (S3 Fig). DTT is smaller than TCEP (the radius of gyration of DTT is �0.26 nm, while it is �0.37 nm for TCEP). For this reason, a smaller cut-off (0.6 nm, instead of 0.8 nm) was used to estimate proximity of the DTT molecules. We found that the distance-based criterion could not distinguish between Cys77-Cys97 and Cys28-Cys129 in CD44 (S3A Fig), while adding the energy-based requirement succeeded in doing so, identifying Cys77-Cys97 as redox labile with 60% probability. The difference between S3A Fig (CD44-DTT interaction) and 4A (CD44-TCEP interaction) indicates that the selected probe, DTT or TCEP, plays a role (as different probes preferentially interact with different patches on the protein surface), although the ultimate conclusions on redox lability are not affected, when the distance+energy-based criterion is used. S3B Fig further shows that both the distance and distance+energy-based criteria succeeded in identifying Cys160-Cys209 in CD132 as redox labile also in the case of DTT as molecular probe.
We next applied our in silico method to FVIII. The relatively large size of FVIII presents a challenge for in silico modelling with regards to computational costs, so we ran simulations for the light chains only (Fig 4F) in this instance. FVIII light chains convey most of the functionality and we had already identified this part of the molecule as being of interest due to the considerable increase in free thiols following thioredoxin or TCEP treatment (Fig 2B). Modelling of the FVIII light chains reveals that a larger number of TCEP molecules can get close to Cys1899-Cys1903 than to Cys1832-Cys1858, Cys2021-Cys2169 or Cys2174-Cys2326, therefore suggesting that Cys1899-Cys1903 is the most likely to be redox labile (70.0% probability according to the distance-based criterion, plain black bars). Applying the energy-based condition further augments the selectivity, indicating Cys1899-Cys1903 as redox-labile with an 81.8% degree of confidence. Finally, it is also interesting to note, looking at S2 Fig, that there is no evident correlation between the lability of a disulfide bond and its RMSD during the MD trajectories.

Experimental validation of in silico predicted redox-labile Cys1899-Cys1903
Where suitable structures are available, some prediction can be made with regards to whether a disulfide bond may be redox labile or not. Using structural parameters like solvent accessibility or the distance between the alpha-carbons of the constituent cysteines and their 3D conformation [56], one can get an indication of how redox labile a given disulfide bond may be. Allosteric disulfide bonds tend to have high torsional bond strain as depicted by shorter alpha carbon distances; the −RHStaple and the −/+RHHook atomic arrangements are particularly prevalent conformations for allosteric disulfide bonds and more likely to undergo cleavage [57]. Using an online tool [58] we assessed redox labile disulfide bonds in published FVIII structures and their structural parameters are shown in Table 1. There is a degree of variability between FVIII structures due to the facts that crystallographic structures are snapshots of a flexible molecule-proteins are dynamic structures that have intrinsic flexibility to permits movement. The light chain disulfide bond in the A3 domain (Cys1899-Cys1903) predicted to be labile by our MD simulation is absent from three of the four structures and shown as reduced. As labile disulfide bonds are often broken by the reducing environment in the x-ray beam during data collection the lack of this disulfide bonds in the three structures lends some support to the MD prediction. Although only one set of parameters are available for this bond, when these are compared to the other bonds (average parameters over the four structures) it is the only one with high (>30 A 2 ) solvent accessibility and a short (<5 Å) distance between its alpha carbons. The other bonds are either poorly solvent exposed or have a long alpha carbon distance.
We applied a quantitative mass spectrometry and kinetic trapping of cysteine redox state method previously used to quantitate labile disulfide bonds in monoclonal antibodies [50] to examine the redox state of each disulfide bond in FVIII treated with chemical or enzymic reducing agents (Fig 5). The redox state of the disulfide bonds was determined from peptides containing at least one cysteine from each disulfide bond (details shown in S1 Table) except for the Cys630-Cys711 bond as there were no suitable tryptic peptides to monitor (the cysteine numbering from the mass spectrometry includes the 19 amino acid signal peptide so the numbers are 19 higher than the numbering in the parentheses which is consistent with the rest of the data). None of the disulfide bonds showed reduction in the control FVIII with the exception of Cys1899-Cys1903, which was approximately 20% reduced. However, following DTT, TCEP and Trx treatment, the Cys1899-Cys1903 disulfide bond was further reduced to~100% Table 1. Structural parameters of disulfide bonds using online prediction analysis. Four published FVIII structures were evaluated in parallel using the UNSW DSB prediction tool [58]. Allosteric disulfide bonds tend to be solvent accessible and have short alpha carbon distances. Stable structural disulphide bonds tend to have a -LHspiral conformation; these have the most common DSB conformation and have the lowest dihedral strain energy.

Solvent Accessibility Å2
Cα indicating it is labile and supporting the results obtained from the MD simulations. Interestingly, Trx treatment also led to the reduction of Cys248-Cys329, which sits in the heavy chain of FVIII. This disulfide bond was also fully reduced revealing that it is labile towards enzymatic Trx reduction but not chemical (DTT or TCEP) reduction.

Effect of Cys1899-Cys1903 reduction on FVIII structure
The data shown in Fig 5 supports the FVIII light chain in silico modelling that predicted Cys1899-Cys1903 is redox labile. We next returned to in silico approaches to try and shed light on how reduction of Cys1899-Cys1903 disulfide bond could affect FVIII structure and function, this time running simulations for the heavy and light chains of FVIII in its oxidized and reduced (broken Cys1899-Cys1903 bond) forms. It is important to note that enhanced sampling methods were impracticable for these simulations due to the very large size of FVIII, and we therefore conducted conventional MD runs. We extended the simulation time until the

PLOS ONE
computed properties (radius of gyration, internal hydrogen bonding network, solvent accessible surface area and backbone RMSD compared to the crystal structure, as shown in Fig 6) converged to stable values. The results of the conventional MD simulations presented are therefore a good indication of what may happen immediately after reduction of the labile disulfide bond; however, because of the impossibility to use enhanced sampling techniques, the reader should be aware that they may be not completely representative of long-time shifts in the protein conformation. The radius of gyration R g decreased for both the oxidized and reduced forms compared to the crystal (static) structure (Fig 6A), with the reduced molecule having a slightly more expanded structure (R g � 3.58 nm compared to R g � 3.53 nm for the oxidized conformation) by the end of the simulation. This corresponded with a modest decrease in intra-molecular hydrogen bonds HB pp (� 938 for the oxidized conformation and � 915 when the disulfide is reduced) and an increase in the solvent-accessible surface area (SASA) (� 519 nm 2 for the oxidized protein compared to � 531 nm 2 with the reduced molecule).
The most sampled protein conformations (according to the Daura algorithm{Daura, 1999 #809}) for the equilibrated trajectories (last 100 ns) are shown in Fig 6B. Fig 6C shows an alignment of the different domains of FVIII in the oxidized and reduced forms. The protein regions Evolution of FVIII radius of gyration R g , number of intra-molecular hydrogen bonds HB pp , solvent accessible surface area SASA, and backbone RMSD with respect to the crystal structure, as function of simulation time, for both fulloxidized FVIII (black curve) and reduced FVIII (red curve). (B-C) Cartoon representation and superimposition of fulloxidized and reduced (Cys1899-Cys1903) FVIII. The alignment has been performed either on the whole structure (B), or on the separate domains (C). The cartoons are coloured according to the Q parameter, which is 1 for identical structures; the blue areas indicate structural similarity between full-oxidized and reduced FVIII, where Q = 1. Areas of poor structural similarity (Q = 0.1-0.3) are shown in red. The figure was obtained using the VMD software [59] and the STAMP algorithm [60], that works by minimizing the C α distance between aligned residues of each molecule. The Cys1899-Cys1903 disulfide bond is represented as yellow beads.
https://doi.org/10.1371/journal.pone.0262409.g006 subject to the largest conformational change upon reduction of the disulfide seemed to be residues 200-230 (A1 domain) and 555-570 (A2 domain) in the case of the heavy chain, and residues 1706-1732 (A3 domain), 1790-1804 (A3 domain) and 1880-1920 (A3 domain, this last region includes the reduced disulfide) in the light chain. Some slight conformational changes were also noted for the C1 domain (Fig 6C), while the C2 domain remained mostly unaltered. Using a panel of patient-derived FVIII neutralising antibodies developed in-house that bind to different regions of the FVIII molecule, we next evaluated FVIII-mAb binding in solid state using an ELISA. Fig 7 shows that, in support of the in silico modelling, reduction of recombinant full-length human FVIII (rFL-FVIII) decreases binding of the A2-and A3-specific anti-FVIII antibodies (NB11B2 and NB41, respectively). We also observe a decrease in the binding of the C1 domain-specific antibody NB33. No change in binding of the C2-specific antibody NB2C11 is observed. These data suggest that reduction of the Cys1899-Cys1903 disulfide bond results in a change in the conformation of FVIII through the A2, A3, and C1 domains.

Effect of DSB reduction on FVIII activity
We next evaluated the effect of reduction on the activity of a number of FVIII products, including recombinant B domain-deleted (BDD) porcine FVIII (rpFVIII), three different recombinant full-length human FVIII (rFL-FVIII), recombinant human (rBDD FVIII), and a human plasma-derived FVIII product (pdFVIII), which also contains von Willebrand factor (VWF). In agreement with the mass spectroscopy data, FVIII is redox labile with an increase in free thiols upon reduction (Fig 8). The effect of reduction on FVIII activity was next assessed using the chromogenic assay, which is used clinically to measure FVIII activity in patient plasma. FVIII is not an enzyme so its activity (function) is measured indirectly by quantifying the activity of the tenase complex (factor IXa and factor Xa) whose formation it catalyses. The chromogenic assay is comprised of separate kit components (factor IXa, factor X, Ca 2+ , phospholipids, and a tenase substrate that produces a coloured product) to which FVIII is addedproduct formation is proportional to the amount of 'active' (functional) FVIII in the sample. As shown in Fig 9, not all products behaved the same way when reduced. Recombinant porcine (BDD, Fig 9A) had marginally lower activity when reduced with DTT or TCEP, but reduction with Trx inhibited FVIII activity by~80%. In contrast, recombinant human BDD FVIII showed an increase in activity Fig 9B. We looked at three full length (FL) recombinant rFVIII products, two expressed in BHK call (Fig 9C and 9D) and one manufactured in CHO cells (E) and found that, again, there were differences in how they behaved upon reduction, with two products showing between 50-80% decrease in activity with Trx but only modest decreases in activity when exposed to chemical reducing agents (C and E). The increased effects seen with Trx over the chemical reducing agents can be attributed to the additional reduction of Cys248-Cys329 in the heavy chain (Fig 5). However, production in different cells lines gives rise to proteins with different glycosylation patterns and other post-translational modifications, and additional differences in production/manufacturing, processing, and excipients could all have an effect on the state of the final product, but it is not known how the data shown here relate to these processes. The plasma derived product was inhibited by reducing agents but there was little difference between enzymic and chemical ( Fig 9F). In summary, Trx has the capacity to reduce FVIII disulfides but the effect of these changes on FVIII activity varies from product to product (summarised in Fig 9G).

Discussion
Molecular dynamics is an extremely powerful tool for studying protein structure/function, especially for proteins that are difficult to work with due to their inherent instability. Molecular dynamics (MD) [59,60] can also provide a detailed picture of events often not directly measurable by existing experimental approaches. In the example here, we used MD to inform and guide experimental approaches to understand how reduction/oxidation affects the activity of coagulation factor(F) VIII. FVIII is an important biotherapeutic used to manage bleeding in haemophilia A (HA) patients but it is a labile protein and known to induce the formation of anti-drug antibodies in approximately one third of severe HA patients. It is known that FVIII is redox labile but there is some ambiguity about how reduction/oxidation affects this essential clotting factor (see below); identifying labile and allosteric disulfide bonds can help understand structure-function relationships and provide a basis for rational, hypothesis-driven protein engineering to improve stability and activity and reduce drug immunogenicity. We therefore applied a molecular dynamics approach to expand upon previously reported observations, resolve contradictory findings, and develop a hypothesis relating to how reduction of FVIII leads to changes in its activity. We applied a proximity criterion that correlates labile disulfide bonds with their accessibility, adding to previous studies that showed accessibility is a key factor [50]. To identify the proximity-based criterion, we made use of metadynamics [49], an enhanced sampling method that dramatically reduces the amount of computational time required to obtain a significant sampling of configuration space. We also showed the importance of adding an energy-based criterion for the correct identification of labile disulfides. This simple computational approach permits simulations for whole proteins at a detailed, atomistic level, with a reasonable computational cost. We then used classical MD simulations to develop hypotheses relating to how the change in protein conformation that results from reduction of the identified labile disulfide bond affects FVIII activity.
We found the concentration of TCEP at a distance of <0.8 nm from the centre of mass (COM) of the disulfide bond was a good predictor of disulfide bond lability, especially when combined with an energy-based condition, requiring that the disulfide reduction was energetically favourable. This distance+energy-based condition was able to successfully identify experimentally validated allosteric disulfide bonds in a number of proteins. Applying this parameter to the FVIII light chains indicated Cys1899-Cys1903 disulfide bond was the most likely to undergo reduction. This particular disulfide bond has been the subject of previous studies and it is also noteworthy that it is the only disulfide bond that is unique to FVIII compared to FV. Disruption of this bond (by mutation of Cys1903 to either glycine or serine) was reported to improve expression of recombinant B domain deleted (BDD) FVIII by 2-fold, with no impact on activity [61,62], which is a little surprising considering it resides in a part of the FVIII molecule that mediates FIXa binding [63]. Indeed, the engineering of disulfide bonds at the A2-A3 interface, despite having little impact on (predicted) domain architecture, did impart changes in FVIII-FIXa affinity and FX activation kinetics [64] suggesting this area is indeed important for the structure/function of FVIII. Early work [26] also highlighted a decrease in FVIII cofactor activity (FVIII:C) as well as a dissociation from its carrier protein von Willebrand Factor (VWF), another redox-labile protein that increases FVIII half-life in the circulation by~6-fold [65]. In contrast to this, infusion of N-acetylcysteine (NAC) was associated with a short-lived rise in measurable FVIII activity, although this increase in FVIII activity was only observed in patients with an adverse reaction to the NAC infusion and a resultant rise in VWF release from the endothelium [66,67]. There is also contradictory evidence in the literature regarding the effect of reduction on FIXa and the tenase complex (FIXa/FXa), whose formation is catalysed by FVIII. A more recent study investigating allosteric disulfide bond regulation of FIXa showed there was no change in FVIII or FIXa disulfide bond status when complexed, but did show a decrease in the low inherent activity of the tenase (FIXa/FX) complex when thioredoxin or PDI were present [67], again suggesting that reduction of FVIII supresses FVIII activity.
Our data lends weight to the argument that reduction of FVIII results in a decrease in activity and provides a structure-function rationale for this observation, i.e. reduction of the Cys1899-Cys1903 disulfide bond. Our modelling and experimental work indicates reduction of Cys1899-Cys1903 changes the structure of FVIII structure in a way that decreases its activity. In particular, the modelling suggests the largest structural changes are in the A1, A2, and A3 domains and this was experimentally validated with domain-specific antibodies. As FX and FIX bind these regions (C1 and C2 are predominantly for membrane association) it logically follows that the FVIII-FIX and FVIII-FX interactions are impacted by reduction-induced structural (allosteric) changes in FVIII. We were surprised to find that reduction of FVIII products resulted in both decreases and increases in activity measurements depending on which product was being evaluated. These observations may explain why there are contradictory reports in the literature about how FVIII reduction affects its activity. It is not clear at this stage what underpins these observed differences, but the FVIII products are all subtly different -different cellular expression systems (Chinese hamster ovary cells versus baby hamster kidney cells) and posttranslational modifications (glycosylation, sulfation, methionine oxidation, etc), different production and purification methods, and they are lyophilised and stored with slightly different excipients-one product includes reduced glutathione for instance. Another possibility is that the structural changes in FVIII caused by the reduction of Cys1899-Cys1903 decreases FVIII activity but increases its association with phospholipids or its ability such that there is a decrease in FVIII-FIX complex formation but an overall increase in tenase complex formation. A similar effect has been shown in therapeutic monoclonal antibodies where reduction of labile dilsufides increases ligand binding but decreases Fc (fragment crystallisable) function, resulting in a decrease in drug efficacy. Further work will be needed to understand why different FVIII products behave differently, but it is a noteworthy observation that reduction/oxidation affects FVIII activity as conditions such as cancer and inflammation, in which thioredoxin levels are elevated, are associated with increased risk of thrombosis. As FVIII and many other haemostasis-related plasma proteins are thioredoxin substrates, it is possible that increased Trx activity may underpin the increased thrombosis risk in these pathophysiological conditions and further supports the development of oxidoreductase inhibitors to control thrombotic risk. and CINECA under the ISCRA initiative (Cys-Surf-HP10CSOLZQ), for the availability of high-performance computing and support.