while anticoagulation therapy is widely used , it has certain undesirable side effects such as the potential to cause life - threatening hemorrhages . such bleeding complications can be mitigated , in the event of an overdose of anticoagulants , by the administration of antidotes which neutralize the anticoagulants while still avoiding thrombosis [ 1 , 2 ] . the most commonly used anticoagulants are heparin - derived drugs , which include unfractionated heparin ( ufh ) , low molecular weight heparins ( lmwhs ) , and the synthetic pentasaccharide derivatives fondaparinux and idraparinux [ 47 ] . because of its predictable dose response , almost complete bioavailability [ 4 , 7 ] , increased half - life , and no occurrence of heparin - induced thrombocytopenia , fondaparinux is becoming increasingly important in clinical medicine ; however , its widespread use is limited by a lack of a specific antidote . administration of protamine , the antidote for ufh and lmwhs , does not reverse the anticoagulant effect of fondaparinux , and hemodialysis only reduces fondaparinux plasma levels by 20% . hence , the development of a clinically safe antidote for this anticoagulant has become critical . currently , only limited experimental work has been reported for the development of an antidote to fondaparinux . it has been shown that heparinase i and the recombinant factor vii ( rviia ) can partially reverse fondaparinux in vitro ; however , these studies were limited in scope : there is no clinical data for heparinase i , and there is only one volunteer study and one clinical case for rviia . have experimentally developed antithrombin ( at ) variants as potential antidotes for heparin derivatives , including fondaparinux [ 9 , 10 ] . although the first of these was shown to neutralize fondaparinux in vitro and in vivo , its production was severely limited . to overcome this problem , a new chemically modified at variant has been produced but so far it lacks critical clinical data such as pharmacokinetics , safety , and immunogenicity . the experimental design of a polymeric antidote for fondaparinux is a daunting challenge due to the multitude of structures that need to be synthesized to arrive at the molecule with optimum binding . given the large number of possible structural configurations for a polymer antidote , a traditional trial - and - error approach to the development of a novel antidote would be a very expensive , labour - intensive , and time - consuming process . moreover , obtaining experimental data on the interactions between these polymer structures and fondaparinux is also very difficult . computer simulations therefore provide the only feasible method to screen putative polymer structures that show promise to be effective antidotes for fondaparinux , even though a rigorous experimental validation of these in silico predictions is not viable due to the arduousness of polymer synthesis and characterization . while fondaparinux has been studied in silico to some extent , the lack of experimental data for the polymer presents a challenge to the application of computational modeling techniques to antidote polymer design . computer programs for structure - based design strategies require the use of 3d structures , which are typically generated by x - ray crystallography . however , producing an x - ray crystallographic structure of the polymer candidates is difficult because they do not crystallize under normal conditions . the aim of this work is therefore to use molecular dynamic ( md ) simulations to gain a deeper insight  at a microscopic level  into the interactions between fondaparinux and individual polymer 's cationic binding groups . this information will guide the selection of favourable binding groups that will promote improved binding between the polymer antidote and fondaparinux . furthermore , the knowledge gained from these md simulations will allow for the improvement of the electrostatic model that we have previously reported to characterize the polymer - fondaparinux complex formation but which , due to the lack of interaction data , contained binding simplifications that consequently overpredicted ka values . in the next sections , the md simulations and the calculation of the free energy as well as the main equations of the electrostatic model is explained . then , the selection of the most promising binding groups based on the results obtained from md simulations and free - energy calculations is discussed . this is followed by a description of the modifications made to the previously published electrostatic model   and a discussion of the impact of these changes on the model 's predictions . all md simulations in this work were performed using the commercially available software package materials studio 5.5 ( ms ) . each md simulation system consisted of one individual fondaparinux molecule interacting with one individual cationic binding group surrounded by water molecules . for each md simulation , the system under study was first prepared , and then its free energies were calculated . the preparation of each model system was performed on an intel i5 2400 quad - core , 3.1  ghz computer and took approximately 44  h. the md calculations of the free energies of the prepared systems were then run using 48 2.66  ghz processors from the bugaboo cluster maintained by westgrid and compute / calcul canada ; these calculations took , on average , 90  h to complete for each system studied .   because electrostatic interactions drive the binding between fondaparinux and the cationic binding groups , five different binding groups ( figure 1 ) were chosen in order to determine the effect of valency on complex formation . based on our previous calculations , the r4 - 1 binding group structure was chosen as the basis for modification . therefore , the range of theoretically predicted cationic charges was chosen to have values lower and higher than that calculated for r4 - 1 ( i.e. , + 3 ) . specifically , binding groups with one , three , four , and six nitrogen , n , atoms connected by ch2ch2 linkages were selected . to observe the structural impact on the binding to fondaparinux , an additional binding group was constructed based on the structure of the binding group r4 - 1 but consisting of ch2ch2ch2 linkages between the four n atoms . the na atoms were deleted from the original 3d atomistic structure in the visualizer module of ms in order to assign a net charge of 10 to fondaparinux as determined at physiological conditions with the chemaxon pka calculator plugin ( marvin 5.5.5 , 2011 ) . the partial charges and force field types were assigned with the force field compass using the discover module . compass is an ab initio force field designed for use with a broad range of organic and inorganic molecules and polymers   and is optimized for the simulation of condensed phases . as with all current force fields , compass does not properly describe heparin [ 18 , 19 ] , in particular some of the sulfonamide functional groups found within heparin derivatives such as fondaparinux . therefore , the force field types and charges assigned by compass to the three sulfonamide oxy anions were modified to match those of the undissociated analog and the sulfonylmethoxy oxy anion , respectively . the charges were then adjusted based on a desired net charge of 10 for fondaparinux . lastly , the structure of fondaparinux was minimized after 575 iterations in the forcite module of ms using the force field compass . similarly , the 3d atomistic structures of the different binding groups were sketched using marvinsketch from chemaxon ( marvin 5.5.5 , 2011 ) ; their protonation state was calculated with chemaxon pka calculator plugin ( marvin 5.5.5 , 2011 ) , and their force field types and partial charges were assigned with the discover module using the force field compass . each model system was constructed using the amorphous cell module of ms and consisted of one deprotonated fondaparinux molecule and one protonated binding group randomly dispersed in 2,600 water molecules . a cubic simulation box was constructed with periodic boundary conditions in all directions to avoid surface effects [ 21 , 22 ] , a density of 1  g  cm since the system consists mainly of water molecules , and a side length of approximately 43   . a distance object was created between the center of mass ( centroid ) of fondaparinux and the centroid of the binding group . then , a harmonic restraint with a harmonic force constant of 100  kcal  mol   and a harmonic minimum of 21.65   was applied to this distance . energy minimization and md simulations were performed to equilibrate the system using the force field compass in the forcite module of ms . the md simulations were carried out under nvt conditions , with temperature held at 298  k by the nose - hoover thermostat [ 16 , 2224 ] . the time step was 1  fs , and the simulation time was 200  ps . this simulation time proved to be sufficient to obtain equilibrium conditions , namely , the potential energy and temperature of the system as shown in figure 2 for the representative system of fondaparinux and an r4 - 1 molecule . once the model system was relaxed , the restraint on the centroid - centroid distance was removed from the system .   using restraint forces , the intermolecular separation of fondaparinux and a binding group was sampled at equal intervals of 0.5     from an initial separation of 21.5   to a final separation of 1.5   . at each interval , a commercially available code distributed by accelrys was used to solve the following equation : 
 ( 1)fr(r)=rorkr(r(r)ro)r , rdr2kbtln(rro)+fr(ro ) , 
 where k is the harmonic force constant that enhances the restraint of the system to the r - coordinate value ro . the reaction coordinate r for each system was defined as the distance between the centroid of the fondaparinux molecule and the centroid of the binding group . a harmonic restraint constant of 100  kcal  mol   was applied over this distance . the simulation time at each interval was 60  ps , which consisted of 10  ps of equilibration followed by 50  ps of trajectory generation   for a total simulation time of 2.5  ns . the equilibration time was found to be sufficient for properties to equilibrate at each interval . the same nvt conditions , temperature , and time step previously described were used for the free - energy calculations , and each simulation was repeated 5 , 4 , 7 , 8 , and 9 times for r1 - 1 , r3 - 1 , r4 - 1 , r4 - 2 , and r6 - 1 , respectively . we have previously described an electrostatic model that provided an indication of the binding affinity between fondaparinux and a polymer structure with r4 - 1 groups . since the association rate constant , ka , is determined by diffusion and can be increased by favourable electrostatic forces [ 2831 ] whereas the dissociation rate constant , kd , is determined by short - range interactions between the molecules and is independent of long - range electrostatic forces [ 28 , 32 ] , the overall association constant , ka , and thus the affinity of a complex can be increased by optimizing the electrostatic interactions between the molecules . our model predicts ka based on these electrostatic interactions and on the following equations , first derived by schreiber and coworkers [ 28 , 33 , 34 ] : 
 ( 2)lnka = lnkaoukbt(11+ka ) , 
 where ka and ka are the association rate constants in the presence and absence of long - range electrostatic forces , respectively ; u is the electrostatic energy of interaction ; kb is the boltzmann constant ; t is the temperature of the solution ; a is the minimal distance of approach between the molecules ; and k is the debye - hckel parameter . the debye - hckel parameter is defined as 
 ( 3)k=2f2iorrt , 
 where f is the faraday constant , i is the ionic strength of the solution , o is the vacuum permittivity , r is the dielectric constant of the solution , and r is the gas constant . the electrostatic energy of interaction is defined as [ 28 , 32 ] 
 ( 4)u = ucomplexumolecule  aumolecule  b , 
 where u , the debye - hckel energy of a molecule , can be calculated from 
 ( 5)u=12i , jqiqj4orrijek(rija)1+ka . 
 in this equation , qi and qj are the charges of the atoms in the molecules , and r is the distance between the charges . as we have shown previously , we extended the empirically proven model to determine ka for the interactions of fondaparinux and a polymer structure with r4 - 1 groups . this model considers the polymer and fondaparinux , due to their structures , as a sphere and a rod , respectively . based on the smoluchowski limit for the diffusion - controlled association of two uniformly reactive molecules with these geometries , ka can be calculated as 
 ( 6)kao=4na(da+db)rx . 
 here , na is the avogadro constant , da and db are the diffusion constants of molecules a and b , respectively , and rx is the interaction radius . rx is defined as 
 ( 7)rx = lln(2l / w ) , 
 where l and w are the major and minor semiaxes of the ellipsoid . the diffusion constants are calculated as 
 ( 8)da = kbt6ra ,  db = kbt6rb , 
 where ra and rb are the hydrodynamic radii and  is the viscosity of the solvent . for modeling purposes , the interaction between the polymer and fondaparinux was assumed to be entirely electrostatic in nature . also , to simulate the functionality of the polymer , cationic binding groups on the polymer were modeled as randomly distributed on the surface of a sphere with a given radius , rh , and a minimum distance between the binding groups , rmin , was introduced . the purpose of this rmin was to account for electrostatic repulsions between the binding groups and to avoid placing these binding groups at the same location , which is physically impossible . in order to represent the solution conditions used in experiments , the following parameters were defined : i = 150  mm , t = 25c ,  = 0.90  10  kgs  m , and r = 80 . to capture the average properties of a large group of individual molecules , each simulation for a given condition consisted of constructing 1,000 unique polymers by randomly attaching the desired number of binding groups over the surface of the polymer . the model equations were then solved for each of the 1,000 polymers , and the average association rate constant was given by the geometric mean of these 1,000 runs . in the case of no polymer - fondaparinux binding , ka was set to a value of 1 . md simulations were performed to calculate the potential of mean force ( pmf)the free - energy profile along the reaction coordinate   that yields the difference in free energy between the two states of interest . these free - energy differences ( ~g ) between the unbound and bound state of the cationic groups are directly related to binding constants .   in order to investigate the effect of the binding group 's charge on fondaparinux binding , md simulations were run to follow the interaction of fondaparinux and each of the cationic binding groups ( figure 1 ) in a solution of na ions . as shown in their respective free - energy profiles , both of the cationic binding groups r1 - 1 and r3 - 1 did not have distinctive energy minima ( figures 3(a ) and 3(b ) ) . the lack of free - energy wells indicates weak binding of both r1 - 1 and r3 - 1 to fondaparinux . the poor binding of r1 - 1 was to be expected since it has been shown that the electrostatic interactions of a single protonated amine with a polyanionic molecule are weak and have to compete with salt binding under physiological conditions . r3 - 1 has an increased charge compared to r1 - 1 and would therefore be assumed to display improved binding ; however , geometry and chemical structure also play a role in complex formation . therefore , the md results suggest that r3 - 1 shows some binding to fondaparinux although it is not sufficient to overcome unfavourable orientations , thus the large degree of variability seen in the pmf for r3 - 1 compared to r1 - 1 . considering the aforementioned weak electrostatic interactions between single protonated amines and polyanionic molecules , it is not surprising that the pmf of both r4 - 1 and r6 - 1 displayed deeper free - energy wells than the lesser charged r1 - 1 and r3 - 1 ( figures 3(c ) and 3(d ) ) . in addition , the energy wells of both r4 - 1 and r6 - 1 were wide , spanning for 10   and 7   , respectively , and plateauing at a centroid - centroid distance of approximately 17.5   in both cases . this indicates that interactions between r4 - 1 and r6 - 1 and fondaparinux were much stronger and therefore felt over a larger distance than those seen in r3 - 1 . moreover , the variabilities in the pmfs of r4 - 1 and r6 - 1 were reduced compared to that of r3 - 1 which indicates that the electrostatic interactions between fondaparinux and the higher charged binding groups were strong enough to overcome unfavourable orientations . since binding affinity can be improved by increasing the cationic charges within the binding group , r6 - 1 could be expected to show a more favourable energy of interaction with fondaparinux than r4 - 1 . however , both of the binding groups had their local minima at a centroid - centroid distance of 8.5   , and both pmf profiles yielded comparable calculated g values of 2.394  kcal  mol for r4 - 1 and 2.768  kcal  mol for r6 - 1 . these results suggest that there is not a significant difference between the binding of individual r4 - 1 and r6 - 1 binding groups to fondaparinux . however , because of its additional cationic charge it could be hypothesized that when many r6 - 1 binding groups are working in concert on the surface of the polymer , the small improvement in binding affinity they show compared to r4 - 1 will be amplified and will provide stronger electrostatic interactions with fondaparinux . the complexes formed by each of these binding groups and fondaparinux can be observed in figure 4 . an alternate method of improving binding affinity between binding groups and fondaparinux is to change the spacing of the cationic charges within the binding groups . shortening the linkage between the n atoms in r4 - 1 from ch2ch2 to ch2 resulted in the protonation of only one of the amines at ph 7.4 ( as calculated with chemaxon pka calculator plugin ( marvin 5.5.5 , 2011 ) ) , and therefore , this molecule was deemed unsuitable for further study . conversely , the theoretical charge of r4 - 1 was maintained at + 3 when the linkage between the n atoms was lengthened from ch2ch2 ( 3.84   ) to ch2ch2ch2 ( 4.97   ) to form the r4 - 2 binding group ( figure 1 ) . the free - energy profile for r4 - 2 ( figure 5 ) displayed a shallow and poorly defined well . comparing the free - energy profiles of r4 - 1 ( figure 3(c ) ) and r4 - 2 suggests that increasing the spacing of cationic charges for this binding group will not generate an improvement in the binding affinity to fondaparinux . therefore , the binding groups r4 - 1 and r6 - 1 are considered the most promising binding groups for the effective neutralization of fondaparinux and will be investigated using the electrostatic model . the electrostatic model first developed in our previous work   is herein refined based on the microscopic information gathered from the md simulations described above . in the simplified model , the minimal distance of approach , a , was set to 6   since this value had been shown to give the best fit for experimental data collected for a variety of protein systems . however , the free - energy profiles of the interaction of binding groups r4 - 1 and r6 - 1 to fondaparinux ( figures 3(c ) and 3(d ) ) showed that the centroid - centroid distance between a binding group and fondaparinux was 8.5   at binding . since the refined model focused on these binding groups , the value for a the binding criterion used in the simplified model was based on the direct contact of fondaparinux with each of the three closely spaced binding groups that formed a binding site . however , the md simulation results suggest that during the binding of fondaparinux to either r4 - 1 or r6 - 1 the ionic sites of fondaparinux remain at a distance from the protonated n atoms of the binding groups . in addition , it was found that the molecules in a bound complex show significant and constant relative motion to one another . a binding site was redefined to be three or more binding groups that are within an area that allows them to come within a prescribed distance ( 8.5   ) of a fondaparinux molecule , centered around one of the binding groups . upon binding , a single charge was then assigned to this binding site ( the sum of the charges of fondaparinux and the associated binding groups ) at the location of the first binding group that formed the binding site . in order to account for bridging effects , excluding associated structural modifications that might occur , the model was altered to consider system neutralization rather than fondaparinux neutralization by allowing fondaparinux molecules to bind without the requirement for almost complete neutralization of a bound fondaparinux molecule . the fondaparinux charges that are not involved in binding to the generated polymer are therefore available for binding to another polymer in solution . the minimum distance between binding groups , rmin , for the refined model was determined as described in our previous publication . the new value of rmin was found to be 9.16   , which is larger than the radius of r4 - 1 ( 6.61   ) and r6 - 1 ( 7.10   ) . since the purpose of introducing an rmin was to account for electrostatic repulsions and steric interactions between binding groups , this new rmin is a more realistic constraint than the previously reported value , which was smaller than the radius of r4 - 1 . the metric for binding affinity used in this work is ka since , as previously described , an increase in ka would increase ka . with the refined model , ka and the number of fondaparinux bound per polymer increased with the number of attached r4 - 1 binding groups in similar fashions to what was observed using the simplified model ( figure 6 ) . however , the range for ka was reduced by 6 orders of magnitude with the refined model compared to the results obtained with the simplified form . these results indicate that although the simplified form of the model predicts overall trends in ka and thus represents an effective but rough tool for antidote discovery , the refined model with the aid of md simulation results yields more quantitatively accurate information . the effect of polymer size and number of r4 - 1 binding groups on the number of fondaparinux molecules bound per polymer and ka was investigated with the improved binding model ( figures 7 and 8) . as with the simplified model , it was observed that as the surface area of the polymer core increases , the probability of finding binding groups sufficiently close to each other to form a binding site decreases . the improvements to the model have a greater impact for smaller polymers and as shown in figure 7 , the refined model does not predict a linear increase in the number of fondaparinux molecules bound per polymer with a radius of 2  nm as was the case for the simplified model . instead , the number of fondaparinux molecules bound per polymer starts to plateau at a high number of binding groups . these results are to be expected as the number of fondaparinux bound to a polymer of a given surface area will eventually reach a maximum due to space limitations .   although the effective charge for r4 - 1 was known a priori , it was not known for the r6 - 1 binding groups . therefore , the impact of the effective charge of attached r6 - 1 binding groups on ka and the number of molecules of fondaparinux bound per polymer was also investigated in silico . for the purposes of this investigation , polymers with 20 r6 - 1 binding groups on their surface were used because charge effects are more easily observed at a high number of binding groups . as seen in figure 9 , the number of fondaparinux molecules bound per polymer was not dependent on the effective charge of the binding groups . this was to be expected as the charge of the binding group does not affect the model 's placement of binding groups on the surface of the polymer nor the binding criterion and , thus , will have no impact on an individual fondaparinux molecule finding an appropriate binding site . conversely , it was expected that ka would increase as the effective charge of the binding groups was increased ( figure 10 ) since ka is dependent on electrostatic forces and can be increased by increasing the magnitude of these forces [ 28 , 31 , 32 ] . the results show that the increase in ka with effective charge was most pronounced at an rh of 2.0  nm and that for polymers with larger radii the effective charge did not have a large impact on ka . this can be explained by the fact that electrostatic energy is inversely proportional to the distance between two charges . as the binding groups become spaced farther apart with increasing polymer radius , the impact of any increase in effective charge becomes dampened . the improved model can also be used to determine the effect of replacing r4 - 1 binding groups with r6 - 1 . as mentioned above , the effective charge of r6 - 1 is unknown ; therefore , model predictions for ka were made over a range of r6 - 1 effective charges ( figures 11(a ) and 11(b ) ) . when r6 - 1 has the same effective charge as r4 - 1 ( figure 11(a ) ) , the model predicts no difference in the calculated values of ka for r4 - 1 and r6 - 1 . this indicates that the slightly larger radius of r6 - 1 does not impact the formation of binding sites to a degree that an effect in the affinity of the polymer to fondaparinux is observed . however , it is likely that the effective charge of the r6 - 1 binding groups would be higher than that of r4 - 1 since the slightly larger r6 - 1 binding groups would most likely experience less charge shielding associated with the polymer 's core . the effect of having r6 - 1 binding groups with higher effective charge is an increase in ka values , especially for polymers with an rh of 2.0  nm and 4.0  nm ( figure 11(b ) ) . this result suggests that the binding of the polymer to fondaparinux can be improved by using r6 - 1 binding groups on small polymer cores . this is consistent with the results of the simplified model , which recommended the synthesis of polymers with low hydrodynamic size in order to increase the charge density of the polymer and , therefore , enhance ka . md simulations have enabled the refinement of a model that has previously been shown to qualitatively but not quantitatively characterize the binding of fondaparinux to a polymer - based antidote molecule . in particular , the free - energy differences calculated from md simulations have been used to determine those potential binding groups that will not bind to fondaparinux with high affinity and those that will provide strong binding . the md simulation results suggest that both r4 - 1 and r6 - 1 binding groups show high affinity for fondaparinux . the refined electrostatic model was also extended to polymers with r6 - 1 binding groups . it was found that if , as would be expected , the r6 - 1 binding groups have a higher effective charge than the r4 - 1 binding groups , improved binding affinity between the polymer and fondaparinux can be achieved . the results of this work therefore indicate that increasing the charge density of the polymer with favourable binding groups will improve complex formation between the polymer and fondaparinux . therefore , it is recommended to synthesize small polymers containing as many r6 - 1 binding groups as possible . the strengths of the model , both in its simplified and refined forms , are twofold : ( a ) minimal experimental data is required  data that is difficult to obtain for these polymers  and ( b ) the polymer structural adjustments can be determined in silico thus avoiding the costly and time - intensive synthesis of new polymers . in fact , because the synthesis and testing of even one of these polymers is such an arduous task , the model represents the only viable method for comprehensive antidote candidate screening ; even without unfeasible , rigorous experimental validation . the research presented in this work is thus a major contribution to the process of finding an antidote for fondaparinux and , therefore , will greatly impact the therapeutic field .
molecular dynamics ( md ) simulations results are herein incorporated into an electrostatic model used to determine the structure of an effective polymer - based antidote to the anticoagulant fondaparinux . in silico data for the polymer or its cationic binding groups 
 has not , up to now , been available , and experimental data on the structure of the polymer - fondaparinux complex is extremely limited . 
 consequently , the task of optimizing the polymer structure is a daunting challenge . 
 md simulations provided a means to gain microscopic information on the interactions of the binding groups and fondaparinux that would have otherwise been inaccessible . 
 this was used to refine the electrostatic model and improve the quantitative model predictions of binding affinity . 
 once refined , the model provided guidelines to improve electrostatic forces between candidate polymers and fondaparinux in order to increase association rate constants .