Applied Science, Faculty of
DSpace
UBCV
Cajiao, Adriana
Doctor of Philosophy - PhD
University of British Columbia
Anticoagulants are an essential clinical tool in modern medicine where they are administered to patients who are undergoing procedures such as cardiac surgery and kidney dialysis. They are also used to treat serious illnesses including venous thromboembolism, unstable angina, and acute myocardial infarction. However, anticoagulants have undesirable side effects, most notably excessive bleeding. Thus, antidotes are required to rapidly reverse the anticoagulant effect and treat life-threatening bleeding complications that might arise from an overdose of anticoagulants or the urgent need to perform an invasive procedure. For instance, while the synthetic anticoagulant fondaparinux has become increasingly important in clinical medicine because of its advantages over other heparin-based drugs, its use is limited by a lack of an antidote. The development of an effective, clinically safe antidote for fondaparinux is therefore critical for its widespread adoption in clinical medicine. The synthetic polymer HBSPCM (Heparin Binding Synthetic Polyvalent Cationic Macromolecule) has been found to be a promising candidate antidote for fondaparinux. However, an optimal design of HBSPCM that binds to fondaparinux with high affinity and neutralizes it has not yet been determined. A robust model is herein developed that describes the electrostatic interactions between fondaparinux and candidate HBSPCMs using data gathered from experimental work as well as molecular dynamics (MD) simulations. This model is then used to characterize the binding affinities of various putative HBSPCM structures for fondaparinux. This enables the prediction of an improved structure for the initial design of an antidote molecule. The synthesis and testing of every possible candidate HBSPCM structure would be costly and time-consuming. By quickly and efficiently predicting molecular structures that show promising binding affinities to fondaparinux, the mathematical model described in this thesis is a vast improvement over a traditional trial-and-error approach to drug design. As such, it represents an essential theoretical tool in the development of fondaparinux as an effective and safe anticoagulation treatment.
Predicting Binding Affinities of a Novel Polymer for the Neutralization of Fondaparinux by Adriana Cajiao B.A.Sc., The University of British Columbia, 2004 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Biomedical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December 2012 © Adriana Cajiao 2012 Abstract Anticoagulants are an essential clinical tool in modern medicine where they are administered to patients who are undergoing procedures such as cardiac surgery and kidney dialysis. They are also used to treat serious illnesses in- cluding venous thromboembolism, unstable angina, and acute myocardial in- farction. However, anticoagulants have undesirable side effects, most notably excessive bleeding. Thus, antidotes are required to rapidly reverse the anti- coagulant effect and treat life-threatening bleeding complications that might arise from an overdose of anticoagulants or the urgent need to perform an invasive procedure. For instance, while the synthetic anticoagulant fonda- parinux has become increasingly important in clinical medicine because of its advantages over other heparin-based drugs, its use is limited by a lack of an antidote. The development of an effective, clinically safe antidote for fonda- parinux is therefore critical for its widespread adoption in clinical medicine. The synthetic polymer HBSPCM (Heparin Binding Synthetic Polyvalent Cationic Macromolecule) has been found to be a promising candidate antidote for fondaparinux. However, an optimal design of HBSPCM that binds to fon- daparinux with high affinity and neutralizes it has not yet been determined. A robust model is herein developed that describes the electrostatic interactions between fondaparinux and candidate HBSPCMs using data gathered from ex- perimental work as well as molecular dynamics (MD) simulations. This model is then used to characterize the binding affinities of various putative HBSPCM structures for fondaparinux. This enables the prediction of an improved struc- ture for the initial design of an antidote molecule. ii Abstract The synthesis and testing of every possible candidate HBSPCM structure would be costly and time-consuming. By quickly and efficiently predicting molecular structures that show promising binding affinities to fondaparinux, the mathematical model described in this thesis is a vast improvement over a traditional trial-and-error approach to drug design. As such, it represents an essential theoretical tool in the development of fondaparinux as an effective and safe anticoagulation treatment. iii Preface A version of Chapter 3 of this dissertation has been published as: A Cajiao, B Gopaluni, E Kwok, JN Kizhakkedathu. 2011. In silico design of polymeric antidote for anticoagulant fondaparinux. Journal of Medical and Biological Engineering 31:129–134. As first author, I was involved in formulating the research design together with the principal investigators B. Gopaluni and E. Kwok. In the course of the research I wrote the necessary program code (see Appendix A) and performed all model simulations. In addition, I prepared the manuscript with subsequent editing from all co-authors. Portions of Chapter 4 and Chapter 5 have been submitted for publication as: A Cajiao, E Kwok, and B Gopaluni. 2012. Use of Molecular Dynamics for the Refinement of an Electrostatic Model for the In Silico Design of a Polymer Antidote for the Anticoagulant Fondaparinux. They are reproduced in part with permission from Macromolecules, submitted for publication. Unpublished work copyright 2012 American Chemical Society, effective on acceptance of the paper for publication. As first author, I was involved in formulating the research design together with one of the principal investigators E. Kwok. I also ran all the molecular dynamics simulations described in Chapter 4. In addition, I wrote the necessary program code (see Appendix D) and performed all model simulations in Chapter 5. Finally, I prepared the manuscript with subsequent editing from all co-authors. iv Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . xii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Blood Coagulation . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Heparin-based Anticoagulants . . . . . . . . . . . . . . . . . . 4 1.3.1 Unfractionated Heparin . . . . . . . . . . . . . . . . . 5 1.3.2 Low Molecular Weight Heparin . . . . . . . . . . . . . 5 1.3.3 Synthetic Pentasaccharides . . . . . . . . . . . . . . . . 6 1.4 Antidotes for Heparin-based Anticoagulants . . . . . . . . . . 8 1.4.1 Current Antidotes for Heparin Derivatives . . . . . . . 8 1.4.2 Heparin Binding Synthetic Polyvalent Cationic Macromolecule (HBSPCM) . . . . . . . . . . . . . . . 9 1.5 Thesis Objectives . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.6 Current Techniques to Study Heparin Interactions . . . . . . . 12 1.7 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 v Table of Contents 2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Electrostatic Interactions . . . . . . . . . . . . . . . . . . . . . 15 2.2.1 Association Rate Constants . . . . . . . . . . . . . . . 15 2.3 Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . . . 18 2.3.1 Potential of Mean Force . . . . . . . . . . . . . . . . . 20 3 In Silico Design of HBSPCM - Electrostatic Interactions Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2 Computational Materials and Methods . . . . . . . . . . . . . 23 3.2.1 Electrostatic Model Extension to Optimize HBSPCM Structure . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2.2 Parameter Selections and Binding Criterion for Compu- tation . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . 26 3.3.1 Characterization of the Binding Between HBSPCM and Fondaparinux . . . . . . . . . . . . . . . . . . . . . . . 26 3.3.2 Effect of Number of Binding Units on ka and Number of Fondaparinux Molecules Bound per HBSPCM . . . . . 29 3.3.3 Effect of HBSPCM Size on ka and Number of Fonda- parinux Molecules Bound per HBSPCM . . . . . . . . 31 3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4 Molecular Dynamics - Interaction between Fondaparinux and Candidate Polymer Binding Units . . . . . . . . . . . . 36 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.2 Computational Materials and Methods . . . . . . . . . . . . . 37 4.2.1 Binding Unit and Fondaparinux 3D System Preparation 39 4.2.2 Free-Energy Calculation . . . . . . . . . . . . . . . . . 41 4.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . 42 4.3.1 Impact of Electrostatic Interactions on Fondaparinux Binding . . . . . . . . . . . . . . . . . . . . . . . . . . 42 vi Table of Contents 4.3.2 Effect of Cationic Charge per Binding Unit on Free Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.3.3 Effect of Binding Unit Structure on Free Energy . . . . 47 4.3.4 Structural Observations . . . . . . . . . . . . . . . . . 49 4.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5 Improved Electrostatic Interactions Model for HBSPCM and Fondaparinux . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.2 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . 54 5.2.1 Experimental Charge Determination of Binding Units . 54 5.2.2 Refinement of Electrostatic Model to Optimize HBSPCM Structure . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . 56 5.3.1 Effective Charge of Binding Unit R4-1 from ζ-Potential Measurements . . . . . . . . . . . . . . . . . . . . . . . 56 5.3.2 Model Predictions for HBSPCM with R4-1 as Binding Unit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.3.3 Model Extension for HBSPCM with R6-1 as Binding Unit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . A MATLAB Code for the Electrostatic Interactions Model . 86 B Modification of Force Field Types and Charges for Fonda- parinux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 C ∆A and ∆G for Systems in the Condensed Phase . . . . . . 117 vii Table of Contents D MATLAB Code for the Improved Electrostatic Interactions Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 viii List of Tables Table B.1 Fondaparinux charge information . . . . . . . . . . . . . 116 ix List of Figures Figure 1.1 Coagulation cascade. . . . . . . . . . . . . . . . . . . . 4 Figure 1.2 Structure of the antithrombin-binding pentasaccharide sequence in heparin, fondaparinux, and idraparinux. . . 7 Figure 1.3 Schematic representation of the HBSPCM structure. . 10 Figure 1.4 Schematic representation of a cationic unit on HBSPCM. 10 Figure 3.1 Computer simulated number of fondaparinux molecules bound to HBSPCM with rigid binding units. . . . . . . 27 Figure 3.2 Computer simulated number of fondaparinux molecules bound to HBSPCM with flexible binding units. . . . . 28 Figure 3.3 Computer simulated ka and number of fondaparinux mole- cules bound to HBSPCM with flexible binding units. . 30 Figure 3.4 Computer simulated number of fondaparinux molecules bound to different sized HBSPCMs with flexible binding units. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Figure 3.5 Computer simulated ka for different sized HBSPCMs with flexible binding units. . . . . . . . . . . . . . . . . 33 Figure 4.1 Structures of the various amines used as binding units. 38 Figure 4.2 Time evolution of the potential energy and temperature for the R4-1 and fondaparinux system. . . . . . . . . . 41 Figure 4.3 Energies of the R4-1 and fondaparinux system at equi- librium. . . . . . . . . . . . . . . . . . . . . . . . . . . 43 Figure 4.4 PMF for the interaction of fondaparinux with different binding units. . . . . . . . . . . . . . . . . . . . . . . . 45 Figure 4.5 Cationic charges for R4-3 at physiological pH . . . . . 47 x List of Figures Figure 4.6 PMF for the interaction of fondaparinux with binding units of equal charge. . . . . . . . . . . . . . . . . . . . 48 Figure 4.7 MD trajectory of the R6-1 and fondaparinux system. . 50 Figure 5.1 ζ-Potential as a function of HBSPCM-2:fondaparinux molar concentration ratio. . . . . . . . . . . . . . . . . 57 Figure 5.2 Computer simulated number of fondaparinux molecules bound to HBSPCM with R4-1 binding units. . . . . . . 58 Figure 5.3 Computer simulated ka and number of fondaparinux mole- cules bound to HBSPCM with R4-1 binding units. . . . 60 Figure 5.4 Computer simulated number of fondaparinux molecules bound to different sized HBSPCMs with R4-1 binding units. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 Figure 5.5 Computer simulated ka for different sized HBSPCMs with R4-1 binding units. . . . . . . . . . . . . . . . . . 62 Figure 5.6 Computer simulated number of fondaparinux molecules bound to different sized HBSPCMs with 20 R6-1 binding units. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 Figure 5.7 Computer simulated ka for different sized HBSPCMs with R6-1 binding units. . . . . . . . . . . . . . . . . . 65 Figure 5.8 Computer simulated ka for different sized HBSPCMs with R4-1 and HBSPCMs with R6-1binding units. . . . 67 Figure B.1 Force field types for the fondaparinux sulfonamide oxy anion. . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 Figure B.2 Partial charges of the atoms in the fondaparinux sulfon- amide groups. . . . . . . . . . . . . . . . . . . . . . . . 115 Figure B.3 Final partial charges of the atoms in the fondaparinux sulfonamide groups. . . . . . . . . . . . . . . . . . . . . 116 Figure C.1 PMF, pressure, and volume for the interaction of fonda- parinux with R4-1. . . . . . . . . . . . . . . . . . . . . 119 xi Acknowledgements This thesis would not have been possible without the contributions of sev- eral people. First and foremost, I would like to express my gratitude to my supervisors Dr. Ezra Kwok and Dr. Bhushan Gopaluni for their unwavering encouragement and constructive critiques. Their expertise and advice were instrumental in the successful completion of this thesis. Special thanks go to Dr. Jayachandran Kizhakkedathu, one of my super- visory committee members, for giving me the opportunity to collaborate on this exciting project. His enthusiasm for this research project was very moti- vational. I would also like to extend my thanks to his lab members Dr. Johan Janzen and Sonja Horte for their assistance with my experiments. Similarly, I would like to thank Dr. Louise Creagh, my other supervi- sory committee member, for her continued support and invaluable insights. I greatly appreciate the time she took for our many fruitful discussions. I must also acknowledge the members of the WestGrid Bugaboo facility for their technical assistance and for generously accommodating their computer resources to allow me to complete all of my molecular dynamics simulations on time. Finally, I would like to thank my family and friends for their unconditional support, understanding, and advice even when it was not clear to them what I was doing these last years. xii Chapter 1 Introduction 1.1 Thesis Overview In modern medicine, the therapeutic use of anticoagulants is an essential clin- ical practice for the treatment and prevention of thrombotic disorders [1]. For example, anticoagulation therapy is the basis for the treatment of venous thromboembolism (VTE), which includes deep vein thrombosis (DVT) and pulmonary embolism (PE) [2, 3]. It is estimated that two million people will develop DVT each year in the United States. In 600,000 of these patients DVT will progress to PE, which will be fatal in 200,000 cases [2]. Immediate anti- coagulation minimizes the risk of thrombus extension and consequently fatal PE [2]. In addition, anticoagulants can also be administered prophylactically during some medical procedures [4]. The most commonly used anticoagulants are heparin-derived drugs [4]. Heparin derivatives include unfractionated heparin (UFH), low molecular weight heparins (LMWHs), and the synthetic pentasaccharide derivatives fondaparinux and idraparinux [2, 5–7]. LMWHs are modified forms of UFH with shorter chain length, more predictable pharmacokinetics and dose-response, and im- proved bioavailability [2, 8]. Fondaparinux and idraparinux are synthetic com- pounds that have predictable dose-responses and an almost complete bioavail- ability [2, 7]. UFH and heparin derivatives are used for the treatment of serious illnesses including VTE, unstable angina, and acute myocardial infarction [9– 11]. They are also used in patients who are undergoing cardiac surgery and kidney dialysis [4]. 1 1.1. Thesis Overview While anticoagulation therapy is widely used, it has certain undesirable side effects such as the potential to cause excessive bleeding. Bleeding compli- cations can be mitigated, in the event of an overdose of anticoagulants, by the administration of antidotes which neutralize the excess anticoagulants while maintaining sufficient levels to avoid thrombosis [12, 13]. Therefore, as new anticoagulants are regularly being created, the development of effective anti- dotes that produce minimal side effects has become a subject of major interest in the therapeutic field [3, 12, 14]. For instance, the synthetic anticoagulant fondaparinux is becoming increasingly important in clinical medicine because of its advantages over UFH and LMWHs; however, it does not yet have a specific antidote [3, 5, 7, 12, 15]. Hence, the development of a clinically safe antidote for this anticoagulant has become critical [14]. A traditional trial-and-error approach to the development of a novel anti- coagulation antidote would be a very difficult, labour-intensive, expensive, and time-consuming process. On average, each drug that came onto the market in the 1990s cost over US$800 million and took 14 years to be developed [16]. This is due to an optimization process that typically requires the synthesis and testing of hundreds or thousands of new molecules. The use of computer-based techniques, however, can speed up this process [17]. The Kizhakkedathu Laboratory in the Department of Pathology and Lab- oratory Medicine at the University of British Columbia is currently developing a promising antidote for fondaparinux, Heparin Binding Synthetic Polyvalent Cationic Macromolecule (HBSPCM). The current HBSPCM designs have been shown to effectively neutralize UFH both in vitro and in vivo and LMWHs under in vitro conditions. However they fail to adequately neutralize fon- daparinux [18]. Consequently, the structural design of HBSPCM must be modified in order to bind fondaparinux with high affinity and neutralize its biological activity. The objective of this thesis is to utilize computer simula- tions to facilitate a rapid and cost-effective development of putative synthetic polymers that can electrostatically neutralize fondaparinux. The ultimate goal 2 1.2. Blood Coagulation is to predict a small number of candidate HBSPCMs that will bind strongly to fondaparinux. The drug discovery process will thus be greatly streamlined resulting in only a few molecules that, although outside the scope of this the- sis project, will require experimental validation of pharmacological properties such as biological activity, bioavailability, and toxicity. 1.2 Blood Coagulation Coagulation is the human body’s major defense mechanism against bleeding [19] and it begins almost immediately after an injury has damaged the lining of a blood vessel (endothelium) [20]. The damage of the endothelium leads to the exposure of the blood to tissue factor which in turn initiates activation of platelets and coagulation factors [8]. The initial response to stop bleeding is the formation of a platelet plug which is produced when platelets adhere to the site of injury and become activated enhancing the initial generation of thrombin, one of the most important constituents of the coagulation cascade as discussed below [20]. The second response to bleeding is the activation of the coagulation cascade [20] (Figure 1.1), which has traditionally been described as consisting of an extrinsic and intrinsic pathway [19]. The extrinsic pathway is activated by the release of tissue factor (TF) at the site of injury while the intrinsic pathway is initiated with the contact of blood with sub-endothelial connective tissues exposed as a result of tissue damage [21]. As coagulation proceeds, both path- ways converge on a common pathway, which is the activation of factor X [21]. Factor Xa (activated factor X) then activates prothrombin to generate a burst of thrombin [5]. Thrombin is the last coagulation factor in the cascade and once generated, catalyzes the formation of fibrin strands from fibrinogen [21]. These fibrin strands bind to the platelet aggregates to form a stable platelet- fibrin thrombus or blood clot [5, 8, 20]. 3 1.3. Heparin-based Anticoagulants Figure 1.1: Schematic representation of the coagulation cascade. Under healthy conditions, there is a balance between the factors that pro- mote blood coagulation and those that inhibit it, thus avoiding pathologic thrombosis (obstruction of blood flow by a thrombus), vascular inflammation, and tissue damage [22]. Antithrombin (AT) is a plasma enzyme that acts as an irreversible inhibitor of thrombin, factor Xa, and other coagulation factors [5] through the formation of irreversible complexes [20]. The inhibitory effect of AT is enhanced via a conformational change when it binds to heparan-sulfate, a heparin-like molecule found on the surface of endothelial cells [6, 22, 23]. When additional anticoagulant activity is needed to maintain blood fluidity, anticoagulant therapy is required [5]. 1.3 Heparin-based Anticoagulants Heparin derivatives are the most commonly used anticoagulant therapies for the treatment of thromboembolic illnesses such as venous thromboembolism, unstable angina, and acute myocardial infarction [2, 9, 10]. They are also 4 1.3. Heparin-based Anticoagulants used in patients who are undergoing surgical procedures such as cardiac [4, 9] and orthopedic surgeries [11, 24] and in patients requiring extracorporeal procedures such as kidney dialysis [4]. Heparin-based anticoagulants interact with proteins (factors) in the blood coagulation cascade and produce their major anticoagulant effect via an AT dependent mechanism (Section 1.2). AT plays a major role as an inhibitor of the coagulation cascade and the conformational changes associated with the AT-heparin interaction accelerate the inhibition of coagulation factors such as thrombin (IIa) and factor Xa [21, 25–27]. Heparin derivatives, which include UFH, LMWHs, and the synthetic pentasaccharides fondaparinux and idraparinux, have their own advantages and disadvantages. 1.3.1 Unfractionated Heparin Heparin is a highly sulfated glycosaminoglycan with the highest negative charge of any biological molecule [4]. Its unfractionated form, UFH, is derived from mucosal tissues of porcine intestine or bovine lung [28]. The molecular weight of UFH ranges from 3,000 Da to 50,000 Da although most commercial forms have an average molecular weight of 15,000 Da [29]. Among the heparin derivatives, UFH is the most commonly used clinical anticoagulant because it is the only anticoagulant drug that has an antidote that completely neutral- izes it [11, 12] (Section 1.4). However, UFH can cause major side effects that include hemorrhage and heparin-induced thrombocytopenia (HIT), which re- sults in the loss of platelets [4, 12]. Frequent coagulation monitoring, a process that is inconvenient for patients and costly for the health care system [6], is therefore critical during the administration of UFH. Furthermore, UFH has a half-life, the time it takes for it to lose half of its pharmacologic activity, of only 1 to 2 hours after intravenous administration [12]. 1.3.2 Low Molecular Weight Heparin LMWHs are chemically or enzymatically modified forms of UFH that con- sist of a shorter chain length and an average molecular weight of less than 5 1.3. Heparin-based Anticoagulants 8,000 Da [28]. Several different LMWHs are commercially available, ranging from 3,600 Da to 6,500 Da [30], and each offers different anti-factor Xa and anti- IIa activities based on the polymerization methods of their preparation [9, 11]. Since protein binding is dependent on chain length and charge, LMWHs bind to less plasma proteins and therefore have more predictable pharmacokinetics and dose-response as well as improved bioavailability as compared to UFH [2, 8]. This predictable dose response obviates the need for laboratory moni- toring in most patients administered with LMWHs [12]. Moreover, compared to UFH, LMWHs have a reduced occurrence of HIT [14, 24] and an increased half-life of 4 hours [31]. Based, in part, on these advantages, LMWHs have become the anticoagulant of choice for the prevention of venous thrombosis following major trauma or orthopedic surgery [9], where they are administered subcutaneously [31]. However, despite their greater predictability with regards to dose and improved bioavailability, the clinical use of LMWHs is limited be- cause their anticoagulant effect cannot be fully reversed which can result in major bleeding events [12]. 1.3.3 Synthetic Pentasaccharides Fondaparinux and idraparinux are the synthetic forms of the pentasaccharide moiety responsible for the anti-factor Xa activity in UFH and LMWHs [2, 7] (Figure 1.2). While idraparinux is still in clinical trials, fondaparinux has re- ceived FDA approval and is in clinical use for the prevention and treatment of venous thromboembolism [24]. Therefore, this thesis focuses on fondaparinux. The anticoagulation activity of fondaparinux is derived from its AT-mediated inhibition of factor Xa which, as shown in Figure 1.1, is positioned at the start of the common pathway of coagulation [31]. Fondaparinux therefore interferes with a central step of the coagulation cascade and effectively inhibits throm- bin formation [7]. It is a chemically defined pentasaccharide with a molecular weight of 1,728 Da that exclusively binds to AT with high affinity. Therefore, fondaparinux has the advantages of a predictable dose-response and an almost complete bioavailability, which means it can be subcutaneously administered 6 1.3. Heparin-based Anticoagulants without coagulation monitoring [2, 7, 22]. Also, it has a half-life four times that of LMWHs (∼ 17 hours). This increased half-life allows for more ease of handling, especially for patients who require assistance with the administra- tion of drugs [12] and, thus, makes fondaparinux an attractive anticoagulant medication [31, 32]. Another significant advantage of fondaparinux is that HIT has not been associated with its use [5]. While fondaparinux shows promise as an anticoagulant with clear advantages over other heparin-derived drugs, its widespread use is hindered by the lack of a specific antidote [14]. Until an antidote is developed there remains the risk of major bleeding events with the use of fondaparinux [3, 12]. (a) (b) (c) Figure 1.2: Structure of: (a) the antithrombin-binding pentasaccharide se- quence in heparin; (b) fondaparinux; and, (c) idraparinux. 7 1.4. Antidotes for Heparin-based Anticoagulants 1.4 Antidotes for Heparin-based Anticoagulants 1.4.1 Current Antidotes for Heparin Derivatives Although anticoagulants are widely used in a variety of clinical applications they have potentially life threatening side effects—most notably excessive bleeding. Antidotes are therefore desired to protect against bleeding com- plications that might arise from an overdose of anticoagulants [12, 13] or the urgent need to perform an invasive procedure [33]. These antidotes neutral- ize the excess anticoagulant while maintaining levels that prevent thrombosis [12, 13]. Currently, the only antidote for UFH and LMWHs is protamine, a highly basic and cationic polypeptide, rich in positively charged arginine residues [12]. The positively charged protamine reacts with the negatively charged UFH to form a stable salt complex thus neutralizing the anticoagu- lant activity of UFH [12, 34]. A dose of 1mg of protamine will completely neutralize 90 United States Pharmacopeia (USP) units of bovine origin UFH [12]. Protamine however does not have a similar effect on LMWHs as it does not neutralize all of their anti-Xa activity [35]; on average protamine neutral- izes only 40% to 50% of the anti-Xa activity of LMWHs [14]. In addition, protamine does not reverse the anticoagulant effects of synthetic pentasaccha- rides such as fondaparinux [5, 7, 12, 14, 22]. Fondaparinux does not currently have a specific antidote [3, 5, 7, 12, 15]. The administration of protamine does not reverse the anticoagulant effect of fondaparinux and hemodialysis only reduces fondaparinux plasma levels by 20% [12]. Although it has been shown that heparinase I and the recombinant factor VII (rVIIa) can partially reverse fondaparinux in vitro, 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 [12]. With its FDA approval for the treatment and prevention of thromboembolic diseases and its advantages over UFH and LMWHs, fondaparinux has the potential to be a 8 1.4. Antidotes for Heparin-based Anticoagulants critical anticoagulant therapy if an antidote can be found. The development of a clinically safe antidote for fondaparinux is therefore of major interest in the therapeutic field [14]. 1.4.2 Heparin Binding Synthetic Polyvalent Cationic Macromolecule (HBSPCM) The Kizhakkedathu Laboratory in the Department of Pathology and Labora- tory Medicine at the University of British Columbia is currently developing HBSPCM as a promising novel antidote for heparin and its derivatives. This co-polymer incorporates binding cationic units attached to a hyperbranched polyglycerol (HPG) core with a protective outer layer of short polyethylene glycol (PEG) chains (Figure 1.3). Thus far, the cationic unit within HBSPCMs has been tris[2-(dimethylamino)ethyl]amine which consists of charged amine groups separated by ethylene groups [18]. These binding units are randomly distributed on the surface of HBSPCM (Figure 1.4). Neutralization of heparin and its derivatives proceeds via electrostatic interactions between the cationic units on HBSPCM and the negatively charged groups of the anticoagulant. While the current HBSPCM designs have been shown to effectively neu- tralize UFH and LMWHs they fail to adequately neutralize fondaparinux [18]. However, Kizhakkedathu and co-workers suspect that the binding affinity of HBSPCM and fondaparinux can be manipulated by modifying the former’s charge density by altering the structure of the cationic units and/or their dis- tribution on the HBSPCM surface. If an HBSPCM structure can be found that neutralizes the anticoagulant effects of fondaparinux, it will be an important breakthrough in clinical practice that will allow for the safe, widespread use of fondaparinux and thus an improvement in current anticoagulation practices. 9 1.5. Thesis Objectives Figure 1.3: Schematic representation of the HBSPCM structure. Figure 1.4: Schematic representation of a cationic unit on HBSPCM. 1.5 Thesis Objectives Although effective for the neutralization of UFH and LMWHs, the current design of HBSPCM does not neutralize the anticoagulant fondaparinux [18]. The effectiveness of HBSPCM to neutralize fondaparinux relies on its ability to form a complex with the anticoagulant, an ability which can be improved by increasing the rate of association of the molecules. Since in vivo conditions 10 1.5. Thesis Objectives are difficult to manipulate, it is the HBSPCM itself that must be modified in order to improve this rate of association [36]. In particular, the electrostatic interactions between the interacting molecules must be modified in order to improve their binding affinity. Given the large number of possible structural configurations for HBSPCM, an experimental approach to find the optimal structure that will provide the most stable complex with fondaparinux would be expensive and time-consuming. Moreover, experimental data on the inter- actions between HBSPCM and fondaparinux is very limited due to the nature of HBSPCM, and in silico studies on these interactions are not available (Sec- tion 1.6). It is therefore the objective of this thesis to improve the current design of HBSPCM by using computer simulations to streamline the develop- ment of a structure that will electrostatically neutralize fondaparinux. In this thesis a mathematical model is developed and used to modify, in silico, the structure of HBSPCM so as to enhance binding with fondaparinux. As previously mentioned, electrostatic interactions play a major role in the HBSPCM-fondaparinux complex formation. Therefore, the structural modifi- cations of HBSPCM carried out in this work focus on the charge characteristics of HBSPCM. Specifically, the effect of overall charge density of the HBSPCM on complex formation is investigated. In addition, in order to properly charac- terize the system, MD simulations are used to gain information at the micro- scopic level on the interactions of the cationic binding units of HBSPCM with fondaparinux. The electrostatic model, in conjunction with MD simulations, is then used to provide candidate HBSPCM structures that show potential to be effective fondaparinux antidotes. The mathematical model developed in this thesis represents a framework for the rapid in silico fine-tuning of an HBSPCM structure that will improve affinity for fondaparinux and thus provide a greater capability for neutraliza- tion. Since this is achieved with minimal experimental data, the model allows for far fewer candidate HBSPCM molecules to be synthesized and evaluated experimentally. The model therefore greatly streamlines the traditional, time- 11 1.6. Current Techniques to Study Heparin Interactions intensive trial-and-error approach to drug discovery and is thus a valuable tool in the development of an antidote for fondaparinux, a critical requirement in the therapeutic field. 1.6 Current Techniques to Study Heparin Interactions The development of specific antidotes for anticoagulants requires the char- acterization of the antidote-anticoagulant interactions [1]. Many established techniques, both experimental and computational, exist to study biomolecular interactions and some of these techniques have been shown to be effective in characterizing interactions involving heparin-based anticoagulants. For exam- ple, heparin-protein interactions have been studied experimentally using such quantitative analyses as: surface plasmon resonance (SPR) [4], isothermal titration calorimetry (ITC), affinity coelectrophoresis (ACE) [4, 37], optical biosensors, and complex trapping [37]. Computational techniques have also become established tools to study biomolecules such as proteins, DNA [4], and peptides [38]. These techniques have also been used to study heparin-protein interactions. Molecular modeling in particular has been used to: predict the heparin-binding domains of four proteins [39]; determine the importance for interaction of the distance between basic amino acids in the heparin-binding site of proteins [40]; study the poten- tial interactions between synthetic peptides and the heparin pentasaccharide unit structure that binds to antithrombin [41]; and, investigate the molecular structure of fondaparinux [1]. 12 1.6. Current Techniques to Study Heparin Interactions In addition to the computational techniques mentioned above there ex- ist some analytical models of protein interactions that can be extended to heparin and its derivatives. For instance, Schreiber and co-workers have de- veloped a protein-protein interaction model that predicts rates of association by determining the electrostatic attraction between the proteins [42–44]. This model uses a relationship between the electrostatic energy of interaction (∆U) and the rate enhancement of the formation of a protein complex [43, 44]. As such, this model is well suited to characterize HBSPCM-fondaparinux com- plex formation since electrostatic interactions are critical in heparin binding [4, 39, 41] and, due to the highly cationic nature of HBSPCM, the driving force for HBSPCM-fondaparinux complex formation. This model is developed further in Section 2.2. While fondaparinux has been studied in silico to some extent [1], it is the lack of experimental data that presents a challenge to the application of com- putational modeling techniques to study HBSPCM. Computer programs for structure-based design strategies require the use of 3D structures. These are typically generated by X-ray crystallography [45]; however, producing an X- ray crystallographic structure of HBSPCM candidates is difficult because they do not crystallize under normal conditions as, compared to other polymers, they posses very low glass transition and melting temperatures and at room temperature are highly viscous liquids [18]. Similarly, generation of a 3D struc- ture from nuclear magnetic resonance (NMR) spectroscopy experiments is not possible for HBSPCMs since, unlike proteins, they do not have a well-defined and finite structure [18]. However, one promising technique to aid in the de- termination of HBSPCM structures is molecular dynamics (MD) which has previously been used to study cationic polymers [46] and thus shows potential for the analysis of fondaparinux interactions. Specifically, the commercially available software package Materials Studio® enables MD simulations of user- generated molecular structures. This information can then be used to predict binding free energies of putative cationic binding units for HBSPCM. Further description of this technique can be found in Section 2.3.1. 13 1.7. Thesis Outline 1.7 Thesis Outline This thesis concerns the development of a mathematical model that describes the electrostatic interactions between putative HBSPCMs and the anticoagu- lant fondaparinux. This model provides a framework for the in silico modi- fication of the structure of HBSPCM candidates to find the design that will provide the most stable HBSPCM-fondaparinux complex and, thus, will effec- tively neutralize fondaparinux. The model will therefore enable the stream- lining of the costly and time-intensive experimental antidote discovery process. As a first step, a model that characterizes the binding between HBSPCM structures and fondaparinux is developed in Chapter 3, “In Silico Design of HBSPCM—Electrostatic Interactions Model”. The model determines the as- sociation rate constants as a function of polymer size and number of binding units for a range of potential HBSPCM structures and uses these constants as a metric for binding affinity. In Chapter 4, “Molecular Dynamics Simulations of Fondaparinux and Can- didate Polymer Binding Units in Water”, molecular dynamic (MD) simulations are performed to gain a deeper understanding—at the microscopic level—of the interactions between individual cationic binding units and fondaparinux. The information garnered from the work discussed in this chapter allows for the refinement of the electrostatic model developed in Chapter 3. In Chapter 5, “Improved Electrostatic Interactions Model for HBSPCM and Fondaparinux”, the MD results from Chapter 4 as well as experimental re- sults are integrated into the model developed in Chapter 3. This refined model is then extended to HBSPCMs containing more positively charged cationic binding units. Specifically, the model identifies increased charge density of the HBSPCM as the critical factor in the neutralization of fondaparinux. More- over, the model provides guidelines for the development of HBSPCM structures with the potential to be effective antidotes for fondaparinux. 14 Chapter 2 Theory 2.1 Introduction Drug discovery relies on the development of high-affinity molecules [47]. It is the goal of this thesis to develop mathematical models to streamline the discovery process and therefore determine an effective antidote molecule that has high affinity to the anticoagulant fondaparinux. As mentioned in Section 1.6, current techniques to study molecular interactions involving heparin-based anticoagulants include an empirical model based on electrostatic interactions as well as MD simulations. This chapter introduces the theoretical under- pinnings of these tools. Specifically, the role of electrostatic interactions in binding affinity is presented and an empirically proven model, which predicts the increase in rates of association and affinities of protein complexes with in- creasing Coulombic complementarity, is described and extended to HBSPCM- fondaparinux complexes. Additionally, an overview of molecular dynamics (MD) is presented with emphasis on its application to the calculation of ther- modynamic properties such as binding affinity. In particular, the potential of mean force (PMF), one of the different methods to estimate free energies from structural ensembles generated with MD simulations, is described. 2.2 Electrostatic Interactions 2.2.1 Association Rate Constants The relative affinity that two or more molecules have to form a complex is defined by the overall association constant (Ka) which is a function of the association (ka) and dissociation (kd) rate constants and given by the simple 15 2.2. Electrostatic Interactions relation [48] Ka = ka kd (2.1) Consequently, the affinity of a complex can be improved by increasing the rate at which the molecules associate and/or by decreasing the rate at which the complex dissociates into its constituent molecules. It has been observed that different forces contribute to ka and kd. The dissociation rate constant is de- termined by short-range interactions between the molecules (such as hydrogen bonds, hydrophobic interactions, and van der Waals interactions) and is inde- pendent of long-range electrostatic forces [44, 49]. In contrast, the association rate constant is determined by diffusion and can be increased by favourable Coulombic electrostatic forces [42, 44, 48, 50]. Therefore, the affinity of a com- plex can be increased by optimizing the electrostatic interaction between the reactive molecules [44]. For instance, favourable electrostatic forces between molecules can increase ka from 10 9-1010 M−1 s−1 to 1011 M−1 s−1 [42]. Such an increase would then increase Ka and thus improve affinity. It has been shown that the association rate constant for protein-protein interactions depends on ionic strength following a Debye-Hückel-like approxi- mation [44, 51, 52] ln ka = ln k o a − ∆U kBT ( 1 1 + ka ) (2.2) where ka and k o a 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 so- lution; a is the minimal distance of approach between the molecules; and, k is the Debye-Hückel parameter. A value of 6Å for a was used by Schreiber and co-workers, as it gave the best fit for the experimental data collected [42– 44, 48, 53]. 16 2.2. Electrostatic Interactions The Debye-Hückel parameter is defined as [54] k = √ 2F 2I orRT (2.3) 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. Schreiber and co-workers [44, 49] developed a model wherein the electro- static energy of interaction was calculated from the difference between the electrostatic energy of the complex and the electrostatic energy of each pro- tein ∆U = Ucomplex − UmoleculeA − UmoleculeB (2.4) where U , the Debye-Hückel energy of a molecule, can be calculated from U = 1 2 ∑ i,j qiqj 4piorrij e−k(rij−a) 1 + ka (2.5) In this equation qi and qj are the charges of the atoms in the molecules and r is the distance between the charges. For simplicity, these equations will be referred to as the Schreiber model throughout this thesis. The basal rate constant, koa, can be obtained from the intercept of a plot of the salt dependence (1/(1 + ka)) versus the experimentally measured ln ka. Schreiber and co-workers derived koa experimentally [43, 53]; however, since the purpose of this thesis is to reduce the cost and time associated with a trial-and-error approach, a theoretical approach will be used instead. Based on the Smoluchowski limit for the diffusion-controlled association of two uniformly reactive molecules where one is a small, rod-like molecule and 17 2.3. Molecular Dynamics the other is a large, spherical molecule, koa can be calculated as [50] koa = 4piNA(DA +DB)Rx (2.6) 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 Rx = l ln (2l/w) (2.7) where l and w are the major and minor semi-axes of the ellipsoid of the rod-like molecule. The diffusion constants are calculated as DA = kBT 6piηr′A and DB = kBT 6piηr′B (2.8) where, r ′ A and r ′ B are the hydrodynamic radii and η is the viscosity of the solvent. 2.3 Molecular Dynamics Molecular dynamics (MD) is an important tool in the study of systems of biomedical interest [55] and has become increasingly useful in understanding the properties of biological macromolecules in terms of their structure and microscopic interactions [56, 57]. Generally, MD consists of constructing the reactive molecules in silico and placing them in a 3D domain. In a classical all-atom MD simulation, the time-dependent information (position, velocity, and acceleration) of each atom is then obtained by integrating the Newto- nian equation of motion [55, 58, 59]. The main result of MD simulations is a trajectory file that stores at a sequence of time steps the atomic positions, velocities, and accelerations of individual atoms [58, 60]. From these trajecto- ries, average values of different properties of the system can be determined [58]. These properties can then be used to calculate, for example, the changes in 18 2.3. Molecular Dynamics binding free energy of a drug candidate [58]. MD simulations, therefore, serve as a complement to experimental results allowing the prediction of molecular properties that are either difficult [57] or impossible to measure experimentally [61]. Since MD simulations provide detail at the microscopic level, including atomic position and velocities [57, 58], they act as a bridge between micro- scopic information and the macroscopic properties of a system [56]. An important aspect of any MD simulation is the force field. A force field is the combination of a mathematical formula and a set of parameters that is used to calculate the relative energy of a molecule as a function of its atomic coordinates [60, 62]. The parameters associated with the mathematical formula describe the geometric and energetic properties of inter-particle inter- actions [62] and therefore, force fields need to be calibrated to specific sets of molecules [60, 63]. Consequently, when performing MD simulations the selec- tion of a force field depends, among other factors, on the molecules for which the force field was developed and the molecular structures being investigated [59, 60, 62]. For example, the consistent-valence force field (cvff) provides parameters for amino acids, water, and various functional groups and is thus applicable to small organic crystals and gas-phase structures [60] whereas the force field COMPASS (Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies) was parameterized for most common organic molecules, small inorganic molecules, and polymers [60] and optimized for the simulation of condensed phases [64]. Once a force field has been selected and the system of molecules under study optimized, or minimized, with respect to the potential energy a dy- namics simulation can be run [60]. The classical equations of motion are the basis of this simulation [55] which are modified, when necessary, to deal with the effects of pressure and temperature of the system [60]. Integrating Newton’s equation of motion allows the study of a system in which energy is conserved [58, 60]. However, most systems are exposed to heat exchange with the environment which results in their total energy no longer being conserved 19 2.3. Molecular Dynamics and, therefore, extended forms of MD are used [60]. These extended forms of MD use several different methods to control pressure and temperature and as a result generate different statistical ensembles [60]: the collection of the various microscopic states of a system that result from different initial condi- tions (atomic position and velocities) and that together correspond to a single macroscopic, or thermodynamic, state [60, 65]. For example, the canonical or constant-temperature, constant-volume ensemble (NVT) is obtained by con- trolling the thermodynamic temperature with a fixed number of atoms, N, and a fixed volume, V [58, 60]. 2.3.1 Potential of Mean Force Free energy is a fundamental concept in thermodynamics [66, 67] and a mea- surement of the probability of finding a system in a given state [68]. Therefore, free-energy differences are directly related to important quantities such as ad- sorption constants, solubilities, and binding constants [69]. The calculation of free-energy differences is a major task in computational science with appli- cations ranging from biochemical processes to rational drug design [67, 70]. However, this task is complicated as it requires sampling of the most probable configurations of the two states under study. These states might be separated by significant free-energy barriers and may not allow all important configura- tions to be sampled [68, 69]. Different sampling methods have been proposed to determine relative free energies [66, 71]. One such method is the potential of mean force (PMF) [72] . In this method, the two states of interest are con- nected by a reaction path, often referred to as the reaction coordinate, along which the process proceeds. Free energy is then calculated at each position along the path [69, 72]. This free-energy profile ultimately yields the difference in free energy between the two end states [69] and is called a PMF as it can be obtained by integrating the mean force,−dFR(R”)/dR”, along the reaction coordinate, R [73]. The PMF of a system at a particular value, R′, of the 20 2.3. Molecular Dynamics reaction coordinate, R, is given by [73]: FR(R ′) = FR(∞) + ∫ R′ ∞ dFR(R”) dR” dR” (2.9) The computation of relative free energies as a function of R’ from standard simulations can be subject to sampling errors. To overcome these errors a restraining or umbrella potential term can be added to focus sampling on a particular range of R′ values [69]: FR(R ′) = − ∫ R′ R′o 〈Kr(R(r)−Ro)〉r,R′′dR′′ − 2kBT ln(R′/R′o) + FR(R′o) (2.10) where Kr is the harmonic force constant that enhances the restraint of the system to the R-coordinate value Ro. Since this equation is derived in the canonical (NVT) ensemble, it calculates the Helmholtz free energy (A) along R [69]. If pressure instead of volume is kept constant (NPT ensemble), the Gibbs free energy (G) is obtained. However, for systems in the condensed phase the change in volume is negligible and the free-energy changes ∆A and ∆G approximate one another [67, 70] (see Appendix C). 21 Chapter 3 In Silico Design of HBSPCM - Electrostatic Interactions Model 3.1 Introduction Although HBSPCM structures have been developed that effectively neutralize both UFH and LMWHs [74], a structure with specific neutralizing activity towards fondaparinux has yet to be discovered. Given the large number of possible HBSPCM structural configurations, an experimental approach to find the one that will provide the most stable complex with fondaparinux would be expensive and time-consuming. To streamline the drug discovery process, a mathematical model is presented in this chapter that will allow for the char- acterization of the binding between putative HBSPCM structures and fonda- parinux. Specifically, the association rate constants for potential HBSPCM structures are determined by the model and used as a metric for binding affin- ity. These results will then represent a framework for the rapid in silico fine- tuning of the structure of an HBSPCM that will display increased affinity for fondaparinux and thus a greater capability for neutralization. Due to the very limited information available on the molecular structure of HBSPCM (Section 1.6) and on the interactions between these molecules and fondaparinux, the model will provide a qualitative description of binding affinities. This qualita- tive data will be sufficient to identify the HBSPCM structure(s) that promises to give the most stable HBSPCM-fondaparinux complex. 22 3.2. Computational Materials and Methods 3.2 Computational Materials and Methods 3.2.1 Electrostatic Model Extension to Optimize HBSPCM Structure The Schreiber model (Equations 2.2, 2.4, and 2.5) has been used to calculate the rate of association of protein-protein binding based on a certain basal rate (i.e. rate of association in the absence of electrostatic forces) as well as on the Coulombic complementarity between the charged residues on the surface of the proteins [44, 49]. As with protein-protein interactions, it has been demon- strated that the cationic segments of polymers similar to HBSPCM interact electrostatically with the negative charges of bio-macromolecules and, there- fore, form stable complexes [75]. Furthermore, it has been shown that the binding between antithrombin and the heparin pentasaccharide unit structure responsible for binding to antithrombin—almost identical sequence of sugars as those found in fondaparinux [2, 6, 14]—is mediated primarily through elec- trostatic contacts [41]. Therefore, it is clear that electrostatic interactions play a major role in the formation of the HBSPCM-fondaparinux complex. It is thus the approach of this work to use the empirically proven Schreiber model to determine ka for the interactions between HBSPCM and fondaparinux. As mentioned in Section 2.2.1, kd is not affected by electrostatic interactions; therefore, changes in ka will be an accurate representation of changes to the equilibrium constant, Ka, and thus an appropriate metric for binding affinity. 3.2.2 Parameter Selections and Binding Criterion for Computation Since it is only possible to theoretically predict binding constants for basic cases [76] a simplified approach is encouraged for design purposes. Therefore, in the binding model HBSPCM and fondaparinux were considered, based on their structures, as a sphere and a rod respectively (Section 1.4.2 and Section 1.3.3). Also, the interaction between HBSPCM and fondaparinux was assumed to be entirely electrostatic in nature. This is because under physiological con- 23 3.2. Computational Materials and Methods ditions, both HBSPCM and fondaparinux are highly charged molecules [74]. Thus, as explained in Section 2.2.1, electrostatic interactions are considered the major binding force while short-range interactions between the molecules be- come insignificant in comparison. To simulate the functionality of HBSPCM, cationic binding units on HBSPCM were modeled as randomly distributed on the surface of a sphere with a given radius, rH . This is a reasonable and conve- nient representation of HBSPCM since binding units are randomly attached at the time of HBSPCM synthesis and these macromolecules exhibit a globular shape [74]. Also, because of the relatively small size of fondaparinux compared to HBSPCM, it can be assumed that the size and spherical shape of HBSPCM will not be altered by the binding of fondaparinux. The binding units were represented as having a charge of +3, as this was their charge at physiological conditions (pH 7.4) calculated using the ChemAxon pKa Calculator Plugin (Marvin 5.5.5, 2011) [77]. A minimum distance between the binding units, rmin, was introduced. The purpose of this rmin was to account for electrostatic repulsions between the binding units and to avoid placing charges at the same location, which is physically impossible. This rmin was not known a priori nor was the flexibility of the binding units. Both of these characteristics were determined by comparing simulation results to the experimentally measured number of fondaparinux molecules bound to a specific, undisclosed due to pro- prietary reasons, HBSPCM structure. The experimental value was obtained from isothermal titration calorimetry (ITC) experiments [74]. A criterion for the binding of a fondaparinux molecule to HBSPCM was required in order to predict when complex formation would occur. The charge of a fondaparinux molecule under physiological conditions was calculated to be -10 using the ChemAxon pKa Calculator Plugin (Marvin 5.5.5, 2011) [77]. Complete charge neutralization of fondaparinux would thus require four bind- ing units. However, simulations run using a criterion of four closely spaced binding units did not yield any binding of fondaparinux. This was because the binding units were never close enough to allow a fondaparinux molecule to contact each one. Therefore, the binding criterion was defined as three binding 24 3.2. Computational Materials and Methods units that were spaced sufficiently close so that a fondaparinux molecule could stretch among them (see Appendix A). This binding definition was proposed in order to come as close as possible to complete charge neutralization of a bound fondaparinux molecule. In addition, it was assumed that fondaparinux binds to only one HBSPCM, so bridging effects were neglected. In order to determine the Debye-Hückel energy of a fondaparinux molecule, its structure was studied using the 3D chemical structure viewer Jmol [78]. The Debye-Hückel energy of the complex was found by considering each bound fon- daparinux molecule and its three associated binding units as one binding site. A binding site was assigned a single charge (the sum of the charges of fon- daparinux and the associated binding groups) located at the center point of the binding site. The minimal distance of approach, a, was set to 6Å since this value has been shown to give the best fit for experimental data collected for a variety of protein systems [44]. In order to represent the solution condi- tions used in experiments the following parameters were defined: I = 150 mM, T = 25o C, η = 0.90× 10−3 kgs−1 m−1, 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 HBSPCMs by randomly attaching the desired number of binding units over the surface of HBSPCM. The model equations were then solved for each of the 1,000 HBSPCMs and the average association rate constant was given by the geometric mean of these 1,000 runs. In the case of no HBSPCM-fondaparinux binding ka was set to a value of 1. The MATLAB code for the model programs used in this chapter is found in Appendix A. 25 3.3. Results and Discussion 3.3 Results and Discussion 3.3.1 Characterization of the Binding Between HBSPCM and Fondaparinux In order to construct the model, certain characteristics of the HBSPCM- fondaparinux binding that were not known a priori needed to be determined; in particular, the nature of the binding units and how they were distributed on the surface of the HBSPCM core. Specifically, it was required to determine: a) whether the binding units should be approximated as point charges or whether they had flexibility about the point of attachment that allowed them to ex- tend the range of their electrostatic influence, and b) the appropriate value of rmin. Simulations were first performed for an HBSPCM of set size rH = 4.0 nm (based on a medium molecular weight globular polymer [74]) over a range of rmin and a range of number of binding units under the assumption that the binding units were rigid. These simulations were then repeated with flexible binding units. Figure 3.1 and Figure 3.2 show the calculated number of fon- daparinux molecules bound per HBSPCM for these two cases. These results were then compared to the experimentally derived number of fondaparinux molecules bound per HBSPCM to determine both the nature and distribution of the binding units. Under the assumption of rigid binding units, and for a specific number of attached binding units, the calculated number of fondaparinux molecules bound per HBSPCM was well below the experimental value (undisclosed due to proprietary reasons) for every rmin tested. Conversely, when the bind- ing units were assumed to be flexible, a sufficient number of fondaparinux molecules were bound at the appropriate number of binding units. This find- ing indicates that the binding units have an inherent flexibility and cannot be treated as rigid molecules. This flexibility was therefore incorporated into the mathematical model for all future simulations. Furthermore, it was deter- mined from Figure 3.2 and the undisclosed ITC experimental data [74] that 26 3.3. Results and Discussion Figure 3.1: Computer simulated number of fondaparinux molecules bound to one HBSPCM for HBSPCMs with an rH of 4.0 nm that dif- fered in the number of rigid binding units attached to the surface. HBSPCM contained 3 to 20 rigid binding units. Simulations for each HBSPCM tested were run at different rmin values. Each sim- ulation was repeated 1,000 times; the resulting arithmetic means are shown. 27 3.3. Results and Discussion Figure 3.2: Computer simulated number of fondaparinux molecules bound to one HBSPCM for HBSPCMs with an rH of 4.0 nm that differed in the number of flexible binding units attached to the surface. HBSPCM contained 3 to 20 flexible binding units. Simulations for each HBSPCM tested were run at different rmin values. Each sim- ulation was repeated 1,000 times; the resulting arithmetic means are shown. 28 3.3. Results and Discussion an rmin of 5.33 Å was required; this was also incorporated into the model going forward. These results were expected since the methylated form of the cationic binding unit used in HBSPCM, tris(2-aminoethyl)-amine, has been previously reported as being a flexible molecule [79]. Also, because of the charged nature of these binding units it is expected that there would be some spread in their distribution due to repulsive forces. 3.3.2 Effect of Number of Binding Units on ka and Number of Fondaparinux Molecules Bound per HBSPCM After characterizing the HBSPCM-fondaparinux binding, the effect of the number of binding units on ka and on the number of fondaparinux bound per HBSPCM was studied. Figure 3.3 shows that as the number of binding units on the surface of HBSPCM increases from 3 to 20, ka increases from 100 M−1 s−1 to 1020 M−1 s−1 while the average number of bound fondaparinux molecules increases from 0.1 to almost 3. This is expected as an increase in binding units increases the likelihood that the binding criterion is met. These results also indicate that an HBSPCM with a radius, rH , of 4.0 nm has poor binding efficiency. For example, the binding of HBSPCM and fondaparinux was found to be very unlikely with three binding units and almost twelve bind- ing units were required to ensure that at least one molecule of fondaparinux binds to an HBSPCM. In fact, regardless of the number of binding units fewer than half are involved in binding. Since kd is independent of long-range electrostatic forces [44, 49] and can therefore be considered constant for the purposes of this model, an increase in ka will necessarily increase the overall association constant and thus the affinity of the HBSPCM-fondaparinux complex. Although high ka values of up to 1010 M−1 s−1 have been reported for different protein pairs [43], Figure 3.3 shows ka values as high as 10 20 M−1 s−1. It should be noted that due to the simplifications assumed in the mathematical model, predicted values for 29 3.3. Results and Discussion Figure 3.3: Computer simulated ka and number of fondaparinux molecules bound per HBSPCM with rH of 4.0 nm that differed in the number of flexible binding units separated by an rmin of 5.33 Å. ka and number of fondaparinux molecules bound to an HBSPCM are the geometric and arithmetic means respectively of 1,000 calculated values. The error bars represent the 95% confidence intervals. 30 3.3. Results and Discussion ka are overestimated and, therefore, cannot be considered absolute ka values. However, the model predicts a trend in ka values with electrostatic enhance- ment and thus is not only useful to investigate the nature and distribution of binding units on HBSPCM (Section 3.3.1) but can also be used to study the effect of other HBSPCM structural parameters on binding affinity, such as HBSPCM size (Section 3.3.3). 3.3.3 Effect of HBSPCM Size on ka and Number of Fondaparinux Molecules Bound per HBSPCM With the ka and the number of bound fondaparinux molecules per HBSPCM determined for an HBSPCM of size rH = 4.0 nm, the model was used to inves- tigate the effect of HBSPCM structural changes on fondaparinux binding. Of particular interest was the investigation of the size of the HBSPCM core and number of binding units attached to its surface. As can be seen in Figure 3.4, it was found that for a given number of binding units, increasing rH decreases the number of fondaparinux molecules that will bind to an HBSPCM. As ob- served during HBSPCM synthesis, as the surface area of the HBSPCM core increases so does the spacing of the randomly placed binding units [74]. This results in a decreased probability of having binding units sufficiently close to each other to satisfy the binding criterion. Consequently, the amount of fon- daparinux bound to an HBSPCM decreases as rH increases. The size of the HBSPCM core was also seen to affect ka which, at a given number of binding units, increased with decreasing HBSPCM size (Figure 3.5). This was due to the increased energy associated with more closely spaced charged groups on the smaller polymers. As can be seen from Equation 2.5, the Debye-Hückel energy of a molecule is dependent on the inverse of the distance between its charged atoms. Decreasing this distance will therefore increase U and in turn ka. 31 3.3. Results and Discussion Figure 3.4: Computer simulated number of fondaparinux molecules bound per HBSPCM with an rH of 2.0 nm, 4.0 nm, 6.0 nm, and 10.0 nm. All different sized HBSPCMs were tested with a number of flexible binding units ranging from 3 to 20 using an rmin of 5.33 Å. Number of fondaparinux molecules bound to an HBSPCM is the arithmetic mean of 1,000 calculated values. The error bars represent the 95% confidence interval. 32 3.3. Results and Discussion Figure 3.5: Computer simulated ka for HBSPCMs with an rH of 2.0 nm, 4.0 nm, 6.0 nm, and 10.0 nm. All different sized HBSPCMs were tested with a number of flexible binding units ranging from 3 to 20 using an rmin of 5.33 Å. ka is the geometric mean of 1,000 calcu- lated values. The error bars represent the 95% confidence interval. 33 3.4. Conclusions The results from Figure 3.4 and Figure 3.5 indicate that decreasing the size of the HBSPCM increases the efficiency of fondaparinux binding. With smaller HBSPCMs a greater proportion of the binding units are involved in the binding of fondaparinux. Moreover, ka for small HBSPCMs was found to be significantly greater compared to larger HBSPCMs with equivalent numbers of bound fondaparinux molecules. These results suggest that decreasing the size of the HBSPCM core will improve HBSPCM-fondaparinux binding. Similarly to the results shown in Figure 3.3, the ka values seen in Figure 3.5 are unrealistically high especially at an rH of 2 nm. These high ka val- ues are the result of the charge groups being clustered closely together over a small surface area and indicate that the predicted rmin is artificially low and the binding criterion too constrained. In order to predict more realistic values of ka a more detailed structural analysis of the interactions between fondaparinux and the binding units is required. However, the trends predicted by the simplified model and seen in Figure 3.5 provide valuable insight into how to ameliorate fondaparinux binding in a rapid and efficient manner in the absence of computationally expensive molecular dynamic simulations. 3.4 Conclusions The empirical model developed in this chapter has been used to characterize the binding of HBSPCM to fondaparinux. By first fitting the model to experimental results the appropriate nature and distribution of the cationic binding units was determined. It was found that the binding units on HBSPCM are flexible and are spaced a minimal distance of 5.33 Å apart. Once this had been determined, the model then predicted the effect on binding affinities of certain structural parameters of HBSPCM. This provided a simple guideline for the improvement of HBSPCM-fondaparinux com- plex formation: synthesis of HBSPCMs having relatively low hydrodynamic size will increase charge density and will enhance ka and binding unit efficiency. 34 3.4. Conclusions The development in this chapter of a mathematical model provides an ex- cellent tool with which to refine putative HBSPCM structures in order to max- imize fondaparinux binding. With minimal experimental data the model has already been able to provide guidelines to improve the structure of HBSPCM. Taking into account synthesis procedures as well as physiological consider- ations, the HBSPCM size and number of binding units that maximize both binding unit efficiency and binding affinity has been determined. By producing these guidelines in silico far fewer candidate HBSPCM need to be synthesized and evaluated experimentally. The model therefore has greatly streamlined the traditional trial-and-error approach to antidote development thus decreas- ing costs and saving time. The model, in its current form, relies on simplifications and assumptions of certain important HBSPCM structural parameters. Most critically, the value of the minimal distance of approach was not known for this system nor were the binding conformations of fondaparinux and the binding units. Consequently, absolute values for the association rate constant were overestimated. In the following chapter, MD simulations are used to provide essential data for the binding of fondaparinux to different cationic binding units. This data will en- able the refinement of the model and, therefore, the quantitative improvement of the model’s predictions. 35 Chapter 4 Molecular Dynamics - Interaction between Fondaparinux and Candidate Polymer Binding Units 4.1 Introduction In Chapter 3 an electrostatic model was developed that characterized the inter- action between fondaparinux and potential HBSPCM candidates. The model predicted qualitative trends in ka as a function of HBSPCM size and num- ber of binding units. These trends provided valuable insight into the effect of HBSPCM structure on fondaparinux binding affinity and indicated that increasing HBSPCM’s charge density would ameliorate neutralization of the anticoagulant. However, due to simplifications the model quantitatively over- predicted ka values. The aim of this chapter is therefore to use MD simula- tions to gain a deeper insight—at a microscopic level—into the interactions between fondaparinux and individual cationic binding units. This information will guide the selection of favourable binding units that will promote improved binding between HBSPCM and fondaparinux. Furthermore, the knowledge gained from MD simulations will allow for the improvement of the electro- static model. As a first step, MD simulation results are used to corroborate the funda- mental assumption of the aforementioned model that electrostatic interactions play the key role in HBSPCM-fondaparinux complex formation. Subsequently, free-energy profiles (PMF) are generated for a variety of putative cationic bind- 36 4.2. Computational Materials and Methods ing units. These profiles are then used to discriminate the binding units that show the most promise to bind to fondaparinux with high affinity. In addition, MD trajectories are inspected to gain a better understanding of the structural configurations of the individual binding units and fondaparinux during their interaction. 4.2 Computational Materials and Methods Five different binding units were studied and their structures are shown in Figure 4.1. Because electrostatic interactions drive the binding between fon- daparinux and HBSPCM, the candidate binding units were chosen in order to determine the effect of valency on complex formation. Since the R4-1 binding unit has already been shown to be partially effective at neutralizing fonda- parinux, its structure was chosen as the basis for modification. Therefore, the range of theoretically predicted cationic charges was chosen to have val- ues lower and higher than that calculated for R4-1 (i.e. +3). The binding units thus differ in their number of constituent nitrogen, N , atoms. Specif- ically, binding units with one, three, four, and six N atoms connected by −CH2 − CH2− linkages were selected. To observe the structural impact on the binding to fondaparinux, an additional binding unit was constructed based on the structure of the binding unit R4-1 but consisting of −CH2−CH2−CH2 linkages between the four N atoms. Each amine in the binding units was N - methylated in order to reduce reactivity and toxicity [18]. The molecular weights of the binding units ranged from 63 g mol−1 for the single amine (R1- 1) to 316 g mol−1 for the hexa-amine (R6-1). It was desired to perform the MD simulations at the physiological pH of 7.4. In order to maintain this pH the required protonation states were set for fondaparinux and the N atoms within the binding units. These protonation states were determined prior to the simulations using the ChemAxon pKa Cal- culator Plugin (Marvin 5.5.5, 2011) [77] and can be seen in Figure 1.2(b) for fondaparinux and Figure 4.1 for the binding units. 37 4.2. Computational Materials and Methods (a) (b) (c) (d) (e) Figure 4.1: Structures of the various amines used as binding units: (a) R1-1; (b) R3-1; (c) R4-1; (d) R4-2; and, (e) R6-1. The nitrogen atoms are connected by −CH2 − CH2− linkages in (b), (c), and (e) and by −CH2−CH2−CH2 linkage in (d). The protonation state is for physiological pH of 7.4 and was calculated with ChemAxon pKa Calculator Plugin (Marvin 5.5.5, 2011) [77]. 38 4.2. Computational Materials and Methods All MD simulations in this work were performed using the commercially available software package Materials Studio® 5.5 (MS). The computational expense of MD simulations increases with system complexity [60, 65]. There- fore, due to resource and time constraints each MD simulation performed in this work consisted of one individual fondaparinux molecule interacting with one individual cationic binding unit surrounded by water molecules. For each MD simulation the system under study was first prepared and then its free en- ergies 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 simulations to calculate the free energies of the prepared systems were then run using 48 2.66 GHz processors from the Bugaboo cluster main- tained by WestGrid and Compute/Calcul Canada; these calculations took, on average, 90 h to complete for each system studied. 4.2.1 Binding Unit and Fondaparinux 3D System Preparation The initial 3D atomistic structure of fondaparinux was obtained from the DrugBank database [80]. Since the protonation state of a residue is kept con- stant in MD simulations [81] and fondaparinux was found to be deprotonated at physiological pH (see Section 3.2.2), the Na+ atoms were deleted from the original 3D atomistic structure in the Visualizer module of MS in order to assign the appropriate net charge to fondaparinux. The partial charges and force field types were assigned with the force field COMPASS using the Dis- cover module. COMPASS is an ab initio force field designed for use with a broad range of organic and inorganic molecules and polymers [60] and is op- timized for the simulation of condensed phases [64]. As with all current force fields, COMPASS does not properly describe heparin [82, 83], in particular some of the sulfonamide functional groups found within heparin derivatives such as fondaparinux [41]. Therefore, some modifications were made to the force field types and charges assigned by COMPASS to the three sulfonamide 39 4.2. Computational Materials and Methods oxy anions. First, the undissociated analog of the sulfonamide oxy anion was built with the Visualizer module and its force field types were properly as- signed using COMPASS. These force field types were then manually assigned to the three dissociated sulfonamide oxy anions in fondaparinux. Next, since COMPASS properly assigned the sulfonylmethoxy oxy charges, these charges were assigned to the three sulfonamide oxy anions as a first approximation. These charges were then adjusted based on a desired net charge of -10 for fondaparinux. Details of this procedure can be found in Appendix B. Because sketching a structure can result in a high energy configuration of a molecule, which can lead to false simulation results [60], 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 units un- der study were sketched using MarvinSketch from ChemAxon (Marvin 5.5.5, 2011) [77] and their protonation state was calculated with ChemAxon pKa Calculator Plugin (Marvin 5.5.5, 2011) [77]. The protonated structures were transferred to the Visualizer module of MS and their force field types and partial charges were assigned with the Discover module using the force field COMPASS. The structures were then minimized in the Forcite module using COMPASS. Each model system was constructed using the Amorphous Cell module of MS and consisted of one deprotonated fondaparinux molecule and one proto- nated binding unit randomly dispersed in 2,600 water molecules. Appropriate amounts of sodium counterions (Na+) were added to achieve charge neutrality. A cubic simulation box was constructed with periodic boundary conditions in all directions to avoid surface effects [65, 84], a density of 1 g cm−3 since the system consists mainly of water molecules, and a side length of approximately 43 Å. A distance object was created between the center of mass (centroid) of fondaparinux and the centroid of the binding unit. Then, a harmonic restraint with a harmonic force constant of 100 kcal mol−1 Å2 and a harmonic minimum 40 4.2. Computational Materials and Methods (a) (b) Figure 4.2: Time evolution of (a) the potential energy (kcal mol−1) and (b) the temperature (K) for a model system consisting of 1 R4-1 molecule, 1 fondaparinux molecule, 7 Na+ atoms, and 2,600 water molecules. of 21.65 Å was applied to this distance to ensure fondaparinux and the binding unit were separated from each other by at least 21.5 Å, the desired initial dis- tance for the PMF simulations. 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 thermo- stat [60, 84–86]. 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 the temperature of the system as shown in Figure 4.2 for the representative system of fondaparinux and an R4-1 molecule. Once the model system was relaxed, the restraint on the centroid-centroid dis- tance was removed from the system. 4.2.2 Free-Energy Calculation Using restraint forces, the intermolecular separation of fondaparinux and a binding unit was sampled at equal intervals of 0.5 Å [87] from an initial sep- aration of 21.5 Å to a final separation of 1.5 Å and Equation 2.10 was solved at each interval using a commercially available code distributed by accelrys®. 41 4.3. Results and Discussion 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 unit. A harmonic restraint constant of 100 kcal mol−1 Å2 was applied over this dis- tance. The simulation time at each interval was 60 ps, which consisted of 10 ps of equilibration followed by 50 ps of trajectory generation [88], for a total sim- ulation time of 2.5 ns. The equilibration time was found to be sufficient for properties to equilibrate at each interval. The same NVT conditions, temper- ature, and time step described in Section 4.2.1 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. 4.3 Results and Discussion 4.3.1 Impact of Electrostatic Interactions on Fondaparinux Binding In Chapter 3 it was argued that electrostatic interactions play the major role in the HBSPCM-fondaparinux complex formation while short-range interac- tions, such as van der Waals forces, are insignificant in comparison. The validity of this argument was examined as a first step in using MD simulations to investigate the interaction between fondaparinux and HBSPCM. Figure 4.3 shows the attractive forces between the cationic binding unit R4-1 and fon- daparinux at the representative distance of 10 Å. As seen in this figure, the energy of electrostatic interaction between the two molecules is an order of magnitude greater than the other forces. Thus, electrostatic interactions are indeed the driving force behind the binding of the anionic fondaparinux and the cationic HBSPCM. This confirms that the model developed in Chapter 3 is suitable to describe HBSPCM-fondaparinux binding and that the affin- ity of the HBSPCM-fondaparinux complex can be improved by increasing the electrostatic interactions between the molecules. 42 4.3. Results and Discussion Figure 4.3: Energies of the R4-1 and fondaparinux system in the presence of water molecules and Na+ atoms at equilibrium. Values are averages of 7 MD simulations at a representative distance of 10 Å between R4-1 and fondaparinux. Similar results were found for R1-1, R3-1, R4-2, and R6-1. 4.3.2 Effect of Cationic Charge per Binding Unit on Free Energy Since the binding of fondaparinux to HBSPCM is mediated primarily by elec- trostatic interactions, the charge and structure of the cationic binding units will have a large impact on affinity. MD simulations were therefore performed on putative binding units in an attempt to determine the structure that would bind most strongly to fondaparinux. More specifically, MD simulations were performed to calculate the PMF. The PMF is the free-energy profile along the reaction coordinate [72] that yields the difference in free energy between the two states of interest [69]. These free-energy differences (∼∆G, see Appendix C) between the unbound and bound state of the binding units are directly re- 43 4.3. Results and Discussion lated to binding constants [69]. In order to investigate the effect of the binding unit’s charge on fondaparinux binding, MD simulations were run to follow the interaction of fondaparinux and each of the binding units containing different numbers of protonated amines (Figure 4.1(a)–(c) and Figure 4.1(e)) in a solu- tion of Na+ ions. The MD simulation results for these four different cationic units (R1-1, R3-1, R4-1, and R6-1) are shown in Figure 4.4 in the form of PMF profiles as a function of centroid-centroid distance. As shown in their respective free-energy profiles both of the cationic bind- ing units R1-1 and R3-1 did not have distinctive energy minima (Figure 4.4(a) and Figure 4.4(b)). The lack of free-energy wells indicates weak binding of 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 proto- nated amine with a polyanionic molecule are weak and have to compete with salt binding under physiological conditions [89]. R3-1 has an increased charge compared to R1-1 and would therefore be assumed to display improved bind- ing as the interaction between the cationic binding units and fondaparinux is electrostatic in nature. 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 but 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 [89], it is not surpris- ing that the PMF of both R4-1 and R6-1 (molecules with theoretical charges at physiological conditions of +3 and +4 respectively) displayed deeper free- energy wells than the lesser charged R1-1 and R3-1 (Figure 4.4(c) and Fig- ure 4.4(d)). In addition, the energy wells of both R4-1 and R6-1 were wide, spanning for 10 Å and 7 Å respectively, and plateauing at a centroid-centroid distance of approximately 17.5 Å in both cases. This indicates that interac- tions between these binding units and fondaparinux were much stronger and 44 4.3. Results and Discussion (a) (b) (c) (d) Figure 4.4: Free-energy profile (PMF) for the interaction of fondaparinux with: (a) R1-1, (n= 5); (b) R3-1, (n= 4); (c) R4-1, (n= 7); and, (d) R6- 1, (n= 9). Calculated values are the Helmholtz free energy which approximates the Gibbs free energy for systems in the condensed phase as explained in Section 2.3.1 and Appendix C. Simulations were performed using a step size of 0.5 Å and a simulation time of 60 ps at each interval. The error bars represent the 95% confidence intervals for n number of replicates. 45 4.3. Results and Discussion therefore felt over a larger distance than those seen with 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 units were strong enough to overcome un- favourable orientations. A design of HBSPCM with R4-1 binding units, HBSPCM-2, has been syn- thesized and tested by the Department of Pathology and Laboratory Medicine at the University of British Columbia. The results from these studies have shown that this HBSPCM structure neutralizes the biological activity of fon- daparinux with a 50% success rate [90]. Since binding affinity may be improved by increasing the cationic charges within the binding unit, R6-1 could be ex- pected to show a more favourable energy of interaction. Comparing the PMF of R4-1 to that of R6-1 though showed little difference. Both of the binding units had their local minima at a centroid-centroid distance of 8.5 Å and both PMF profiles yielded comparable calculated ∆G values of -2.394 kcal mol−1 and -2.768 kcal mol−1 respectively. However, because of its additional cationic charge it could be hypothesized that when many R6-1 binding units are work- ing in concert on the surface of HBSPCM, the small improvement in binding affinity they show compared to R4-1 will be amplified and will provide stronger electrostatic interactions with fondaparinux. This will be investigated further in Chapter 5. Also, it should be noted that the ∆G values calculated from the PMF profiles for R4-1 and R6-1 were similar to one obtained from ITC measurements for fondaparinux binding to MPEG-(R4-1), which is a methoxy polyethylene glycol (molecular weight of 350 g mol−1) conjugated with a single R4-1 binding unit developed to study the binding properties of R4-1 to heparin derivatives [91]. 46 4.3. Results and Discussion 4.3.3 Effect of Binding Unit Structure on Free Energy As mentioned above, the geometry and chemical structure of HBSPCM (and its binding units) are of great importance to the HBSPCM-fondaparinux com- plex formation [92]. Therefore, an alternate method of improving binding affinity between binding units and fondaparinux is to change the spacing of the cationic charges within the binding units [92]. As shown in Figure 4.5, shortening the linkage between the N atoms in R4-1 from −CH2 − CH2− 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) [77]. Since the objective of the structural modification is to increase binding affinity, the reduced number of cationic charges in the modified binding unit R4-3 was not desirable; therefore, it 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 −CH2−CH2− (3.84 Å) to −CH2−CH2−CH2− (4.97 Å) to form the R4-2 binding unit (Figure 4.1(d)). R4-2 was thus considered an appropriate structure for further study as it would provide insight into the effect of changing the structure but not the charge of the binding units. Figure 4.5: Structure of the R4-3 binding unit—a modified form of R4-1. The protonation state is for a pH of 7.4 as calculated by the ChemAxon pKa Calculator Plugin (Marvin 5.5.5, 2011) [77]. 47 4.3. Results and Discussion The free-energy profile for R4-2 is shown in Figure 4.6(b). As can be seen in this figure, the free-energy profile for R4-2 displayed a shallow and poorly defined well. Comparing the free-energy profiles of R4-1 and R4-2 suggests that increasing the spacing of cationic charges for this binding unit will not generate an improvement in the binding affinity to fondaparinux. In fact, such a structural change is shown to inhibit fondaparinux binding. (a) (b) Figure 4.6: Free-energy profiles (PMF) for the interaction of fondaparinux with: (a) R4-1, with nitrogen atoms connected by −CH2−CH2− linkages (n= 7) and (b) R4-2, with nitrogen atoms connected by −CH2−CH2−CH2− linkages (n= 8). Calculated values are the Helmholtz free energy which approximates the Gibbs free energy for systems in the condensed phase as explained in Section 2.3.1 and Appendix C. Simulations were performed using a step size of 0.5 Å and a simulation time of 60 ps at each interval. The er- ror bars represent the 95% confidence intervals for n number of replicates. 48 4.3. Results and Discussion 4.3.4 Structural Observations As well as producing free-energy profiles that can be used to determine the binding affinities of fondaparinux to individual binding units, MD simula- tions can also provide insight into the 3D structural configurations of these complexes. For instance, Figure 4.7 shows the structure and orientation of fondaparinux and the binding unit R6-1 at various points along their binding trajectory. Since the MD simulations for R4-1 and R6-1 showed similar trends for the structural configurations of the molecules, only the results for R6-1 are discussed herein. As shown in Figure 4.7, the MD simulations showed a planar structural configuration of the anionic fondaparinux when it was sep- arated from the cationic binding unit and no electrostatic interactions were evident from the free-energy profile—i.e. at centroid-centroid distances larger than 15 Å (Figure 4.4(d) and Figure 4.7(a)–(c)). This is in agreement with reports that have shown fondaparinux to exhibit a planar structure in its an- ionic form [82]. More surprising was the observation from MD simulations of a near planar configuration of fondaparinux when binding to the cationic units (Figure 4.7(d)). It was expected that the interaction of the cationic N atoms of the binding unit with the ionic sites of fondaparinux would pro- duce a bent fondaparinux structure as reported for fondaparinux sodium (i.e. fondaparinux bound to 10 Na+ atoms) [82]. The planar configuration of fon- daparinux observed in MD simulations is hypothesized to be due to the fact that not all 10 anionic sites were interacting with the positively charged atoms of the binding unit R6-1 as was the case for fondaparinux sodium. As the molecules are forced closer and closer together by the restraint on the centroid- centroid distance, more of the anionic sites on fondaparinux interact with the cationic N atoms on R6-1 resulting in a bent configuration of fondaparinux (Figure 4.7(e)–(f)). This, however, was not the most stable configuration for the fondaparinux/(R6-1) system as shown by the large free energies in the PMF at centroid-centroid separations of less than 6 Å (Figure 4.4(d)). In addition to providing details of the structural configuration of fonda- parinux during complex formation, MD simulations also showed that R6-1 49 4.3. Results and Discussion (a) (b) (c) (d) (e) (f) Figure 4.7: MD trajectory for R6-1 and fondaparinux with 2,600 water molecules and 6 Na+ atoms. The structure and orientation of the molecules is shown at centroid-centroid distances of: (a) 20.5 Å, 180 ps; (b) 16.5 Å, 660 ps; (c) 13.5 Å, 1.2 ns; (d) 8.5 Å, 1.62 ns; (e) 6.0 Å, 1.92 ns; and, (f) 1.5 Å, 2.46 ns. Water molecules and Na+ atoms are deleted for clarity. 50 4.4. Conclusions had a high degree of mobility and its carbon backbone angles were constantly changing. This suggests that these binding units are flexible molecules as has been reported for R4-1 [79]. The distance at which the binding affinity between fondaparinux and a binding unit is maximized (thermodynamically favoured) corresponds to the minimum of the free-energy profiles. For both R4-1 and R6-1 the minimum of the free-energy profile occurs when the centroids of the binding unit and fondaparinux are separated by a distance of 8.5 Å (Figure 4.4(d)). This value as well as the information ascertained about the structural configurations of the molecules during binding represent valuable insights that will enable the refinement of the electrostatic binding model in the following chapter. 4.4 Conclusions In silico data for HBSPCM or its binding units has not, up to now, been available. Moreover, there has also been a lack of experimental information on the structure of the HBSPCM-fondaparinux complex. Consequently, the task of optimizing the structure of HBSPCM and its binding units is a daunting challenge. Although time-consuming and very computationally expensive, MD simulations provide a means to gain microscopic information on the interac- tions of HBSPCM binding units and fondaparinux that would have otherwise been inaccessible. As an initial step, MD simulations provided quantitative data that showed that electrostatic interactions are the key contributors to the binding of HBSPCM binding units to fondaparinux. Once this had been determined, the binding affinities of fondaparinux for five cationic binding units were compared using free-energy profiles (PMF). The binding units R4-1 and R6-1 were found to have the deepest free-energy wells suggesting that their binding to fondaparinux is stronger and more stable than that of the other cationic binding units. Additionally, the PMF results indicated that increasing the 51 4.4. Conclusions separation of cationic charges of the binding unit used in the current design of HBSPCM-2 (R4-1) did not improve binding affinity for fondaparinux. Lastly, MD simulations supplied information into the structural configurations exhibited by fondaparinux and the binding units. Specifically, it was found that the complex formed by fondaparinux and a binding unit consisted of a planar fondaparinux molecule at a distance of 8.5 Å from the binding unit and that the molecules were in constant motion relative to one another. The valuable microscopic information gathered from the MD simulations will be used in the next chapter to refine the electrostatic model initially developed in Chapter 3 and, therefore, improve the quantitative model pre- dictions of binding affinity. Once refined, the model will be used to investi- gate HBSPCMs consisting of R4-1 and R6-1 binding units as these molecules were shown to have the most promise for the effective neutralization of fon- daparinux. A model that can accurately predict binding affinities for putative HBSPCMs, without the repeated need for the computationally intensive rou- tines required with MD simulation, is of critical importance as it will greatly streamline the development of an antidote for fondaparinux. 52 Chapter 5 Improved Electrostatic Interactions Model for HBSPCM and Fondaparinux 5.1 Introduction MD simulations were performed in Chapter 4 in order to gain fundamental information of the HBSPCM-fondaparinux interactions. These simulations provided valuable insight that can be used in this chapter, along with experi- mental work, to improve the electrostatic model developed in Chapter 3. This refined model will not only produce more accurate predicted values of ka than were previously calculated but will also represent a useful benchmark against which to compare the simplified model and determine whether its simplifica- tions negatively affect its usefulness. Firstly, the results from the MD simulations presented in Chapter 4 are integrated into the electrostatic model. Experimental ζ-potential results are then used in conjunction with the refined model to determine the appropriate rmin. The predictions of the refined model are then compared to those from the simplified form in order to determine the impact of the latter’s assump- tions. Finally, the model is extended to HBSPCMs with R6-1 binding units to determine if HBSPCM-fondaparinux complex formation can be improved by altering the structure of the cationic binding units. 53 5.2. Materials and Methods 5.2 Materials and Methods 5.2.1 Experimental Charge Determination of Binding Unit R4-1 ζ-Potential measurements were performed at T = 25oC using a Coulter® Delsa 440SX Zeta Potential Analyzer (Coulter Co. Instrument Co.) equipped with a 5 mW helium-neon laser operating at a wavelength of 632 nm. The fondaparinux stock solution was prepared by diluting commercially available ARIXTRA® (GSK) in 5 mM NaCl (ACS grade, Fisher Scientific) to a con- centration of 1 mg mL−1. The HBSPCM stock solution was a dilution of the HBSPCM-2 (synthesized by the Kizhakkedathu Laboratory in the De- partment of Pathology and Laboratory Medicine at the University of British Columbia) in 5 mM NaCl (ACS grade, Fisher Scientific) to a concentration of 5 mg mL−1. The conductivity standard used was Traceable Conductivity Stan- dard (Fisher Scientific) and the mobility standard was Aqueous Suspension of Polystyrene Latex (Otsuka Electronics Co. Ltd.). At pH 7.0, the ζ-potential of the HBSPCM-fondaparinux complex was measured at HBSPCM to fonda- parinux concentration ratios of: 0.0, 1.7, 2.5, 5.0, 10.0, 15.0, and 40.0. The reported results represent the average of 24 measured values at each concen- tration. 5.2.2 Refinement of Electrostatic Model to Optimize HBSPCM Structure The basic structure of the electrostatic model is described fully in Section 3.2. The model is herein refined based on experimental work as well as on the mi- croscopic structural and binding characteristics of the HBSPCM binding units and fondaparinux found in Chapter 4. In the simplified form of the model developed in Chapter 3, the binding units (R4-1) were represented as having a charge of +3, as this was the charge that was calculated for a free binding unit at physiological conditions using the ChemAxon pKa Calculator Plugin (Marvin 5.5.5, 2011) [77]. However, the results of ζ-potential measurements 54 5.2. Materials and Methods discussed in Section 5.3.1 showed that these units have an effective charge of +1.29. Therefore, the model was altered to utilize this value for the charge of the R4-1 binding units. As a result, the refined model generates one HBSPCM and allows 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 HBSPCM are there- fore available for binding to another HBSPCM in solution. Consequently, the model now considers system neutralization rather than fondaparinux neutral- ization and will implicitly account for bridging effects, except for any associ- ated structural modifications that might occur. In the simplified model, the minimal distance of approach, a, was set to 6 Å since this value had been shown to give the best fit for experimental data collected for a variety of protein systems [44]. However, the free-energy pro- files of the interaction of binding units R4-1 and R6-1 to fondaparinux (Figure 4.4(c) and Figure 4.4(d)) showed that the centroid-centroid distance between a binding unit and fondaparinux was 8.5 Å at binding (Section 4.3.4). Since the refined model focused on these binding units, the value for a was changed to 8.5 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 units that formed a binding site. However, the MD simulation results presented in Chapter 4 suggest that during the binding of fondaparinux to either R4- 1 or R6-1 the ionic sites of fondaparinux remain at a certain distance from the protonated N atoms of the binding units. In addition, it was found that the molecules in a complex show significant and constant relative motion to one another. Therefore, the model was refined to incorporate a more realistic binding criterion. A binding site was defined as three or more binding units that are within an area that allows them to come within a prescribed distance (8.5 Å) of a fondaparinux molecule which is centered around one of the binding units. 55 5.3. Results and Discussion As with the simplified model, the Debye-Hückel energy of the complex in the refined model was found by considering a bound fondaparinux molecule and its associated binding units to be one binding site. Upon binding, a single charge was then assigned to this binding site (the sum of the charges of fonda- parinux and the associated binding groups) at the location of the first binding unit that formed the binding site. This was due to the geometrical consider- ations of the new binding criterion in which a randomly selected binding unit is chosen as the starting point to look for other binding units that are within the required radius. The MD simulation results showed that the R4-1 binding unit was flexi- ble in nature (Section 4.3.4) as has been reported in the literature [79] and as was proved by calculations performed in Section 3.3.1. The R6-1 bind- ing unit was also found to exhibit similar flexibility during the MD studies (Section 4.3.4). Therefore, binding unit flexibility was incorporated into the refined model a priori. Simulation calculations were performed in the same manner as in Chapter 3 except for the modifications explained in this section. The MATLAB code for the model programs used in this chapter is found in Appendix D. 5.3 Results and Discussion 5.3.1 Effective Charge of Binding Unit R4-1 from ζ-Potential Measurements The ζ-potential for the complex formed between HBSPCM-2 and fondaparinux was measured and is shown in Figure 5.1. As expected, the ζ-potential in- creased with the addition of HBSPCM-2 to the fondaparinux solution. Since the magnitude of the ζ-potential is directly correlated to the complex stability [93], a ζ-potential greater than or equal to zero represents a complete electro- static charge neutralization of fondaparinux by HBSPCM. By interpolation, 56 5.3. Results and Discussion the molar concentration ratio that produces a ζ-potential of 0 mV was deter- mined. This value was then used to yield an effective charge for the individual binding units on HBSPCM-2 which was found to be +1.29. This is lower than the charge predicted by the ChemAxon pKa Calculator Plugin (Marvin 5.5.5, 2011) [77] and suggests that to neutralize fondaparinux, more binding unit charges are required than would be prescribed theoretically. This may be the result of the attachment of the binding units to the HBSPCM core. For instance, the protonation state of the binding units may be altered by the attachment and/or the PEG layer of the HBSPCM may shield some of the N atoms within the binding units. In either case, the electrostatic model was refined to take into account the effective charge of the binding units. Figure 5.1: ζ-Potential as a function of molar concentration ratio for mix- tures of HBSPCM-2 and fondaparinux. The concentration of fon- daparinux was kept constant at 1 mg mL−1. Measurements were performed at T = 25oC in a Coulter® Delsa 440SX Zeta Poten- tial Analyzer (Coulter Co. Instrument Co.) equipped with a 5 mW helium-neon laser operated at a wavelength of 632 nm. ζ-Potential is the arithmetic mean of 24 measured values. The error bars rep- resent the 95% confidence interval. 57 5.3. Results and Discussion 5.3.2 Model Predictions for HBSPCM with R4-1 as Binding Unit Characterization of Binding Between HBSPCM and Fondaparinux Following the same procedure as in Section 3.3.1, the results from the ζ- potential measurements (Section 5.3.1) were used to determine the minimum distance between binding units, rmin, for HBSPCMs with R4-1 as binding units. From the simulation results shown in Figure 5.2 the new value of rmin was found to be 9.16 Å, which is larger than the radius of R4-1 (6.61 Å) and R6-1 (7.10 Å). Since the purpose of introducing an rmin was to account for electrostatic repulsions and steric interactions between binding units, this new rmin is a more realistic constraint than the previous value which was smaller than the radius of R4-1. Figure 5.2: Computer simulated number of fondaparinux molecules bound per HBSPCM for HBSPCMs with 3 to 20 attached R4-1 bind- ing units. The HBSPCMs had an rH of 4.0 nm and an rmin of 9.16 Å. Each simulation was repeated 1,000 times; the resulting arithmetic means are shown. 58 5.3. Results and Discussion Effect of Number of Binding Units on ka and Number of Molecules of Fondaparinux Bound per HBSPCM After determining the appropriate rmin the improved model was used, as in Section 3.3.2, to determine the effect of number of binding units on the binding of HBSPCM to fondaparinux. Once again ka is used as a metric for binding affinity since, as previously described in Section 2.2.1, kd is unaffected by elec- trostatic interactions. Results from the refined model showed similar trends to those seen in the simplified form. For instance, as seen in Figure 5.3, ka and the number of fondaparinux bound per HBSPCM increased with the number of attached binding units in similar fashions to what was observed using the sim- plified model (Figure 3.5). In addition, both the simplified and refined models predicted similar values for the number of fondaparinux molecules bound to HBSPCM. However, the range for ka was reduced by 6 orders of magnitude in 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 anti- dote discovery, the refined model with the aid of MD simulation results yields more quantitatively accurate information. While the ka values produced by the refined model are more realistic than those previously reported in Chap- ter 3, they may still be slightly overestimated compared to known values for protein-protein interactions [43] as the improved model does not take into ac- count any repulsion that could exist between a fondaparinux already bound to HBSPCM and an incoming fondaparinux molecule. 59 5.3. Results and Discussion Figure 5.3: Computer simulated ka and number of fondaparinux molecules bound per HBSPCM for HBSPCMs with 3 to 20 attached R4-1 binding units. The HBSPCMs had an rH of 4.0 nm and an rmin of 9.16 Å. ka and number of fondaparinux molecules bound to an HBSPCM are the geometric and arithmetic means respectively of 1,000 calculated values. The error bars represent the 95% confi- dence intervals. Effect of HBSPCM Size on ka and Number of Molecules of Fondaparinux Bound per HBSPCM The effect of number of binding units and HBSPCM size on the number of fondaparinux molecules bound per HBSPCM and on ka was investigated with the improved binding model and is shown in Figure 5.4 and Figure 5.5. As with the simplified model, it was observed that as the surface area of the HBSPCM core increases, the probability of finding binding units sufficiently close to each other to form a binding site decreases. By comparing the num- ber of fondaparinux molecules bound per HBSPCM predicted by the simplified model (Figure 3.4) and that predicted by the improved model (Figure 5.4), it 60 5.3. Results and Discussion can be observed that the improvements to the model have a greater impact for smaller HBSPCMs. As shown in Figure 5.4, the refined model does not predict a linear increase in the number of fondaparinux molecules bound per HBSPCM with a radius of 2 nm as was the case for the simplified model (Figure 3.4). Instead, the number of fondaparinux molecules bound per HBSPCM starts to plateau at a high number of binding units. Similar trends were observed for ka. These results are to be expected as the number of fondaparinux bound to an HBSPCM of a given surface area will eventually reach a maximum due to space limitations. Figure 5.4: Computer simulated number of fondaparinux molecules bound per HBSPCM for HBSPCMs with radii, rH , of 2.0 nm, 4.0 nm, 6.0 nm, and 10.0 nm. All different sized HBSPCMs were tested with a number of R4-1 binding units ranging from 3 to 20 using an rmin of 9.16 Å. The number of fondaparinux molecules bound to an HBSPCM is the arithmetic mean of 1,000 calculated values. The error bars represent the 95% confidence interval. 61 5.3. Results and Discussion Figure 5.5: Computer simulated ka for HBSPCMs with radii, rH , of 2.0 nm, 4.0 nm, 6.0 nm, and 10.0 nm. All different sized HBSPCMs were tested with a number of R4-1 binding units ranging from 3 to 20 using an rmin of 9.16 Å. ka is the geometric mean of 1,000 calculated values. The error bars represent the 95% confidence interval. 5.3.3 Model Extension for HBSPCM with R6-1 as Binding Unit The MD simulations from Chapter 4 indicated that the R6-1 binding unit had a slightly lower ∆G for binding than R4-1. It was argued that this slight difference might be amplified when multiple binding units are working together to bind fondaparinux. Therefore, the improved model was used to explore the HBSPCM-fondaparinux complex formation with R6-1 binding units in place of R4-1. In order to use the model with R6-1, the appropriate effective charge was required. Due to the time and effort associated with the synthesis of HBSPCMs, experimental data for this value was not available. However, the model allows for predictions of ka and of the number of fondaparinux molecules bound per HBSPCM to be made for a range of effective charges. Because R6-1 62 5.3. Results and Discussion has a theoretical charge +1 higher than R4-1, the effective charges investigated herein ranged up to +2.29 (+1 higher than the effective charge of R4-1). It was assumed that the rmin of 9.16 Å required for the R4-1 binding units was a good approximation for R6-1 since both binding units bind to fondaparinux at the same centroid-centroid distance and have a similar size. Impact of Binding Unit Effective Charge on ka and Number of Molecules of Fondaparinux Bound per HBSPCM The impact of the effective charge of the R6-1 binding units on ka and on the number of molecules of fondaparinux bound per HBSPCM was investigated for varying sizes of HBSPCMs. For the purposes of this investigation, HBSPCMs with 20 R6-1 binding units on their surface were used because the effects being studied are more easily observed at a high number of binding units. As seen in Figure 5.6, the number of fondaparinux molecules bound per HBSPCM was not dependent on the effective charge of the binding units. This was to be ex- pected as the charge of the binding units does not affect the model’s placement of these units on the surface of HBSPCM nor the binding criterion and, thus, will have no impact on whether or not an individual fondaparinux molecule will find an appropriate binding site. Conversely, as explained in Section 2.2.1, ka is dependent on electrostatic forces and can be increased by increasing the magnitude of these forces [44, 48, 49]. Therefore, it was expected that ka would increase as the effective charge of the binding units was increased (Figure 5.7). The results show that the in- crease in ka with increasing binding unit effective charge was most pronounced at an rH of 2.0 nm. For this size of HBSPCM, ka increased from 10 17 M−1 s−1 to 1024 M−1 s−1 at effective charges of +0.1 and +2.29, respectively. For HB- SPCMs 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 units be- come spaced farther apart with increasing HBSPCM radius, the impact of any increase in effective charge becomes dampened. 63 5.3. Results and Discussion Figure 5.6: Computer simulated number of fondaparinux molecules bound per HBSPCM for HBSPCMs with radii, rH , of 2.0 nm, 4.0 nm, 6.0 nm, and 10.0 nm. All different sized HBSPCMs were tested with 20 R6-1 binding units with an effective charge ranging from +0.1 to +2.29 using an rmin of 9.16 Å. The number of fondaparinux molecules bound to an HBSPCM is the arithmetic mean of 1,000 calculated values. The error bars represent the 95% confidence interval. 64 5.3. Results and Discussion Figure 5.7: Computer simulated ka for HBSPCMs with radii, rH , of 2.0 nm, 4.0 nm, 6.0 nm, and 10.0 nm. All different sized HBSPCMs were tested with 20 R6-1 binding units with an effective charge ranging from +0.1 to +2.29 using an rmin of 9.16 Å. ka is the geometric mean of 1,000 calculated values. The error bars represent the 95% confidence interval. 65 5.3. Results and Discussion Effect of Binding Units on ka and Number of Molecules of Fondaparinux Bound per HBSPCM The improved model can also be used to determine the effect of replacing R4-1 binding units with R6-1. As mentioned, the effective charge of the R6-1 binding units is not known; therefore, model predictions for ka were made with R6-1 effective charges of +1.29 and +2.29 (Figure 5.8(a) and Figure 5.8(b)). These values represent an effective charge that is equivalent to that of R4-1 and one that is increased by +1. As shown in Figure 5.8(a), when R6-1 is assumed to have an equivalent effective charge to that of R4-1 the model predicts no difference in the calculated values of ka. This indicates that the slightly larger radius of R6-1 does not impact the formation of binding sites to a degree that an effect is observed in the affinity of HBSPCM to fondaparinux. However, the slightly larger R6-1 binding units would most likely experience reduced charge shielding associated with the PEG layer of HBSPCM. Therefore, it is likely that the effective charge of the R6-1 binding units would be higher than that of R4-1. The effect of an increased R6-1 effective charge can be seen in Figure 5.8(b) which shows increased ka values with R6-1 compared to R4-1, especially for HBSPCMs with an rH of 2.0 nm and 4.0 nm. This result suggests that the binding of HBSPCM to fondaparinux can be improved by using R6-1 binding units on small HBSPCM cores. This is consistent with the results of the simplified model, which recommended the synthesis of HBSPCMs with low hydrodynamic size in order to increase the charge density of HBSPCM and, therefore, enhance ka. 66 5.3. Results and Discussion (a) (b) Figure 5.8: Comparison of computer simulated ka values for HBSPCMs con- sisting of R4-1 binding units and R6-1 binding units with effective charges of: (a) +1.29 and (b) +2.29. HBSPCMs of radii, rH , of 2.0 nm, 4.0 nm, 6.0 nm, and 10.0 nm were tested with the number of binding units ranging from 3 to 20 and using an rmin of 9.16 Å. ka is the geometric mean of 1,000 calculated values. The error bars represent the 95% confidence interval. 67 5.4. Conclusions 5.4 Conclusions The electrostatic interaction model developed in Chapter 3 provided valuable insight into the effect of HBSPCM structure on fondaparinux binding affinity and indicated that improved fondaparinux neutralization could be achieved by increasing the HBSPCM’s charge density. Because of a lack of knowledge of the microscopic structural characteristics of fondaparinux and HBSPCM this model, however, relied on certain assumptions which caused it to quantitatively overpredict ka values. In an attempt to improve accuracy, this simplified model was refined in this chapter by incorporating the results of MD simulations of the interactions between individual cationic binding units and fondaparinux. The refined model was first used in conjunction with ζ-potential experi- ments to determine the effective charge of the R4-1 binding units attached to the HBSPCM surface. This effective charge was found to be lower than that of the isolated binding unit suggesting that the attachment of the binding units reduced or limited the accessibility to the charged atoms. With an appropriate effective charge defined, the predictions of the simplified and refined forms of the model were compared. It was found that the qualitative trends in complex formation predicted by the simplified model closely approximated those from the more refined form. However, the values of ka predicted by the refined model were found to be more realistic, most significantly for small HBSPCM radii. The refined electrostatic model was then extended to HBSPCMs consist- ing of R6-1 binding units. Although the effective charge of these binding units was not known, different hypothetical charges based on MD observations of R6-1 were used with the model to test the effect of polymer size and num- ber of binding units on the binding of HBSPCM and fondaparinux. It was found that if, as would be expected, the R6-1 binding units have a larger ef- fective charge than the R4-1 binding units, improved binding affinity between HBSPCM and fondaparinux is seen. The results of this chapter therefore in- 68 5.4. Conclusions dicate that increasing the charge density of HBSPCM with favourable binding units will improve complex formation between HBSPCM and fondaparinux. Taking into account physiological considerations such as toxicity and clear- ance, it is recommended to synthesize small HBSPCMs containing as many R6-1 binding units as possible. Although the refined model has been shown in this chapter to predict more realistic values for ka than the simplified form, both models provided the same recommendations concerning HBSPCM candidate structures. The initial model thus represents, even though it relies on a simplified concept of the molecular structures involved in binding, a useful tool for the study of HBSPCM-fondaparinux complex formation in the absence of MD simulation data. Therefore, the electrostatic model, both in its simplified and refined forms, produces guidelines on HBSPCM candidate selection that are critical for the rapid and efficient discovery of an effective antidote for fondaparinux. 69 Chapter 6 Conclusions As discussed in Chapter 1, anticoagulants are widely used for the treatment of a variety of blood-clot related illnesses such as venous thromboembolism, unstable angina, and acute myocardial infarction [9–11]. They are also ad- ministered during different procedures such as cardiac [4, 9] and orthopedic surgeries [11, 24] and kidney dialysis [4]. The most commonly used antico- agulants are heparin derived drugs [4] which include UFH, LMWHs, and the synthetic pentasaccharide derivatives fondaparinux and idraparinux [2, 5–7]. Fondaparinux in particular shows much promise as an anticoagulant because of its advantages over current therapies such as UFH and LMWHs; however, the lack of a specific antidote for fondaparinux is a great impediment to its widespread use since the inability to neutralize anticoagulants can lead to hemorrhaging [3, 5, 7, 12, 15]. There is therefore an urgent need to develop an antidote that will effectively neutralize fondaparinux [14]. This will not only reduce the health risks for patients but it will also help reduce the costs associated with blood transfusions and longer hospital stays [13]. A promising antidote for fondaparinux, HBSPCM, is currently under de- velopment; however, a molecular structure still needs to be found that will reverse the anticoagulant effect of fondaparinux [18]. Traditionally, the opti- mization process for a new drug is both time-consuming and labour intensive as it requires the synthesis and testing of thousands of new molecules [17]. Both costs and time associated with drug discovery can be greatly reduced with the implementation of computer simulations that allow for the rapid testing of multiple candidate molecules in silico [17]. The objective of this thesis therefore was to construct a framework of modeling tools that would allow for the characterization and refinement of the HBSPCM-fondaparinux 70 Chapter 6. Conclusions complex formation in order to streamline the discovery process of an HBSPCM structure that will bind to fondaparinux with high affinity and neutralize it. Firstly, an empirical electrostatic model was developed and used to char- acterize the binding of HBSPCM to fondaparinux. This model predicted the binding affinities for different HBSPCM structures and indicated that im- proved fondaparinux neutralization could be achieved with small HBSPCMs containing as many cationic binding units as possible. However, because of the lack of experimental data on the interaction between HBSPCMs and fonda- parinux, the model relied on certain assumptions. These simplifications caused the model to quantitatively overpredict ka values. In an attempt to improve the accuracy of this simplified model, MD simulations were performed. The information garnered on the microscopic structure of the complex formed by fondaparinux and a binding unit was then incorporated into the model. This refined model not only produced more realistic values for ka but also repre- sented a useful benchmark with which to demonstrate that the simplifications made in the initial model did not negatively affect its usability. In addition, the refined model was extended to HBSPCMs that contained a more posi- tively charged cationic binding unit. The model results, like those from the simplified form, showed that the ka between HBSPCM and fondaparinux can be enhanced by increasing charge density on HBSPCMs of small radii. With the aid of the model developed in this thesis, Kizhakkedathu and co-workers have recently found an HBSPCM structure that initial results have shown to effectively neutralize fondaparinux in vitro. Although the specific molecular structure of this HBSPCM cannot be disclosed for proprietary rea- sons, as predicted by the electrostatic model described in this thesis, charge density was crucial in its development. While outside the scope of this thesis, the experimental characterization of this or other candidate HBSPCMs will involve in vitro cytotoxicity and efficacy studies followed by in vivo studies. These tests include activated partial thromboplastin time (aPTT), thromboe- lastography (TEG), anti-Xa assay, reversal of inhibition of thrombin activity, 71 Chapter 6. Conclusions complement activation, and platelet aggregation and activation by HBSPCM and/or HBSPCM-fondaparinux complexes [18]. At the time of printing this thesis, in vitro studies were being performed on the aforementioned candidate HBSPCM. The tentative discovery of an effective HBSPCM structure has proven that the combined approach of computer simulations and experimental validation reduces the time and expense associated with the conventional trial-and-error design approach for drug discovery. Specifically, the development, in this the- sis, of an electrostatic model to refine the structure of HBSPCM was an es- sential first step towards the discovery of an antidote for fondaparinux. This model, both in its simplified (Chapter 3) and refined (Chapter 5) form, pro- vided guidelines for the selection of an HBSPCM candidate which were critical for the rapid and cost efficient discovery of an effective antidote for fonda- parinux. The strengths of the model are twofold: a) minimal experimental data is required—data that is difficult to obtain for HBSPCMs [18]—and, b) the HBSPCM structural adjustments can be determined in silico thus avoid- ing the costly and time-intensive synthesis of new HBSPCMs. The research presented in this thesis thus represents a major contribution to the process of finding an antidote for fondaparinux and, therefore, will greatly impact the therapeutic field. Although the model presented in this thesis represents a useful framework to streamline the discovery process for an effective antidote for fondaparinux, modeling cannot fully replace experimental data and there still exists a lack of experimental results in the literature for the binding of fondaparinux and HBSPCM. While outside the scope and resources of this thesis project, exper- imental measurement of ka and kd for different HBSPCMs and fondaparinux would be useful not only for the therapeutic development of HBSPCM but also to determine the quantitative accuracy of the model predictions. One method of experimentally determining rates of association is surface plasmon resonance (SPR). In this technique, the candidate HBSPCM would be immo- 72 Chapter 6. Conclusions bilized on the surface of a chip over which a solution of fondaparinux would be passed. As binding occurs, the mass and therefore the refractive index of the solution immediately above the chip changes and is registered in a real-time sensorgram from which ka and kd can be determined [4, 94]. Furthermore, us- ing ITC technology to experimentally determine Ka for the reaction between each candidate HBSPCM and fondaparinux would be useful to assess the va- lidity of the model predictions. In this technique, the candidate HBSPCM would be injected in a stepwise manner into a solution of fondaparinux and the heat released on binding after each injection would be measured. From these measurements, Ka can then be obtained [4]. This would allow to com- pare the experimentally measured Ka trends to the ka trends predicted by the model for different HBSPCMs. In addition to experimental work, another area of interest concerning the interaction between HBSPCM and fondaparinux would be further MD simu- lations. For instance, the electrostatic model did not take into account any re- pulsion or attraction effects that could exist between a fondaparinux molecule already bound to HBSPCM and an incoming fondaparinux molecule since the MD studies presented in this thesis were, out of necessity, performed for one to one interactions of a binding unit and a fondaparinux molecule. To further investigate the behaviour of multiple fondaparinux molecules binding to an HBSPCM, free-energy profiles (PMF) of systems at different concentrations of binding units and fondaparinux could be carried out. The results of these MD simulations could then be used to further refine, if needed, the electro- static model. As mentioned above, a promising HBSPCM has been found to effectively neutralize fondaparinux. However, if the in vitro or in vivo experiments per- formed on this HBSPCM were to reveal its ineffectiveness, the electrostatic model would be required to screen for another promising molecular structure. Most likely, this would entail the study of additional binding units. In order to make the electrostatic model applicable to more binding units, ζ-potential 73 Chapter 6. Conclusions experiments will be required to determine the effective charge of at least two other binding units. These charges will allow for the generation of an effective charge calibration curve that can be incorporated into the model to appropri- ately characterize the binding of fondaparinux to different binding units. The electrostatic model also shows promise for applications other than fondaparinux antidote discovery. In fact, the model could be useful to describe any complex formation that is mediated primarily by electrostatic interactions between molecules. One such application is the delivery of therapeutic genes to specific targets, which still remains a challenge for researchers [92]. A promising approach to achieve successful gene transfer is the use of DNA- polycation complexes (or polyplexes). Polyplexes are formed as a result of the electrostatic interactions between the negatively charged groups in DNA and the positively charged groups of a polymer [95–98]. The spacing of the cationic charges and the number of these charges can determine the size of the polyplexes, which in turn could play a key role in the efficiency of gene transfer [95]. Because of the electrostatic nature of these polyplexes, the mathematical model developed in this thesis can potentially be used to aid in their successful design. 74 Bibliography [1] Milan Remko. Molecular structure, lipophilicity, solubility, absorption, and polar surface area of novel anticoagulant agents. Journal of Molec- ular Structure: THEOCHEM, 916:76–85, 2009. [2] Jeffrey I. Weitz. Emerging anticoagulants for the treatment of venous thromboembolism. Journal of Thrombosis and Haemostasis, 96:274–284, 2006. [3] Tracy Minichiello and Patrick F. Fogarty. Diagnosis and management of venous thromboembolism. The Medical Clinics of North America, 92(2):443–465, 2008. [4] Ishan Capila and Robert J. Lindhardt. Heparin-protein interactions. Angewandte Chemie International Edition, 41(3):390–412, Feb. 2002. [5] Jeffrey I. Weitz and Lori-Ann Linkins. Beyond heparin and warfarin: The new generation of anticoagulants. Expert Opinion on Investigational Drugs, 16(3):271–282, 2007. [6] Lori-Ann Linkins and Jeffrey I. Weitz. New anticoagulants. Seminars in Thrombosis and Hemostasis, 29(6):619–631, 2003. [7] S. Alban. From heparins to factor Xa inhibitors and beyond. European Journal of Clinical Investigation, 35(1):12–20, 2005. [8] John F. Canales and James J. Ferguson. Low-molecular-weight heparins: Mechanisms, trials, and role in contemporary interventional medicine. American Journal Of Cardiovascular Drugs: Drugs, Devices, And Other Interventions, 8(1):15–25, 2008. 75 Bibliography [9] Jack Hirsh, Stephen G. Shaughenessy, Jonathan L. Halperin, Christo- pher Granger, E. Magnus Oluman, and James E. Dalen. Heparin and low-molecular-weight heparin: Mechanism of action, pharamcokinetics, dosing, monitoring, efficacy, and safety. Chest, 119:64S–94S, 2001. [10] Michael D. Freedman. Pharmacokinetics, clinical indications, and ad- verse effects of heparin. The Journal of Clinical Pharmacology, 32:584– 596, 1992. [11] Jeffrey I. Weitz. Low-molecular-weight heparins. The New England Journal of Medicine, 337(10):688–698, 1997. [12] Sam Schulman and Nick R. Bijsterveld. Anticoagulants and their rever- sal. Transfusion Medicine Reviews, 21(1):37–48, Jan. 2007. [13] Jerrold H. Levy and Kenichi A. Tanaka. The anticoagulated patient: Strategies for effective blood loss management. Surgery, 142:S71–S77, Oct. 2007. [14] Brian K. Adler. Unfractionated heparin and other antithrombin medi- ated anticoagulants. Clinical Laboratory Science, 17(2):113–117, 2004. [15] J. Ingerslev, T. Vanek, and S. Culic. Use of recombinant factor viia for emergency reversal of anticoagulation. Journal of Postgraduate Medicine, 53(1):17–22, 2007. [16] Roderick E. Hubbard. Structure-Based Drug Discovery, chapter 1. RSC Publishing, 2006. [17] C. L. Propst and Thomas J. Perun. Computer-Aided Drug Design, chap- ter 1. Marcel Dekker, Inc., 1989. [18] Jayachandran N. Kizhakkedathu. Personal Communications, 2008. Uni- versity of British Columbia, Vancouver, BC, Canada. [19] Jack Hirsh. Current anticoagulant therapy – unmet clinical needs. Thrombosis Research, 109:S1–S8, 2003. 76 Bibliography [20] Lawrence L.K. Leung. Overview of hemostasis. www.uptodate.com, November 2009. [21] I. Björk and U. Lindhal. Mechanism of the anticoagulant action of hep- arin. Molecular and Cellular Biochemistry, 48:161–182, 1982. [22] Shannon M. Bates and Jeffrey I. Weitz. Emerging anticoagulant drugs. Arteriosclerosis, Thrombosis, and Vascular Biology, 23:1491–1500, 2003. [23] D.J. Perry. Antithrombin and its inherited deficiencies. Blood Reviews, 8:37–55, 1994. [24] Nicole T. Ansani. Fondaparinux: The first pentasaccharide anticoagu- lant. P & T, 27(6):310–317, Jun. 2002. [25] Earl W. Davie. Biochemical and molecular aspects of the coagulation cascade. Thrombosis and Haemostasis, 71(1):1–6, 1995. [26] Umesh R. Desai, Maurice Petitou, Ingemar Björk, and Steven T. Ol- son. Mechanism of heparin activation of antithrombin. The Journal of Biological Chemistry, 273(13):7478–7487, Mar. 1998. [27] Lei Jin, Jan Pieter Abrahams, Richard Skinner, Maurice Petitou, and Robert N. Pike. The anticoagulant activation of antithrombin by hep- arin. Proceedings of the National Academy of Sciences of USA, 94:14683– 14688, Dec. 1997. [28] R. J. Linhardt and N. S. Gunay. Production and chemical processing of low molecular weight heparins. Seminars in Thrombosis and Hemostasis, 25(3):5–16, 1999. [29] Alfonso Bentolila, Israel Vlodavsky, Christine Haloun, and Abraham J. Domb. Synthesis and heparin-like biological activity of amino acid-based polymers. Polymers for Advanced Technologies, 11(8–12):377–387, 2000. [30] Elaine Gray, Barbara Mulloy, and Trevor W. Barrowcliffe. Heparin and low-molecular-weight heparin. Thrombosis Haemostasis, 99:807– 818, 2008. 77 Bibliography [31] Kenneth A. Bauer. New anticoagulants. Hematology, pages 450–456, Jan. 2006. [32] Jack Hirsh. New anticoagulants. American Heart Journal, 142(2):S3–S8, Aug. 2001. [33] Mark A. Crowther and Theodore E. Warkentin. Bleeding risk and the management of bleeding complications in patients undergoing anticoagu- lant therapy: Focus on new anticoagulant agents. Blood, 111:4871–4879, 2008. [34] Jan Charles Horrow. Protamine: A review of its toxicity. Anesthesia & Analgesia, 64:348–361, 1985. [35] Mark A. Crowther, Leslie R. Berry, Paul T. Monagle, and Anthony K. Chan. Mechanisms responsible for the failure of protamine to inactivate low-molecular-weight heparin. British Journal of Haematology, 116:178– 186, 2002. [36] Brian W. Pontius. Close encounters: Why unstructured, polymeric do- mains can increase rates of specific macromolecules association. TIBS, 18:181–186, May 1993. [37] Andrew K Powell, Edwin A. Yates, David G. Fernig, and Jeremy E. Turnbull. Interactions of heparin/heparan sulfate with proteins: Ap- praisal of structural factors and experimental approaches. Glycobiology, 14(4):17R–30R, Jan. 2004. [38] Gregory V. Nikiforovich. Computational molecular modeling in pep- tide drug design. International Journal of Peptide & Protein Research, 44:513–531, 1994. [39] Wolfgang Bitomsky and Rebecca C. Wade. Docking of glycosamino- glycans to heparin-binding proteins: Validation for afgf, bfgf, and an- tithrombin and application to il-8. Journal of the American Chemical Society, pages 3004–3013, 1999. 78 Bibliography [40] Hanah Margalit, Nurit Fischer, and Shmuel A. Ben-Sasson. Compar- ative analysis of structurally defined heparin binding sequences reveals a distinct spatial distribution of basic residues. Journal of Biological Chemistry, 268(26):19228–19231, Sept. 1993. [41] Ruth Tyler-Cross, Michael Sobel, Louis E. McAdory, and Robert B. Har- ris. Structure-function relations of antithrombin iii-heparin interactions as assessed by biophysical and biological assays and molecular modeling of peptide-pentasaccharide-docked complexes. Archives of Biochemistry and Biophysics, 334(2):206–213, Oct. 1996. [42] Gideon Schreiber and Alan Fersht R. Rapid, electrostatically assisted association of proteins. Nature Structural Biology, 3(5):427–431, May 1996. [43] Tzvia Selzer and Gideon Schreiber. Predicting the rate enhancement of protein complex formation from the electrostatic energy of interaction. Journal of Molecular Biology, 287:409–419, 1999. [44] Tzvia Selzer, Shira Albeck, and Gideon Schreiber. Rational design of faster associating and tighter binding protein complexes. Nature Struc- tural Biology, 7(7):537–541, Jul. 2000. [45] Giovanna Scapin. Structural biology and drug discovery. Current Phar- maceutical Design, 12(17):2087–2097, 2006. [46] Lachelle Arnt and Gregory N. Tew. New poly(phenyleneethynylene)s with cationic, facially amphiphilic structures. Journal of the American Chemical Society, 124:7664–7665, 2002. [47] Tommy Liljefors and Ingrid Petterson. Textbook of Drug Design, chap- ter 4. Taylor & Francis, 2002. [48] Gideon Schreiber, Gilad Haran, and Huan-Xiang Zhou. Fundamen- tal aspects of protein-protein association kinetics. Chemical Reviews, 109(3):839–860, 2009. 79 Bibliography [49] Tzvia Selzer and Gideon Schreiber. New insights into the mechanism of protein-protein association. Proteins: Structure, Function, and Genet- ics, 45:190–198, 2001. [50] Otto G. Berg and Peter H. von Hippel. Diffusion-controlled macro- molecular interactions. Annual Review of Biophysics and Biophysical Chemistry, 14:131–160, 1985. [51] Huan-Xiang Zhou. Disparate ionic-strength dependencies of on and off rates in protein-protein association. Biopolymers, 59:427–433, 2001. [52] M. Vijayakumar, Kwan-Yin Wong, Gideon Schreiber, Alan R. Fersht, Attila Szabo, and Huan-Xiang Zhou. Electrostatic enhancement of diffusion-controlled protein-protein association: Comparison of theory and experiment on barnase and barstar. Journal of Molecular Biology, 278:1015–1024, 1998. [53] Gideon Schreiber, Yosi Shaul, and Kay E. Gottschalk. Protein Design: Methods and Applications, chapter 11, pages 250–264. Humana Press, 2006. [54] J.TH.G. Overbeek and B.H. Bijsterbosch. Electrokinetic Separation Methods, chapter The electrical double layer and the theory of elec- trophoresis, pages 1–33. Elsevier/North-Holland Biomedical Press, 1979. [55] Freddie R. Salsbury Jr. Molecular dynamics simulations of protein dy- namics and their relevance to drug discovery. Current Opinion in Phar- macology, 10:1–7, 2010. [56] Michael P. Allen. Introduction to molecular dynamics simulation. Com- putational Soft Matter: From Synthetic Polymers to Proteins, 23:1–28, 2004. [57] Martin Karplus and J. Andrew McCammon. Molecular dynamics sim- ulations of biomolecules. Nature Structural Biology, 9(9):646–652, 2002. Review. 80 Bibliography [58] Roland Stote, Annick Dejaegere, Dmitry Kuznetsov, and Lau- rent Falquet. Molecular dynamics simulations - CHARMM. www.ch.embnet.org/MD tutorial/, 1999. [59] Robert J. Woods and Matthew B. Tessier. Computational glycoscience: Characterizing the spatial and temporal properties of glycans and glyca- protein complexes. Current Opinion in Structural Biology, 20:575–583, 2010. [60] Accelrys Software Inc., San Diego, USA. Materials Studio® Online Help, Release 5.5, 5.5 edition, 2010. [61] Wilfred van Gunsteren, Dirk Bakowies, Bürgi, Indira Chandrasekhar, Markus Christen, Xavier Daura, Peter Gee, Alice Glaätti, tomas Hans- son, Chris Oostenbrink, Christine Peter, Jed Pitera, Lukas Schuler, Thereza Soares, and Haibo Yu. Molecular dynamics simulation of biomolecular systems. CHIMIA, 55:856–860, 2001. [62] Olgun Guvench and Alexander D. Jr. MacKerell. Comparison of protein force fields for molecular dynamics simulations. Methods in Molecular Biology, 443:63–88, 2008. [63] F. Sato and S. Hojo. On the transferability of force field parameters - with an ab initio force field developed for sulfonamides. The Journal of Physical Chemistry A, 107:248–257, 2003. [64] H. Sun. COMPASS: An ab initio forcefield optimized for condensed- phase applications - overview with details on alkane and benzene com- pounds. Journal of Physical Chemistry B, 102:7338, 1998. [65] M. P. Allen and D. J. Tildesley. Computer Simulation of Liquids. Oxford University Press, 1991. [66] David Kofke. Free energy methods in molecular simulation. Fluid Phase Equilibria, 228-229:41–48, 2005. 81 Bibliography [67] Johannes Kästener. Umbrella integration in two or more reaction coor- dinates. The Journal of Chemical Physics, 131:034109–1–8, 2009. [68] Clara D. Christ, Alan E. Mark, and Wilfred F. van Gunsteren. Basic ingredients of free energy calculations: A review. Journal of Computa- tional Chemistry, 31(8):1569–1582, 2009. [69] Daniel Trzesniak, Anna-Pitschna E. Kunz, and Wilfred F. van Gun- steren. A comparison of methods to compute the potential of mean force. ChemPhysChem, 8:162–169, 2007. [70] Johannes Kästener. Umbrella sampling. Wiley Interdisciplinary Reviews: Computational Molecular Science, 1(6):932–942, 2011. [71] Bjørn O. Brandsdal, Fredrik Österberg, Martin Almlöf, and Isabella Feierberg. Free energy calculations and ligand binding. Advances in Protein Chemistry, 66:123–158, 2003. [72] Sanghyun Park and Klaus Schulten. Calculating potentials of mean force from steered molecular dynamics simulations. Journal of Chemical Physics, 120(13):5946–5961, 2004. [73] John G. Kirkwood. Statistical mechanics of fluid mixtures. Journal of ch, 3:300–313, The Journal of Chemical Physics. [74] Jayachandran N. Kizhakkedathu. Unpublished data, 2010. [75] Rajesh Kumar Kainthan, Muthiah Gnanamani, Munia Ganguli, Tanay Gosh, Donald E. Brooks, Souvik Maiti, and Jayachandran N. Kizhakke- dathu. Blood compatibility of novel water soluble hyperbranched polyglycerol-based multivalent cationic polymers and their interaction with DNA. Biomaterials, pages 5377–5390, 2006. [76] Johan Åqvist, Carmen Medina, and Jan-Erik Samuelson. A new method for predicting binding affinity in computer-aided drug design. Protein Engineering, 7(3):385–391, 1994. 82 Bibliography [77] ChemAxon. Marvin 5.5.5. www.chemaxon.com, 2011. [78] Jmol. Jmol: an open-source java viewer for chemical structures in 3D. www.jmol.org, 2010. [79] Xiao-Yan Chen, Claire Marchal, Yaroslav Filinchuk, Daniel Imberta, and Marinella Mazzanti. A flexible tripodal ligand linking octametallic ter- bium rings into luminescent polymeric chainsw. ChemComm, (29):3378– 3380, 2008. [80] DrugBank. Fondaparinux sodium. www.drugbank.ca/drugs/DB00569, 2010. [81] Ulf Börjesson and Philippe H. Hünenberger. Explicit-solvent molecular dynamics simulation at constant ph: Methodology and application to small amines. Journal of Chemical Physics, 114(22):9706–9719, 2001. [82] Milan Remko and Claus-Wilhelm von der Lieth. Conformational struc- ture of some trimeric and pentameric structural units of heparin 2007. Journal of Physical Chemistry A, 111:13484–13491, 2007. [83] Milan Remko, Marcel Swart, and F. Matthias Bickelhaupt. Conforma- tional behavior of basic monomeric building units of glycosaminoglycans: Isolated systems and solvent effect. Journal of Physical Chemistry B, 111:2313–2321, 2007. [84] Shiang-Lin Lin, Prabal K. Maiti, and William A. Goddard III. Dynamics and thermodynamics of water in PAMAM dendrimers at subnanosecond time scales. Journal of Physical Chemistry B, 109:8663–8672, 2005. [85] Prabal K. Maiti, Tahir Caǧin, Shiang-Tai Lin, and William A. God- dard III. Effect of solvent and ph on the structure of PAMAM den- drimers. Macromolecules, (38):979–991, 2005. [86] Yong Cui. Parallel stacking of caffeine with riboflavin in aqueous so- lutions: The potential mechanism for hydrotropic solubilization of ri- boflavin. International Journal of Pharmaceutics, 2010. 83 Bibliography [87] Turgut Baştuğ and Serdar Kuyucak. Free energy simulations of single and double ion occupancy in gramidicin a. The Journal of Chemical Physics, 126:105103–1–12, 2007. [88] Toby W. Allen, Turgut Baştuğ, Serdar Kuyucak, and Shin-Ho Chung. Gramicidin a channel as a test ground for molecular dynamics force fields. Biophysical Journal, 84:2159–2168, 2003. [89] M.A. Kostiainen, J.H. Hardy, and D.K. Smith. High-affinity multiva- lent DNA binding by using low-molecular-weight dendrons. Angewandte Chemie International Edition, 44:2556–2559, 2005. [90] Jayachandran N. Kizhakkedathu. Unpublished data, 2011. [91] Jayachandran N. Kizhakkedathu and Louise Creagh. Unpublished data, 2011. [92] C.K. Nisha, Sunkara V. Manorama, Munia Ganguli, Souvik Maiti, and Jayachandran N. Kizhakkedathu. Complexes of poly(ethylene glycol)- based cationic random copolymer and calf thymus DNA: A complete biophysical characterization. Langmuir, 20:2386–2396, 2004. [93] Coulter Co. Instrument Co. COULTER® DELSA 440, pn 4235648a edition. [94] Biacore. www.biacore.com/lifesciences/technology/introduction/data in- teraction/index.html, 2009. [95] Alexander V. Kabanov. Taking polycation gene delivery systems from in vitro to in vivo. Pharmaceutical Science & Technology Today, 2:365–372, 1999. [96] Kazunori Kataoka, Atsushi Harada, and Yukio Nagasaki. Block copoly- mer micelles for drugdelivery: Design, characterization and biological significance. Advanced Drug Delivery Reviews, 47(1):113–131, 2001. 84 [97] Yoshinori Kakizawa and Kazunori Kataoka. Block copolymer micelles for delivery of gene and related compounds. Advanced Drug Delivery Reviews, 54(2):203–222, 2002. [98] Alexander V. Kabanov, Pierre Lemieuxb, Sergey Vinogradova, and Valery Alakhovb. Pluronic® block copolymers: Novel functional molecules for gene therapy. Advanced Drug Delivery Reviews, 54(2):223– 233, 2002. [99] Crc handbook of chemistry and physics. www.hbcpnetbase.com, 2009. 90th ed. [100] Ira N. Levine. Physical Chemistry. McGraw-Hill, 1221 Avenue of the Americas, New York, NY 10020, 6th edition, 2009. [101] David L. Mobley and Ken A. Dill. Binding of small-molecule ligands to proteins: what you see is not always what you get. Structure, 17:489–498, 2009. [102] Huan-Xiang Zhou and Michael K. Gilson. Theory of free energy and entropy in noncovalent binding. Chemical Reviews, 109(9):4092–4107, 2009. 85 Appendix A MATLAB Code for the Electrostatic Interactions Model Program A.1 This program calculates ka using the Equations 2.2-2.8 for HBSPCMs over a range of rmin, the minimal distance between the R4-1 binding units. For each value of rmin, ka is calculated by calling the function findka rigid which calculates ka for HBSPCMs containing a varying number of rigid binding units attached to the molecule. The program executes a set number of iterations and from the results for each HBSPCM tested calculates the geometric average and standard deviation for ka and the arithmetic average and standard deviation for the number of fondaparinux molecules bound (f). Average values of ka and f are needed because the locations of the R4-1 binding units on the surface of HBSPCM are randomly generated for each iteration. Also, the geometric average and not the arithmetic average is needed for ka because the values span a range of orders of magnitude. Functions: findka rigid clear all format short e global rH Lfond g N rR4 iterations = 1000; % Number of iterations used to calculate ka % for each polymer m = 20; % Maximum number of binding units on HBSPCM a = 6e−10; % Minimal distance of approach [m]. Value % from [44] rH = 4.0e−9; % Radius of HBSPCM [m]. Value from [18]. rR4 = 0.704e−9; % Radius of binding unit R4−1 [m]. Calculated % with data from [99] 86 Appendix A. MATLAB Code for the Electrostatic Interactions Model Lfond = 2.522e−9; % Length of fondaparinux assuming 2D % representation of fondaparinux [m]. Value % from [78]. num rmin = 25; % Number of minimal distance between binding % units, R4−1, to be tested dcut = linspace(0,0.5*Lfond+rR4,num rmin); % Minimal distance % between binding % units, R4−1 [m] length dcut = length(dcut); g = 3; % Number of charges grouped for fondaparinux % binding N = [3:m]; % Range of number of R4−1 binding units that % will be used to generate HBSPCMs for i = 1:length dcut for j=1:iterations count iterations = j; [ka iter(j,:),f iter(j,:)] = findka rigid(m,a,dcut(i)); end ka avg = geomean(ka iter); std ka = geostd(ka iter,0); ka dcut(i,:) = ka avg(1,:); std ka dcut(i,:) = std ka(1,:); f avg = mean(f iter); std f = std(f iter); f dcut(i,:) = f avg(1,:); std f dcut(i,:) = std f(1,:); end save rigid ka vary rmin.out ka dcut −ASCII save rigid stdka vary rmin.out std ka dcut −ASCII save rigid f vary rmin.out f dcut −ASCII save rigid stdf vary rmin.out std f dcut −ASCII 87 Appendix A. MATLAB Code for the Electrostatic Interactions Model findka rigid The function findka rigid calculates ka for HBSPCMs with a varying number of rigid binding units attached to their surface. In order to do so the function calls the functions: fonda energy to calculate the internal energy of the fon- daparinux molecule, gen charges and poly energy to construct an HBSPCM and calculate its internal energy respectively, bind rigid charges to determine fondaparinux binding to rigid binding units, and complex energy to calculate the internal energy of the bound complex. In the case that HBSPCM contains only three binding units and fondaparinux binds to it, the vectors containing the x,y,z coordinates of the binding units (i.e. Xm, Ym, and Zm) contain only one element and therefore Ucomplex = 0. Input arguments: m = maximum number of binding units on HBSPCM a = minimal distance of approach [m] dcut = minimal distance between binding units, R4-1 [m] Output arguments: ka = m x i matrix of values for ka, where i is the number of iterations set by the program calling this function f = number of fondaparinux molecules bound to one HBSPCM Functions: fonda energy gen charges poly energy bind rigid charges complex energy function [ka,f] = findka rigid(m,a,dcut) global rH K kao eo er Lfond g N Wfond = 1.044e−9; % Width of fondaparinux assuming 2D % representation of fondaparinux[m]. % Value from [78]. rF = 1.268e−9; % Diagonal of fondaparinux assuming % 2D representation of fondaparinux [m]. % Value from [78]. cH = 3; % Elementary charge of each binding % unit, R4−1 [atomic units]. Value % from [77]. qH = cH*1.602176487e−19; % Charge of each binding % unit [C] 88 Appendix A. MATLAB Code for the Electrostatic Interactions Model cF = −1; % Elementary charge of one % fondaparinux molecule % [atomic units]. Value from % [77]. qF = cF*1.602176487e−19; % Charge of each fondaparinux % molecule [C] cC = (g*cH)+(cF*10); % Elementary charge of each % (3x(R4−1))−fondaparinux % complex qC = cC*1.602176487e−19; % Charge of each (3x(R4−1))− % fondaparinux complex [C] % DEBYE−HUCKEL PARAMETER, K [1/m]: F = 96485.34; % Faraday's constant [C/mol]. Value from [99]. I = 150; % Ionic strength of solution (physiological) % [mM = mol/mˆ3] eo = 8.854e−12; % Vacuum permittivity [Cˆ2/(J*m)]. Value % from [99]. er = 80; % Dielectric constant at T for water % [unitless]. Value from [99]. R = 8.314472; % Gas constant [J/(K*mol)]. Value from [99]. T = 25 + 273.15; % Temperature of solution (experimental) [K] K = sqrt((2*(Fˆ2)*I)/(eo*er*R*T)); % Debye−Huckel parameter [1/m] % BASAL RATE CONSTANT, kao [1/(M*s)]: visc = 0.90e−3; % Viscosity of the solvent (i.e. water at % 25oC) [Pa*s]. % Value from [99]. kB = 1.3806504e−23; % Boltzmann constant [J/K]. Value from [99]. NA = 6.02214179e23; % Avogadro constant [1/mol]. Value from [99]. DH = (kB*T)/(6*pi()*visc*rH); % Diffusion constant of HBSPCM DF = (kB*T)/(6*pi*visc*rF); % Diffusion constant of % fondaparinux D = DH + DF; % Diffusion constant of % complex Rint = Lfond/(log((2*Lfond)/Wfond));% Interaction radius [m] kao = 4*pi*NA*D*Rint; % Basal rate constant % [1/(M*s)] 89 Appendix A. MATLAB Code for the Electrostatic Interactions Model % DEBYE−HUCKEL ENERGIES Ufond = fonda energy(a,qF); n = length(N); ka = zeros(1,n); for k = 1:n numbercharges = N(k); X = []; Y = []; Z = []; [X,Y,Z] = gen charges(N(k),dcut); Upoly(k) = poly energy(X,Y,Z,N(k),a,qH); Xm = []; Ym = []; Zm = []; chargem = []; [Xm,Ym,Zm,chargem,f temp] = ... bind rigid charges(X,Y,Z,N(k),qH,qC); f(k) = f temp; Ucomplex(k) = complex energy(a,Xm,Ym,Zm,chargem); if f(k) ˜= 0 DeltaU(k) = Ucomplex(k) − Upoly(k) − f(k)*Ufond; ka(k) = kao*exp(−(DeltaU(k)/(kB*T))*(1/(1+K*a))); else ka(k) = 1; end end 90 Appendix A. MATLAB Code for the Electrostatic Interactions Model fonda energy The function fonda energy calculates the Debye-Hückel energy of fonda- parinux. Input arguments: a = minimal distance of approach [m] qF = charge of each fondaparinux molecule [C] Output arguments: Ufond = Debye-Hückel energy of fondaparinux [J]. Scalar. function Ufond = fonda energy(a,qF) global K eo er % Location of fondaparinux charges on the x−axis. The x−axis % is the major semi−axis of the ellipsoid (i.e. fondaparinux). % The origin is located so the x−axis is on the first quadrant % only. Coordinates from [78]. Xf = [0.196e−9; 0.388e−9; 0.604e−9; 0.962e−9; 1.380e−9; ... 1.443e−9; 1.803e−9 1.895e−9; 2.337e−9; 2.371e−9]; % Location of fondaparinux charges on the y−axis. The y−axis % is the minor semi−axis of the ellipsoid. The origin is located % so that the y−axis is on the first quadrant only. Coordinates % from [78]. Yf = [0; 1.044e−9; 0; 0; 1.044e−9; 0; 0; 1.044e−9; 1.044e−9; 0]; Ufond = 0; for i = 1:10 for j = i+1:10 rfond = sqrt((Xf(i)−Xf(j))ˆ2 + (Yf(i)−Yf(j))ˆ2); Ufond = Ufond + exp(−K*(rfond−a))/rfond*(qFˆ2); end end Ufond = Ufond/(4*pi*eo*er*(1+(K*a))); 91 Appendix A. MATLAB Code for the Electrostatic Interactions Model gen charges The function gen charges generates random points (i.e. R4-1 charges) on an HBSPCM surface (i.e. sphere) that are separated from each other by a distance (arc length) of at least rmin (i.e. dcut). Input arguments: num = number of charges on HBSPCM. Scalar. dcut = minimum distance between R4-1 binding units on HBSPCM [m] Ouput arguments: X = x-coordinates of charges on HBSPCM - must be vector Y = y-coordinates of charges on HBSPCM - must be vector of same size as X Z = z-coordinates of charges on HBSPCM - must be vector of same size as X function [X,Y,Z] = gen charges(num,dcut) global rH X = []; Y = []; Z = []; X(1) = −rH + 2*rH*rand; Y(1) = −sqrt(rHˆ2−X(1)ˆ2) + 2*sqrt(rHˆ2−X(1)ˆ2)*rand; signt = ((−1))ˆ(round(rand)); numt = sqrt(rHˆ2 − X(1)ˆ2 − Y(1)ˆ2); Z(1) = signt*numt; lx = 1; for i=2:num; countx = 1; rcut = []; while lx < i Xbin = −rH + 2*rH*rand; Ybin = −sqrt(rHˆ2−Xbinˆ2) + 2*sqrt(rHˆ2−Xbinˆ2)*rand; signt = ((−1))ˆ(round(rand)); numt = sqrt(rHˆ2 − Xbinˆ2 − Ybinˆ2); Zbin = signt*numt; 92 Appendix A. MATLAB Code for the Electrostatic Interactions Model for j = 1:lx rcut(j) = sqrt((Xbin−X(j))ˆ2 + (Ybin−Y(j))ˆ2 + ... (Zbin−Z(j))ˆ2); arc cut(j) = rH*2*asin((rcut(j))/(2*rH)); end if min(arc cut) >= dcut X(i) = Xbin; Y(i) = Ybin; Z(i) = Zbin; lx = lx + 1; else countx = countx + 1; end if countx > 1000 error('I cannot do this anymore') end end end 93 Appendix A. MATLAB Code for the Electrostatic Interactions Model poly energy The function poly energy calculates the Debye-Hückel energy of HBSPCM. Input arguments: X = x-coordinates of charges on HBSPCM - must be vector Y = y-coordinates of charges on HBSPCM - must be vector of same size as X Z = z-coordinates of charges on HBSPCM - must be vector of same size as X num = number of charges on HBSPCM. Scalar. a = minimal distance of approach [m] qH = charge of each R group [C] Output arguments: Upoly = Debye-Hückel energy of HBSPCM [J]. Vector. function Upoly = poly energy(X,Y,Z,num,a,qH) global K eo er Upoly = 0; for i = 1:num for j = i+1:num rpoly = sqrt((X(i)−X(j))ˆ2 + (Y(i)−Y(j))ˆ2 + ... (Z(i)−Z(j))ˆ2); Upoly = Upoly + exp(−K*(rpoly−a))/rpoly*(qHˆ2); end end Upoly = Upoly/(4*pi*eo*er*(1+(K*a))); 94 Appendix A. MATLAB Code for the Electrostatic Interactions Model bind rigid charges The function bind rigid charges groups rigid R4-1 binding units in sets of three according to their proximity (arc length) by calling the function best group 3 and then checks if the binding criterion is met. Input arguments: X = x-coordinates of charges on HBSPCM - must be vector Y = y-coordinates of charges on HBSPCM - must be vector of same size as X Z = z-coordinates of charges on HBSPCM - must be vector of same size as X num = number of charges on HBSPCM. Scalar. qH = charge of each R group [C] qC = charge of each (3x(R4-1))-fondaparinux complex [C] Output arguments: Xm = x-coordinates of charges on HBSPCM - must be vector Ym = y-coordinates of charges on HBSPCM - must be vector of same size as Xm Zm = z-coordinates of charges on HBSPCM - must be vector of same size as Xm chargem = charges of Xm,Ym,Zm - must be a vector f = number of fondaparinux molecules bound to one HBSPCM Functions: best group 3 function [Xm,Ym,Zm,chargem,f]=bind rigid charges(X,Y,Z,num,qH,qC) global rH Lfond g % Groups binding units R4−1 in sets of three according to their % proximity (arc length) I = []; I extras = []; Iv = []; Iv extras = []; X group = []; Y group = []; Z group = []; [I, I extras] = best group 3(rH,X,Y,Z); % Converts matrix I and I extras into column vectors Iv = reshape(I,[],1); 95 Appendix A. MATLAB Code for the Electrostatic Interactions Model Iv extras = reshape(I extras,[],1); % Organizes coordinates of charges on sphere according to the % indexes given by I (now a column vector) dim I = size(I); num col I = dim I(2); for i = 1:(g*num col I) X group(i,1) = X(Iv(i,1)); Y group(i,1) = Y(Iv(i,1)); Z group(i,1) = Z(Iv(i,1)); end % Calculates the arc length between the charges within a group % of three R4−1 binding units counter = 0; a counter = 1:g:(length(X group)−(g−1)); b counter = a counter+(g−1); count = 0; while (counter < (length(X group)/g)) counter = counter + 1; for i = a counter(counter):b counter(counter) for j = i+1:b counter(counter) count = count + 1; r(count) = sqrt((X group(i)−X group(j))ˆ2 + ... (Y group(i)−Y group(j))ˆ2 + ... (Z group(i) − Z group(j))ˆ2); d(count) = rH*2*asin((r(count))/(2*rH)); end end end % Checks if the arc lengths between charges in a group of three % R4−1 binding units are less than the length of fondaparinux. % If they are, the coordinates of the first point are transferred % to a new matrix, 'the complex matrix' and a charge of % qC=gxqH+10xqF is assigned to this new point on the % HBSPCM−fondaparinux complex. Otherwise, their coordinates are % transferred to 'the complex matrix' with their charges of qH. count = 0; 96 Appendix A. MATLAB Code for the Electrostatic Interactions Model f = 0; Xm = []; Ym = []; Zm = []; chargem = []; for i = 1:num col I count = count + 1; if (((d(1+g*(i−1)) + d(2+g*(i−1)))<=Lfond)... | | ((d(1+g*(i−1)) + d(3+g*(i−1)))<=Lfond)... | | ((d(2+g*(i−1)) + d(3+g*(i−1)))<=Lfond)) Xm(count) = X group(1+g*(i−1)); Ym(count) = Y group(1+g*(i−1)); Zm(count) = Z group(1+g*(i−1)); chargem(count) = qC; f = f + 1; else Xm(count) = X group(1+g*(i−1)); Xm(count+1) = X group(2+g*(i−1)); Xm(count+2) = X group(3+g*(i−1)); Ym(count) = Y group(1+g*(i−1)); Ym(count+1) = Y group(2+g*(i−1)); Ym(count+2) = Y group(3+g*(i−1)); Zm(count) = Z group(1+g*(i−1)); Zm(count+1) = Z group(2+g*(i−1)); Zm(count+2) = Z group(3+g*(i−1)); chargem(count) = qH; chargem(count+1) = qH; chargem(count+2) = qH; count = count + 2; end end % Adds to 'the complex matrix' those R4 points that did % not make it into the groups of three with the function 97 Appendix A. MATLAB Code for the Electrostatic Interactions Model % best group 3 above. Adds to the chargem vector the % corresponding charges for these points. if num/g ˜= 1 for i = 1:(length(Iv extras)) Xm(count+i) = X(Iv extras(i,1)); Ym(count+i) = Y(Iv extras(i,1)); Zm(count+i) = Z(Iv extras(i,1)); chargem(count+i) = qH; end end 98 Appendix A. MATLAB Code for the Electrostatic Interactions Model best group 3 The function best group 3 calculates the arc length perimeter between all pos- sible groups of three points on the sphere (i.e. charges on HBSPCM) and groups three binding units based on their proximity by finding the lowest perimeter possible. Input arguments: rH = radius of HBSPCM X = x-coordinates of charges on HBSPCM - must be a vector. Y = y-coordinates of charges on HBSPCM - must be a vector of same size as X. Z = z-coordinates of charges on HBSPCM - must be a vector of same size as X. Output arguments: I = matrix with index location of charges in X, Y, and Z extras = index location of any extra charges that can’t be grouped into a group of three R4-1 binding units function [I, extras] = best group 3(rH,X,Y,Z) N = length(X); I = []; tri = []; bin = []; d = []; bin = 1:N; d = inf(N−1,N−1); tri = inf(N−2,N−2,N−2); for f = 1:N−1 for e = f:N−1 r(e,f) = sqrt((X(f)−X(e+1))ˆ2 + (Y(f)−Y(e+1))ˆ2 + ... (Z(f) − Z(e+1))ˆ2); % Produces an array that contains arc distances % between all possible two−point groups d(e,f) = rH*2*asin((r(e,f))/(2*rH)); end end 99 Appendix A. MATLAB Code for the Electrostatic Interactions Model % Creates 3D array with the perimeters of all possible groups % of three points using the distances calculated in the d % matrix above for j = 1:N−2 for k = 1:N−(j+1) for i = k+(j−1):N−2 tri(i,j,k) = d(j+k−1,j)+d(i+1,j)+d(i+1,j+k); end end end counter = 0; % Determines the lowest possible value for the a perimeter % of three binding units. It will continue as long as there % are enough points to make a group of three binding units while (floor(N/3) − counter) >= 1 min value = min(min(min(tri))); pos = find(tri==min value); pos = pos(1,1); [loc1,emp1,emp2] = find(tri==min value); loc1 = loc1(1,1); % Ensures that if there are more than % two instances of the minimum the program % just takes the first one loc3 = floor(pos/((N−2)ˆ2+1)) + 1; loc2 = floor((pos−(loc3−1)*(N−2)ˆ2)/(N−1))+1; loc = [loc1,loc2,loc3]; % Coordinates in tri matrix of minimum % perimeter for a group of three % binding units % loc variables are used to turn the position that the % function 'find' returns into coordinates in the array I = [I [loc2;loc3+loc2;loc1+2;]]; bin(loc2)=0; bin(loc3+loc2)=0; % bin is used to track which points % are not grouped bin(loc1+2)=0; 100 Appendix A. MATLAB Code for the Electrostatic Interactions Model % for loop below rids the tri array of any group of three % binding units that contain one or more of the points % allocated to a group of three binding units above. It % does this by looking at every individual value of the % tri array, turns its coordinates into the indices of % the actual points and compares them to the indices of % the points allocated to I above. for r = 1:N−2 for s = 1:N−(r+1) for q = s+(r−1):N−2 pt1 = r; pt2 = r+s; pt3 = q+2; if pt1==loc2 | | pt1==loc3+loc2 | | pt1==loc1+2 tri(q,r,s) = inf; elseif pt2==loc2 | | pt2==loc3+loc2 | | pt2==loc1+2 tri(q,r,s) = inf; elseif pt3==loc2 | | pt3==loc3+loc2 | | pt3==loc1+2 tri(q,r,s) = inf; end end end end counter = counter + 1; end % Outputs extras and I extras = find(bin˜=0); I; 101 Appendix A. MATLAB Code for the Electrostatic Interactions Model complex energy The function complex energy calculates the Debye-Hückel energy of the HBSPCM-fondaparinux complex. Input arguments: a = minimal distance of approach [m] Xm = x-coordinates of charges on HBSPCM - must be vector Ym = y-coordinates of charges on HBSPCM - must be vector of same size as Xm Zm = z-coordinates of charges on HBSPCM - must be vector of same size as Xm chargem = charges of Xm,Ym,Zm - must be a vector Output arguments: Ucomplex = Debye-Hückel energy of HBSPCM-fondaparinux complex [J]. Vec- tor. function Ucomplex = complex energy(a,Xm,Ym,Zm,chargem) global K eo er Ucomplex = 0; for i = 1:length(Xm) for j = i+1:length(Xm) rcomp = sqrt((Xm(i)−Xm(j))ˆ2 + (Ym(i)−Ym(j))ˆ2 + ... (Zm(i)−Zm(j))ˆ2); Ucomplex = Ucomplex + ... (chargem(i)*chargem(j)*exp(−K*(rcomp−a))/rcomp); end end Ucomplex = Ucomplex /(4*pi*eo*er*(1+(K*a))); 102 Appendix A. MATLAB Code for the Electrostatic Interactions Model Program A.2 This program calculates ka using the Equations 2.2-2.8 for HBSPCMs over a range of rmin, the minimal distance between the R4-1 binding units. For each value of rmin, ka is calculated by calling the function findka flexible which cal- culates ka for HBSPCMs containing a varying number of flexible binding units attached to the molecule. The program executes a set number of iterations and from the results for each HBSPCM tested calculates the geometric aver- age and standard deviation for ka and the arithmetic average and standard deviation for the number of fondaparinux molecules bound (f). Average val- ues of ka and f are needed because the locations of the R4-1 binding units on the surface of HBSPCM are randomly generated for each iteration. Also, the geometric average and not the arithmetic average is needed for ka because the values span a range of orders of magnitude. Functions: findka flexible clear all format short e global rH Lfond g N rR4 iterations = 1000; % Number of iterations used to calculate ka % for each polymer m = 20; % Maximum number of binding units on HBSPCM a = 6e−10; % Minimal distance of approach [m]. Value % from [44] rH = 4.0e−9; % Radius of HBSPCM [m]. Value from [18]. rR4 = 0.704e−9; % Radius of binding unit R4−1 [m]. Calculated % with data from [99] Lfond = 2.522e−9; % Length of fondaparinux assuming 2D % representation of fondaparinux [m]. Value % from [78]. num rmin = 25; % Number of minimal distance between binding % units, R4−1, to be tested dcut = linspace(0,0.5*Lfond+rR4,num rmin); % Minimal distance % between binding % units, R4−1 [m] length dcut = length(dcut); g = 3; % Number of charges grouped for fondaparinux 103 Appendix A. MATLAB Code for the Electrostatic Interactions Model % binding N = [3:m]; % Range of number of R4−1 binding units % that will be used to generate HBSPCMs for i = 1:length dcut for j=1:iterations count iterations = j; [ka iter(j,:),f iter(j,:)] = findka rigid(m,a,dcut(i)); end ka avg = geomean(ka iter); std ka = geostd(ka iter,0); ka dcut(i,:) = ka avg(1,:); std ka dcut(i,:) = std ka(1,:); f avg = mean(f iter); std f = std(f iter); f dcut(i,:) = f avg(1,:); std f dcut(i,:) = std f(1,:); end save flex ka vary rmin.out ka dcut −ASCII save flex stdka vary rmin.out std ka dcut −ASCII save flex f vary rmin.out f dcut −ASCII save flex stdf vary rmin.out std f dcut −ASCII 104 Appendix A. MATLAB Code for the Electrostatic Interactions Model findka flexible The function findka flexible calculates ka for HBSPCMs with a varying number of flexible binding units attached to their surface. In order to do so the function calls the functions: fonda energy, gen charges, poly energy, complex energy and bind flex charges to determine fondaparinux binding to flexible binding units. In the case that HBSPCM contains only three binding units and fonda- parinux binds to it, the vectors containing the x,y,z coordinates of the binding units (i.e. Xm, Ym, and Zm) contain only one element and therefore Ucom- plex = 0. Input arguments: m = maximum number of binding units on HBSPCM a = minimal distance of approach [m] dcut = minimal distance between binding units, R4-1 [m] Output arguments: ka = m x i matrix of values for ka, where i is the number of iterations set by the program calling this function f = number of fondaparinux molecules bound to one HBSPCM Functions: fonda energy gen charges poly energy bind flex charges complex energy function [ka,f] = findka flexible(m,a,dcut) global rH K kao eo er Lfond g N Wfond = 1.044e−9; % Width of fondaparinux assuming 2D % representation of fondaparinux[m]. % Value from [78]. rF = 1.268e−9; % Diagonal of fondaparinux assuming % 2D representation of fondaparinux [m]. % Value from [78]. cH = 3; % Elementary charge of each binding % unit, R4−1 [atomic units]. Value % from [77]. qH = cH*1.602176487e−19; % Charge of each binding % unit [C] cF = −1; % Elementary charge of one % fondaparinux molecule 105 Appendix A. MATLAB Code for the Electrostatic Interactions Model % [atomic units]. Value from % [77]. qF = cF*1.602176487e−19; % Charge of each fondaparinux % molecule [C] cC = (g*cH)+(cF*10); % Elementary charge of each % (3x(R4−1))−fondaparinux % complex qC = cC*1.602176487e−19; % Charge of each (3x(R4−1))− % fondaparinux complex [C] % DEBYE−HUCKEL PARAMETER, K [1/m]: F = 96485.34; % Faraday's constant [C/mol]. Value % from [99]. I = 150; % Ionic strength of solution (physiological) % [mM = mol/mˆ3] eo = 8.854e−12; % Vacuum permittivity [Cˆ2/(J*m)]. % Value from [99]. er = 80; % Dielectric constant at T for water % [unitless]. Value from [99]. R = 8.314472; % Gas constant [J/(K*mol)]. Value from [99]. T = 25 + 273.15; % Temperature of solution (experimental) [K] K = sqrt((2*(Fˆ2)*I)/(eo*er*R*T)); % Debye−Huckel parameter [1/m] % BASAL RATE CONSTANT, kao [1/(M*s)]: visc = 0.90e−3; % Viscosity of the solvent (i.e. water % at 25oC) [Pa*s]. % Value from [99]. kB = 1.3806504e−23; % Boltzmann constant [J/K]. Value from [99]. NA = 6.02214179e23; % Avogadro constant [1/mol]. Value from [99]. DH = (kB*T)/(6*pi()*visc*rH); % Diffusion constant of HBSPCM DF = (kB*T)/(6*pi*visc*rF); % Diffusion constant of % fondaparinux D = DH + DF; % Diffusion constant of % complex Rint = Lfond/(log((2*Lfond)/Wfond));% Interaction radius [m] kao = 4*pi*NA*D*Rint; % Basal rate constant % [1/(M*s)] % DEBYE−HUCKEL ENERGIES 106 Appendix A. MATLAB Code for the Electrostatic Interactions Model Ufond = fonda energy(a,qF); n = length(N); ka = zeros(1,n); for k = 1:n numbercharges = N(k); X = []; Y = []; Z = []; [X,Y,Z] = gen charges(N(k),dcut); Upoly(k) = poly energy(X,Y,Z,N(k),a,qH); Xm = []; Ym = []; Zm = []; chargem = []; [Xm,Ym,Zm,chargem,f temp] = ... bind flex charges(X,Y,Z,N(k),qH,qC); f(k) = f temp; Ucomplex(k) = complex energy(a,Xm,Ym,Zm,chargem); if f(k) ˜= 0 DeltaU(k) = Ucomplex(k) − Upoly(k) − f(k)*Ufond; ka(k) = kao*exp(−(DeltaU(k)/(kB*T))*(1/(1+K*a))); else ka(k) = 1; end end 107 Appendix A. MATLAB Code for the Electrostatic Interactions Model bind flex charges The function bind flex charges groups flexible R4-1 binding units in sets of three according to their proximity (arc length) by calling the function best group 3 and then checks if the binding criterion is met. Input arguments: X = x-coordinates of charges on HBSPCM - must be vector Y = y-coordinates of charges on HBSPCM - must be vector of same size as X Z = z-coordinates of charges on HBSPCM - must be vector of same size as X num = number of charges on HBSPCM. Scalar. qH = charge of each R group [C] qC = charge of each (3x(R4-1))-fondaparinux complex [C] Output arguments: Xm = x-coordinates of charges on HBSPCM - must be vector Ym = y-coordinates of charges on HBSPCM - must be vector of same size as Xm Zm = z-coordinates of charges on HBSPCM - must be vector of same size as Xm chargem = charges of Xm,Ym,Zm - must be a vector f = number of fondaparinux molecules bound to one HBSPCM Functions: best group 3 function [Xm,Ym,Zm,chargem,f] = bind flex charges(X,Y,Z,num,qH,qC) global rH Lfond g % Groups binding units R4−1 in sets of three according to their % proximity (arc length) I = []; I extras = []; Iv = []; Iv extras = []; X group = []; Y group = []; Z group = []; [I, I extras] = best group 3(rH,X,Y,Z); % Converts matrix I and I extras into column vectors Iv = reshape(I,[],1); 108 Appendix A. MATLAB Code for the Electrostatic Interactions Model Iv extras = reshape(I extras,[],1); % Organizes coordinates of charges on sphere according to the % indexes given by I (now a column vector) dim I = size(I); num col I = dim I(2); for i = 1:(g*num col I) X group(i,1) = X(Iv(i,1)); Y group(i,1) = Y(Iv(i,1)); Z group(i,1) = Z(Iv(i,1)); end % Calculates the arc length between the charges within a group % of three R4−1 binding units counter = 0; a counter = 1:g:(length(X group)−(g−1)); b counter = a counter+(g−1); count = 0; while (counter < (length(X group)/g)) counter = counter + 1; for i = a counter(counter):b counter(counter) for j = i+1:b counter(counter) count = count + 1; r(count) = sqrt((X group(i)−X group(j))ˆ2 + ... (Y group(i)−Y group(j))ˆ2 + ... (Z group(i) − Z group(j))ˆ2); d(count) = rH*2*asin((r(count))/(2*rH)); end end end % Checks if the arc lengths between charges in a group of three % R4−1 binding units are less than the length of fondaparinux. % If they are, the coordinates of the first point are transferred % to a new matrix, 'the complex matrix' and a charge of % qC=gxqH+10xqF is assigned to this new point on the % HBSPCM−fondaparinux complex. Otherwise, their coordinates are % transferred to 'the complex matrix' with their charges of qH. count = 0; 109 Appendix A. MATLAB Code for the Electrostatic Interactions Model f = 0; Xm = []; Ym = []; Zm = []; chargem = []; for i = 1:num col I count = count + 1; % Adds the arc lengths: AB+AC, AB+BC, and AC+BC len1 = d(1+g*(i−1)) + d(2+g*(i−1)); len2 = d(1+g*(i−1)) + d(3+g*(i−1)); len3 = d(2+g*(i−1)) + d(3+g*(i−1)); M = [len1,len2,len3]; min len = min(M); if min len <= (Lfond + rR4) Xm(count) = X group(1+g*(i−1)); Ym(count) = Y group(1+g*(i−1)); Zm(count) = Z group(1+g*(i−1)); chargem(count) = qC; f = f + 1; else Xm(count) = X group(1+g*(i−1)); Xm(count+1) = X group(2+g*(i−1)); Xm(count+2) = X group(3+g*(i−1)); Ym(count) = Y group(1+g*(i−1)); Ym(count+1) = Y group(2+g*(i−1)); Ym(count+2) = Y group(3+g*(i−1)); Zm(count) = Z group(1+g*(i−1)); Zm(count+1) = Z group(2+g*(i−1)); Zm(count+2) = Z group(3+g*(i−1)); chargem(count) = qH; chargem(count+1) = qH; 110 Appendix A. MATLAB Code for the Electrostatic Interactions Model chargem(count+2) = qH; count = count + 2; end end % Adds to 'the complex matrix' those R4 points that did % not make it into the groups of three with the function % best group 3 above. Adds to the chargem vector the % corresponding charges for these points. if num/g ˜= 1 for i = 1:(length(Iv extras)) Xm(count+i) = X(Iv extras(i,1)); Ym(count+i) = Y(Iv extras(i,1)); Zm(count+i) = Z(Iv extras(i,1)); chargem(count+i) = qH; end end 111 Appendix A. MATLAB Code for the Electrostatic Interactions Model Program A.3 This program calculates ka using the Equations 2.2-2.8, as many times as the value set for the variable iterations. In order to calculate ka, it calls the function findka flexible which calculates ka for HBSPCMs that differ on the number of flexible binding sites. It then calculates the geometric average and standard deviation for ka and the arithmetic average and standard deviation for binding (f) for each HBSPCM tested. An average value of ka and f is needed because the location of the R4-1 binding units on the surface of HBSPCM is randomly generated with each iteration. Also, the geometric average and not the arithmetic average is needed because the ka values span over a wide range. Functions: findka flexible clear all format short e global rH Lfond g N rR4 iterations = 1000; % Number of iterations used to calculate ka % for each polymer m = 20; % Maximum number of binding units on HBSPCM a = 6e−10; % Minimal distance of approach [m]. Value % from [44] rH = 4e−9; % Radius of HBSPCM [18] [m] rR4 = 0.704e−9; % Radius of binding unit R4−1. Calculated % with data from [99] Lfond = 2.522e−9; % Length of fondaparinux assuming 2D % representation of fondaparinux [m]. Value % from [78]. dcut = 4.533e−10; % Minimal distance between binding units, % R4−1 [m]. Found using Program A.1 and % experimental data g = 3; % Number of charges grouped for fondaparinux % binding N = [3:m]; % Range of number of R4−1 binding units % that will be used to generate HBSPCMs count = 0; 112 Appendix A. MATLAB Code for the Electrostatic Interactions Model for i=1:iterations [ka(i,:),f(i,:)] = findka flexible(m,a,dcut); count = count + 1; end ka sort = sort(ka,1,'descend'); f sort = sort(f,1,'descend'); ka avg = geomean(ka); std ka = geostd(ka,0); f avg = mean(f); std f = std(f); save all kas flex.out ka sort −ASCII save all fs flex.out f sort −ASCII save kas flex.out ka avg −ASCII save kas flex std.out std ka −ASCII save f flex.out f avg −ASCII save f flex std.out std f −ASCII 113 Appendix B Modification of Force Field Types and Charges for Fondaparinux As with all current force fields, COMPASS does not properly describe hep- arin [82, 83], in particular some of the sulfonamide functional groups found within heparin derivatives such as fondaparinux [41]. Therefore, in the MD simulations performed in this thesis, some modifications were made to the force field types and charges assigned by COMPASS to the three sulfonamide groups in fondaparinux. This appendix describes the details pertaining to these changes. Fondaparinux has three sulfonamides in which the sulfur atom is connected to an alcohol with the hydrogen removed. Although COMPASS does not have bond parameters for this dissociated sulfonamide it does for the undissociated analog. Therefore, the undissociated analog of the sulfonamide oxy anion was built with the Visualizer module of MS and its force field types were properly assigned using COMPASS. These force field types were then manually assigned to the three dissociated sulfonamide oxy anions in fondaparinux (Figure B.1). 114 Appendix B. Modification of Force Field Types and Charges for Fondaparinux (a) (b) Figure B.1: Force field types for: (a) the undissociated analog and (b) the dissociated analog of the sulfonamide oxy anion in fondaparinux. Force field types were assigned using the force field COMPASS for (a) and then manually assigned to (b). Fondaparinux also contains sulfonyl methoxy oxy groups. Since COMPASS properly assigns the sulfonyl methoxy oxy charges (Figure B.2), these charges were assigned to the three sulfonamide oxy anions as a first approximation. These charges were then adjusted based on a desired net charge of -10 for fondaparinux with the information in Table B.1. The three sulfonamide frag- ments were selected and their total charge was set to -4.625 e. This modified the partial charges of each atom in the fragment and, in turn, the total charge of the deprotonated fondaparinux molecule to the desired -10 (Figure B.3). (a) (b) Figure B.2: Partial charges of the atoms in: (a) the sulfonyl methoxy oxy an- ion and (b) the sulfonamide group (or fragment) in fondaparinux. Charges were assigned using the force field COMPASS for (a) and then manually assigned to (b). 115 Appendix B. Modification of Force Field Types and Charges for Fondaparinux Table B.1: Fondaparinux charge information required to adjust the molecule’s net charge to the appropriate value of -10. Groups Charge (e) Current total charge of fondaparinux -9.440 Current total charge of sulfonamide fragments -4.065 Current total charge per sulfonamide fragment -1.355 Required total charge of fondaparinux -10.000 Required total charge per sulfonamide fragment -1.542 Required total charge of sulfonamide fragments -4.625 Figure B.3: Final partial charges of the atoms in the sulfonamide groups in fondaparinux, adjusted for a fondaparinux net charge of -10. 116 Appendix C ∆A and ∆G for Systems in the Condensed Phase In Chapter 4, the change in free energy upon binding of a cationic binding unit to a fondaparinux molecule was investigated. The method described in Section 2.3.1 to obtain free-energy profiles (PMF) was derived using the canon- ical ensemble (NVT) [69], for which the Helmholtz free energy (A) is minimized at equilibrium [100]. However, the Gibbs free energy (G) is commonly reported in the literature [101] since constant pressure is a typical experimental condi- tion and ∆G is directly related to the equilibrium association constant (Ka) [102]. Thus in this thesis it was beneficial to be able to calculate ∆G. In this appendix the equivalency of ∆A and ∆G for systems in the condensed phase is verified. For a closed system—a system where there is no transfer of matter between system and surroundings—at constant temperature (T ) and constant volume (V ) the Helmholtz free energy is defined as [100] A = U − T S (C.1) where U and S are the internal energy and entropy of the system, respectively. When the system is held at constant temperature and pressure (P ) the Gibbs free energy is given by [100] G = H − T S = U + P V − T S (C.2) From Equation C.1 and Equation C.2 it is clear that the difference be- tween Helmholtz and Gibbs free energies is the P V contribution. For systems 117 Appendix C. ∆A and ∆G for Systems in the Condensed Phase in the condensed phase the change in volume is negligible and therefore the free-energy changes ∆A and ∆G approximate one another [67, 70]. This ap- proximation was verified with one of the systems studied in this thesis. The model system consisted of an R4-1 cationic binding unit, a deproto- nated fondaparinux molecule, 2,600 water molecules, and appropriate amounts of sodium counterions (Na+) to achieve charge neutrality. Three MD simu- lations under NPT conditions were performed in which the change in P V as a function of restraint distance was observed. The temperature was held at 298 K by the Nose-Hoover thermostat and the pressure was controlled at 0.00 Pa by the Berendsen barostat. The time step was 1 fs and the free-energy calculations were performed as described in Section 4.2.2 but using an inter- molecular separation of fondaparinux and the binding unit (R4-1) of 1.0 Å instead of 0.5 Å for computational speed. The MD simulation results in Figure C.1 indicate that the change in volume is indeed negligible for the systems studied in this thesis and thus the free energy changes ∆A and ∆G are numerically very similar. 118 Appendix C. ∆A and ∆G for Systems in the Condensed Phase Figure C.1: Free energy profile (PMF) and pressure (P ) and (V ) changes dur- ing an MD simulation of the interaction between fondaparinux and the binding unit R4-1, (n= 3). Simulations were performed using the NPT ensemble, a step size of 1.0 Å, and a simulation time of 60 ps at each interval for n number of replicates. 119 Appendix D MATLAB Code for the Improved Electrostatic Interactions Model Program D.1 This program calculates ka using the Equations 2.2-2.8 for HBSPCMs over a range of rmin, the minimal distance between the R4-1 binding units. For each value of rmin, ka is calculated by calling the function findka which calculates ka for HBSPCMs containing a varying number of rigid binding units attached to the molecule. The program executes a set number of iterations and from the results for each HBSPCM tested calculates the geometric average and standard deviation for ka and the arithmetic average and standard deviation for the number of fondaparinux molecules bound (f). Average values of ka and f are needed because the locations of the R4-1 binding units on the surface of HBSPCM are randomly generated for each iteration. Also, the geometric average and not the arithmetic average is needed for ka because the values span a range of orders of magnitude. Functions: findka clear all format short e global rH Lfond g N rR4 d iterations = 1000; % Number of iterations used to calculate ka % for each polymer m = 20; % Maximum number of binding units on HBSPCM rH = 4.0e−9; % Radius of HBSPCM [18] [m] a = 8.5e−10; % Minimal distance of approach [m]. Value 120 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model % from Chapter 4. rR4 = 6.613e−10; % Radius of binding unit R4−1 [m]. 3D % representation. Value from Chapter 4. bd = a; % Binding distance between fondaparinux % and binding unit R4−1. Measured with % respect to centroid of the molecules. % Value from Chapter 4. d = sqrt(2*(bd+rR4)ˆ2); % Diagonal of binding site [m]. % 3D representation. Value from MS. Lfond = 19.939e−10; % Length of fondaparinux of fondaparinux [m]. % 3D representation. Value from Chapter 4. num rmin = 25; % Number of minimal distance between binding % units, R4−1, to be tested dcut = linspace(0,1.5e−9,num rmin); % Minimal distance between % binding units, R4−1 [m] length dcut = length(dcut); g = 3; % Number of charges grouped for fondaparinux % binding N = [3:m]; % Range of number of R4−1 binding units that % will be used to generate HBSPCMs for i = 1:length dcut for j=1:iterations count iterations = j; [ka iter(j,:),f iter(j,:)] = findka(m,a,dcut(i)); end ka avg = geomean(ka iter); std ka = geostd(ka iter,0); ka dcut(i,:) = ka avg(1,:); std ka dcut(i,:) = std ka(1,:); f avg = mean(f iter); std f = std(f iter); f dcut(i,:) = f avg(1,:); std f dcut(i,:) = std f(1,:); se f dcut(i,:) = std f dcut(i,:)/sqrt(iterations); bar f dcut(i,:)= 1.96*se f dcut(i,:); end 121 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model save ka vary rmin.out ka dcut −ASCII save stdka vary rmin.out std ka dcut −ASCII save f vary rmin.out f dcut −ASCII save stdf vary rmin.out std f dcut −ASCII 122 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model findka The function findka calculates ka for HBSPCMs with a varying number of R4-1 binding units attached to their surface. In order to do so the function calls the functions: fonda energy to calculate the internal energy of the fondaparinux molecule, gen charges and poly energy to construct an HBSPCM and calculate its internal energy respectively, binding to determine fondaparinux binding to binding units, and complex energy to calculate the internal energy of the bound complex. In the case that HBSPCM contains only three binding units and fondaparinux binds to it, the vectors containing the x,y,z coordinates of the binding units (i.e. Xm, Ym, and Zm) contain only one element and therefore Ucomplex = 0. Input arguments: m = maximum number of binding units on HBSPCM a = minimal distance of approach [m] dcut = minimal distance between binding units, R4-1 [m] Output arguments: ka = m x i matrix of values for ka, where i is the number of iterations set by the program calling this function f = number of fondaparinux molecules bound to one HBSPCM Functions: fonda energy gen charges poly energy binding complex energy function [ka,f] = findka(m,a,dcut) global rH K kao eo er Lfond N R bound temp R bound Wfond = 10.467e−10; % Width of fondaparinux, 3D [m]. Value from % Chapter 4. rF = 11.812e−10; % Radius of fondaparinux, 3D [m]. Average % of distances measured from O19−O53 and % O19−O62. Value from Chapter 4. cH = 1.29; % Elementary charge of each binding % unit, R4−1 [atomic units]. Value % from Chapter 5.3.1. qH = cH*1.602176487e−19; % Charge of each binding % unit [C] 123 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model cF = −1; % Elementary charge of one % fondaparinux molecule % [atomic units]. Value from % [77]. qF = cF*1.602176487e−19; % Charge of each fondaparinux % molecule [C] % DEBYE−HUCKEL PARAMETER, K [1/m]: F = 96485.34; % Faraday's constant [C/mol]. Value from [99]. I = 150; % Ionic strength of solution (physiological) % [mM = mol/mˆ3] eo = 8.854e−12; % Vacuum permittivity [Cˆ2/(J*m)]. Value % from [99]. er = 80; % Dielectric constant at T for water % [unitless]. Value from [99]. R = 8.314472; % Gas constant [J/(K*mol)]. Value from [99]. T = 25 + 273.15; % Temperature of solution (experimental) [K] K = sqrt((2*(Fˆ2)*I)/(eo*er*R*T)); % Debye−Huckel parameter [1/m] % BASAL RATE CONSTANT, kao [1/(M*s)]: visc = 0.90e−3; % Viscosity of the solvent (i.e. water at % 25oC) [Pa*s]. % Value from [99]. kB = 1.3806504e−23; % Boltzmann constant [J/K]. Value from [99]. NA = 6.02214179e23; % Avogadro constant [1/mol]. Value from [99]. DH = (kB*T)/(6*pi()*visc*rH); % Diffusion constant of HBSPCM DF = (kB*T)/(6*pi*visc*rF); % Diffusion constant of % fondaparinux D = DH + DF; % Diffusion constant of % complex Rint = Lfond/(log((2*Lfond)/Wfond));% Interaction radius [m] kao = 4*pi*NA*D*Rint; % Basal rate constant % [1/(M*s)] % DEBYE−HUCKEL ENERGIES Ufond = fond energy(a,qF); n = length(N); ka = zeros(1,n); 124 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model R bound temp = NaN(n,floor(m/3)); for k = 1:n numbercharges = N(k); X = []; Y = []; Z = []; [X,Y,Z] = gen charges(N(k),dcut); Upoly(k) = poly energy(X,Y,Z,N(k),a,qH); Xm = []; Ym = []; Zm = []; chargem = []; [Xm,Ym,Zm,chargem,f temp] = binding(X,Y,Z,N(k),qH,qF); for j = 1:length(R bound) R bound temp(k,j) = R bound(j); end f(k) = f temp; Ucomplex(k) = complex energy(a,Xm,Ym,Zm,chargem); if f(k) ˜= 0 DeltaU(k) = Ucomplex(k) − Upoly(k) − f(k)*Ufond; ka(k) = kao*exp(−(DeltaU(k)/(kB*T))*(1/(1+K*a))); else ka(k) = 1; end end 125 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model binding The function binding checks if binding units, R4-1, attached to the surface of an HBSPCM are within a specified distance (d) from each other. If this is the case, they are considered a binding site and are assigned a charge that reflects binding to fondaparinux. A binding site contains three or more binding units. Input arguments: X = x-coordinates of charges on HBSPCM - must be vector Y = y-coordinates of charges on HBSPCM - must be vector of same size as X Z = z-coordinates of charges on HBSPCM - must be vector of same size as X n = number of charges on HBSPCM. Scalar. qH = charge of each R group [C] qF = charge of each fondaparinux molecule [C] Output arguments: X c = x-coordinates of charges on HBSPCM - must be vector Y c = y-coordinates of charges on HBSPCM - must be vector of same size as Xm Z c = z-coordinates of charges on HBSPCM - must be vector of same size as Xm charge c = charges of X c,Y c,Z c - must be a vector f = number of fondaparinux molecules bound to one HBSPCM function [X c,Y c,Z c,charge c,f] = binding(X,Y,Z,n,qH,qF) global rH d arc length R bound % Calculates the arc length between all the charges on the % surface of an HBSPCM. The diagonal of the nxn matrix is set % to 'inf' since the distance of a point with itself is zero. arc length = zeros(n,n); v = inf*ones(1,n); arc length = arc length − diag(diag(arc length))+diag(v); for i = 1:n for j = i+1:n dist(i,j) = sqrt((X(i)−X(j))ˆ2 + (Y(i)−Y(j))ˆ2 ... +(Z(i)−Z(j))ˆ2); arc length(i,j) = rH*2*asin((dist(i,j))/(2*rH)); arc length(j,i) = arc length(i,j); end end 126 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model % Randomly selects a binding unit, R4−1, on the surface of an % HBSPCM and checks for other R4−1 binding units that are within % a specified distance, d, to form a binding site. It then % creates vectors with the x,y,z coordinates of the R4−1 binding % units that formed a binding site and vectors with the x,y,z % coordinates of the R4−1 binding units that did not form a % binding site. start = randperm(n); Rs = 1; f = 0; nf = 0; X b = []; Y b = []; Z b = []; charge b = []; X nb = []; Y nb = []; Z nb = []; charge nb = []; R used = [1:n]; count = 0; R bound = []; for i = 1:n loc = []; if R used(start(i)) > 0 for j = 1:n if R used(j) > 0 if arc length(start(i),j) <= d R used(start(i)) = 0; R used(j) = 0; loc = [loc j]; Rs = Rs + 1; end end end 127 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model % If a binding site is formed by three or more R4−1 % binding units, it gets assigned the coordinates from % the binding unit used as the reference point to find % the other binding units 'd' or less away from it. % The charge of the binding site will depend on the % number of binding units present. % Rs = number of binding units in a binding site. if Rs >= 3 count = count + 1; R bound(count) = Rs; f = f + 1; X b(f) = X(start(i)); Y b(f) = Y(start(i)); Z b(f) = Z(start(i)); charge b(f) = qF − qH*Rs; else nf = nf + 1; X nb(nf) = X(start(i)); Y nb(nf) = Y(start(i)); Z nb(nf) = Z(start(i)); charge nb(nf) = qH; R used(start(i)) = start(i); R used(loc) = loc; end Rs = 1; end end % Fondaparinux−HBSPCM complex: % Combines the vectors of the x,y,z coordinates of both the % binding units inside the binding sites and those outside % the binding site. X c = [X b X nb]; Y c = [Y b Y nb]; Z c = [Z b Z nb]; charge c = [charge b charge nb]; 128 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model Program D.2 This program calculates ka using the Equations 2.2-2.8 for HBSPCMs by call- ing the function findka which calculates ka for HBSPCMs containing a varying number R4-1 binding units attached to the molecule. The program executes a set number of iterations and from the results calculates the geometric average and standard deviation for ka and the arithmetic average and standard devia- tion for the number of fondaparinux molecules bound (f) for each HBSPCM tested. Average values of ka and f are needed because the locations of the R4-1 binding units on the surface of HBSPCM are randomly generated for each iteration. Also, the geometric average and not the arithmetic average is needed for ka because the values span a range of orders of magnitude. Functions: findka clear all format short e global rH Lfond g N rR4 d iterations = 1000; % Number of iterations used to calculate ka % for each polymer m = 20; % Maximum number of binding units on HBSPCM rH = 4.0e−9; % Radius of HBSPCM [18] [m] a = 8.5e−10; % Minimal distance of approach [m]. Value % from Chapter 4. rR4 = 6.613e−10; % Radius of binding unit R4−1 [m]. 3D % representation. Value from Chapter 4. bd = a; % Binding distance between fondaparinux % and binding unit R4−1. Measured with % respect to centroid of the molecules. % Value from Chapter 4. d = sqrt(2*(bd+rR4)ˆ2); % Diagonal of binding site [m]. % 3D representation. Value from MS. Lfond = 19.939e−10; % Length of fondaparinux of fondaparinux [m]. % 3D representation. Value from Chapter 4. dcut = 9.160e−10; % Minimal distance between binding units [m]. % Found using Program D.1 and % experimental data units, R4−1. g = 3; % Number of charges grouped for fondaparinux 129 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model % binding N = [3:m]; % Range of number of R4−1 binding units that % will be used to generate HBSPCMs bound = NaN(iterations,floor(m/3)+1,m−g+1); count = 0; for i = 1:iterations [ka(i,:),f(i,:)] = findka(m,a,dcut); count = count + 1; % Stores number of R groups per binding site for each % polymer and iteration for k = 1:m−g+1 for j = 1:length(R bound temp(k,:)) bound(i,j,k) = R bound temp(k,j); end end end 130 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model Program D.3 This program calculates ka using the Equations 2.2-2.8 for HBSPCMs over a range of effective charges for the R6-1 binding units. At each effective charge ka is calculated by calling the function findka varycH which calculates ka for HBSPCMs containing a varying number of R6-1 binding units attached to the molecule. The program executes a set number of iterations and from the results calculates the geometric average and standard deviation for ka and the arithmetic average and standard deviation for the number of fondaparinux molecules bound (f) for each HBSPCM tested. Average values of ka and f are needed because the locations of the R4-1 binding units on the surface of HBSPCM are randomly generated for each iteration. Also, the geometric average and not the arithmetic average is needed for ka because the values span a range of orders of magnitude. Functions: findka varycH clear all format short e global rH Lfond g N rR4 d iterations = 1000; % Number of iterations used to calculate ka % for each polymer m = 20; % Maximum number of binding units on HBSPCM rH = 4.0e−9; % Radius of HBSPCM [18] [m] a = 8.5e−10; % Minimal distance of approach [m]. Value % from Chapter 4. rR6 = 7.099e−10; % Radius of binding unit R6−1 [m]. 3D % representation. Value from Chapter 4. bd = a; % Binding distance between fondaparinux % and binding unit R6−1. Measured with % respect to centroid of the molecules. % Value from Chapter 4. d = sqrt(2*(bd+rR6)ˆ2); % Diagonal of binding site [m]. % 3D representation. Value from MS. Lfond = 19.939e−10; % Length of fondaparinux of fondaparinux [m]. % 3D representation. Value from Chapter 4. dcut = 9.160e−10; % Minimum distance between binding units % [m]. Value used for R4−1 in Program D.2. 131 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model g = 3; % Number of charges grouped for fondaparinux % binding N = [3:m]; % Range of number of R6−1 binding units that % will be used to generate HBSPCMs num appcH = 25; % Number of effective charges to be tested appcH = linspace(0.1,2.29,num appcH); % Effective charge of % R6−1 binding unit % [atomic units] length appcH = length(appcH); for i = 1:length appcH count appcH = i; for j = 1:iterations count iterations = j; [ka iter(j,:),f iter(j,:)] = ... findka varycH(m,a,dcut,appcH(i)); end ka avg = geomean(ka iter); std ka = geostd(ka iter,0); ka appcH(i,:) = ka avg(1,:); std ka appcH(i,:) = std ka(1,:); f avg = mean(f iter); std f = std(f iter); f appcH(i,:) = f avg(1,:); std f appcH(i,:) = std f(1,:); se f appcH(i,:) = std f appcH(i,:)/sqrt(iterations); bar f appcH(i,:)= 1.96*se f appcH(i,:); end save ka vary appcH.out ka appcH −ASCII save stdka vary appcH.out std ka appcH −ASCII save f vary appcH.out f appcH −ASCII save stdf vary appcH.out std f appcH −ASCII 132 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model findka varycH The function findka varycH calculates ka for HBSPCMs with a varying number of R6-1 binding units attached to their surface. In order to do so the function calls the functions: fonda energy to calculate the internal energy of the fonda- parinux molecule, gen charges and poly energy to construct an HBSPCM and calculate its internal energy respectively, binding R to determine fondaparinux binding to binding units, and complex energy to calculate the internal energy of the bound complex. In the case that HBSPCM contains only three binding units and fondaparinux binds to it, the vectors containing the x,y,z coordi- nates of the binding units (i.e. Xm, Ym, and Zm) contain only one element and therefore Ucomplex = 0. Input arguments: m = maximum number of binding units on HBSPCM a = minimal distance of approach [m] dcut = minimal distance between binding units, R6-1 [m] appcH = apparent charge of the binding units, R6-1 Output arguments: ka = m x i matrix of values for ka, where i is the number of iterations set by the program calling this function f = number of fondaparinux molecules bound to one HBSPCM Functions: fonda energy gen charges poly energy binding R complex energy function [ka,f] = findka varycH(m,a,dcut,appcH) global rH K kao eo er Lfond N R bound temp R bound Ufond Upoly ... Ucomplex DeltaU Wfond = 10.467e−10; % Width of fondaparinux, 3D [m]. Value from % Chapter 4. rF = 11.812e−10; % Radius of fondaparinux, 3D [m]. Average % of distances measured from O19−O53 and % O19−O62. Value from Chapter 4. cH = 1.29; % Elementary charge of each binding % unit, R4−1 [atomic units]. Value % from Chapter 5.3.1. 133 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model qH = cH*1.602176487e−19; % Charge of each binding % unit [C] cF = −1; % Elementary charge of one % fondaparinux molecule % [atomic units]. Value from % [77]. qF = cF*1.602176487e−19; % Charge of each fondaparinux % molecule [C] % DEBYE−HUCKEL PARAMETER, K [1/m]: F = 96485.34; % Faraday's constant [C/mol]. Value from [99]. I = 150; % Ionic strength of solution (physiological) % [mM = mol/mˆ3] eo = 8.854e−12; % Vacuum permittivity [Cˆ2/(J*m)]. Value % from [99]. er = 80; % Dielectric constant at T for water % [unitless]. Value from [99]. R = 8.314472; % Gas constant [J/(K*mol)]. Value from [99]. T = 25 + 273.15; % Temperature of solution (experimental) [K] K = sqrt((2*(Fˆ2)*I)/(eo*er*R*T)); % Debye−Huckel parameter [1/m] % BASAL RATE CONSTANT, kao [1/(M*s)]: visc = 0.90e−3; % Viscosity of the solvent (i.e. water at % 25oC) [Pa*s]. % Value from [99]. kB = 1.3806504e−23; % Boltzmann constant [J/K]. Value from [99]. NA = 6.02214179e23; % Avogadro constant [1/mol]. Value from [99]. DH = (kB*T)/(6*pi()*visc*rH); % Diffusion constant of HBSPCM DF = (kB*T)/(6*pi*visc*rF); % Diffusion constant of % fondaparinux D = DH + DF; % Diffusion constant of % complex Rint = Lfond/(log((2*Lfond)/Wfond));% Interaction radius [m] kao = 4*pi*NA*D*Rint; % Basal rate constant % [1/(M*s)] % DEBYE−HUCKEL ENERGIES Ufond = fond energy(a,qF); 134 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model n = length(N); ka = zeros(1,n); R bound temp = NaN(n,floor(m/3)); for k = 1:n numbercharges = N(k); X = []; Y = []; Z = []; [X,Y,Z] = gen charges(N(k),dcut); Upoly(k) = poly energy(X,Y,Z,N(k),a,qH); Xm = []; Ym = []; Zm = []; chargem = []; [Xm,Ym,Zm,chargem,f temp] = binding R(X,Y,Z,N(k),qH,qF); for j = 1:length(R bound) R bound temp(k,j) = R bound(j); end f(k) = f temp; Ucomplex(k) = complex energy(a,Xm,Ym,Zm,chargem); if f(k) ˜= 0 DeltaU(k) = Ucomplex(k) − Upoly(k) − f(k)*Ufond; ka(k) = kao*exp(−(DeltaU(k)/(kB*T))*(1/(1+K*a))); else ka(k) = 1; end end 135 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model binding R The function binding R checks if binding units, R6-1, attached to the surface of an HBSPCM are within a specified distance (d) from each other. If this is the case, they are considered a binding site and are assigned a charge that reflects binding to fondaparinux. A binding site contains three or more binding units. Input arguments: X = x-coordinates of charges on HBSPCM - must be vector Y = y-coordinates of charges on HBSPCM - must be vector of same size as X Z = z-coordinates of charges on HBSPCM - must be vector of same size as X n = number of charges on HBSPCM. Scalar. qH = charge of each R group [C] qF = charge of each fondaparinux molecule [C] Output arguments: X c = x-coordinates of charges on HBSPCM - must be vector Y c = y-coordinates of charges on HBSPCM - must be vector of same size as Xm Z c = z-coordinates of charges on HBSPCM - must be vector of same size as Xm charge c = charges of X c,Y c,Z c - must be a vector f = number of fondaparinux molecules bound to one HBSPCM function [X c,Y c,Z c,charge c,f] = binding(X,Y,Z,n,qH,qF) global rH d arc length R bound arc length = zeros(n,n); v = inf*ones(1,n); arc length = arc length − diag(diag(arc length))+diag(v); for i = 1:n for j = i+1:n dist(i,j) = sqrt((X(i)−X(j))ˆ2 + (Y(i)−Y(j))ˆ2 ... +(Z(i)−Z(j))ˆ2); arc length(i,j) = rH*2*asin((dist(i,j))/(2*rH)); arc length(j,i) = arc length(i,j); end end % Randomly selects a binding unit, R6−1, on the surface of an 136 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model % HBSPCM and checks for other R6−1 binding units that are within % a specified distance, d, to form a binding site. It then % creates vectors with the x,y,z coordinates of the R6−1 binding % units that formed a binding site and vectors with the x,y,z % coordinates of the R6−1 binding units that did not form a % binding site. start = randperm(n); Rs = 1; f = 0; nf = 0; X b3 = []; Y b3 = []; Z b3 = []; charge b3 = []; X nb = []; Y nb = []; Z nb = []; charge nb = []; R used = [1:n]; count = 0; R bound = []; nb bin = 0; for i = 1:n loc = []; if R used(start(i)) > 0 for j = 1:n if R used(j) > 0 if arc length(start(i),j) <= d R used(start(i)) = 0; R used(j) = 0; loc = [loc j]; Rs = Rs + 1; end end end 137 Appendix D. MATLAB Code for the Improved Electrostatic Interactions Model % If a binding site is formed by three or more R6−1 % binding units, it gets assigned the coordinates from % the binding unit used as the reference point to find % the other binding units 'd' or less away from it. % The charge of the binding site will depend on the % number of binding units present. % Rs = number of binding units in a binding site. if Rs >= 3 count = count + 1; R bound(count) = Rs; f = f + 1; X b3(f) = X(start(i)); Y b3(f) = Y(start(i)); Z b3(f) = Z(start(i)); charge b3(f) = qF − qH*Rs; else nf = nf + 1; X nb(nf) = X(start(i)); Y nb(nf) = Y(start(i)); Z nb(nf) = Z(start(i)); nb bin(nf) = start(i); charge nb(nf) = qH; R used(start(i)) = start(i); R used(loc) = loc; end Rs = 1; end end % Fondaparinux−HBSPCM complex: % Combines the vectors of the x,y,z coordinates of both the % binding units inside the binding sites and those outside % the binding site. X c = [X b3 X nb]; Y c = [Y b3 Y nb]; Z c = [Z b3 Z nb]; charge c = [charge b3 charge nb]; 138
Thesis/Dissertation
eng
Biomedical Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Predicting binding affinities of a novel polymer for the neutralization of fondaparinux
Text