M A S A R Y K U N I V E R S I T Y F A C U L T Y O F S C I E N C E C E N T R A L E U R O P E A N I N S T I T U T E O F T E C H N O L O G Y Protein—protein interactions important in neurodegenerative diseases D o c t o r a l T h e s i s Petr Louša Supervisor: doc. RNDr. Jozef Hritz, Ph.D. UNI £ * < = ^ ' Brno 2022 BIBLIOGRAPHIC E N T R Y Author Title of the Thesis Degree Programme Study Plan Supervisor Academic Year Number of Pages Keywords Mgr. Petr Lousa Faculty of Science, Masaryk University Central European Institute of Technology Protein-protein interactions important in neurodegenerative diseases Life Sciences Structural Biology doc. RNDr. Jozef Hritz, Ph.D. 2021/2022 xiv +116 protein, protein-protein interactions, 14-3-3, tyrosine hydroxylase, dimer, monomer, fluorescence, FRET, MST, NMR, kinetics, dissociation constant, phosphorylation BIBLIOGRAFICKÝ ZÁZNAM Autor Mgr. Petr Louša Prírodovedecká fakulta, Masarykova univerzita Středoevropský technologický institut Název práce Protein­proteinové interakce významné při neurodegenerativních chorobách Studijní program Vědy o živé přírodě Studijní plán Strukturní biologie Vedoucí práce doc. RNDr. Jozef Hritz, Ph.D. Akademický rok 2021/2022 Počet stran xiv + 116 Klíčová slova protein, protein­proteinove interakce, 14­3­3, tyrosin hydroxylasa, dimer, monomer, fluorescence, FRET, MST, NMR, kinetika, disociacni konstanta, fosforylace Abstract The family of 14-3-3 proteins is present in all eukaryotic organisms. The 14-3-3 proteins serve as important regulatory hubs interacting with hundreds of protein partners typically through phosphorylated intrinsically disordered regions. The family consists of seven isoforms in humans that readily form homodimers and heterodimers differing in function and specificity. Although the dimerization is a crucial feature of 14-3-3 proteins, the biophysical properties of the dimerization were not described in quantitative manner previously. Moreover, the dimer-monomer equilibria of 14-3-3 proteins are significantly affected by phosphorylation of serine 58 (numbering according the human £ isoform) located on the dimeric interface. This thesis presents the design and analysis of a set of assays providing quantitative characterization of the 14-3-3 dimermonomer equilibrium, evaluating the values of rate constants, dissociation constants, and thermodynamic parameters including Gibbs free energy, enthalpy and entropy. The presented experiments were performed on the human £ isoform, its S58-phosphorylated variant, and the a isoform. However, it is readily possible to extend the developed methods and analysis to other 14-3-3 isoforms and also to other protein families forming tight dimers. The theme of protein-protein interactions is further extended in the other projects covered in this thesis. The tyrosine hydroxylase is an enzyme catalysing the first step in a pathway leading to dopamine and adrenaline. It interacts with 14-3-3 proteins through its regulatory domain that gets phosphorylated at multiple sites located within the intrinsically disordered N-terminal region. Using non-uniformly sampled time-resolved nuclear magnetic resonance, we quantitatively described the phosphorylation kinetics. We studied the impact of phosphorylated regulatory domain on the 14-3-3 dimer-monomer equilibrium. The structural features of the regulatory domain characterized by NMR relaxation experiments were employed to benchmark force fields used in molecular dynamics simulations. Abstrakt Proteiny z rodiny 14­3­3 se vyskytují u všech eukaryotických organismů, kde regulují významnou část buněčných pochodů skrze interakce se stovkami vazebných partnerů, obvykle prostřednictvím fosforylovaných neuspořádaných regionů. Tato rodina zahrnuje v lidském organismu sedm isoforem, které interagují za tvorby homodimerů a heterodimerů, jež se vzájemně liší svou funkcí i specificitou. I přes klíčovou roli dimerizace pro proteiny 14­3­3, její biofyzikálni parametry doposud nebyly v literatuře kvantifikovány. Dimermonomerní rovnováha je navíc silně ovlivněna fosforylací serinu 58 (číslováno dle lidské £ isoformy) na dimerním rozhraní. V této práci představujeme design a analýzu sady experimentů, umožňujících kvantitativní analýzu dimer­monomerní rovnováhy proteinů 14­3­3. Určili jsme hodnoty rychlostních a rovnovážných konstant, Gibbsovy volné energie, entropie a entalpie. Hlavním objektem zájmu byla lidská isoforma £ a její varianta fosforylovaná na serinu 58, několik experimentů bylo provedeno s využitím isoformy a. Ukazujeme také, že popsané postupy je možné relativně snadno využít pro další isoformy 14­3­3, ale i jiné rodiny proteinů, tvořící dimery s vysokou afinitou. Téma protein­proteinových interakcí se objevuje i v dalších popsaných projektech. Tyrosin hydroxylasa je enzym, který katalyzuje první krok v syntéze dopaminu a adrenalinu. Interakce s proteiny 14­3­3 je zprostředkována pomocí regulační domény, která může být fosforylovaná na několika místech v neuspořádané části. Pomocí neuniformně vzorkované, časově rozlišené nukleární magnetické resonance jsme kvantitativně charakterizovali kinetiku této fosforylace. Rovněž jsme určovali vliv vazby fosforylované regulační domény na dimer­monomerní rovnováhu 14­3­3. V dalším projektu jsme využili data získaná z relaxačních NM R experimentů, provedených na regulační doméně, pro porovnání biomolekulárních silových polí využívaných v simulacích molekulární dynamiky. (q) Petr Louša, Masaryk University, 2022 Acknowledgement I would like to thank Jožo who has been a great supervisor. You always were there whenever I needed you, full of patience and understanding. I am happy I could join your lab, it has been a wonderful time working in your team. I would like to appreciate Lukáš, his patience and explanatory skills. Thank you for the opportunity to learn from you and also to teach students with both of you. The next big thanks go to all the people in our team. It has been a great time and I really enjoyed it. Among others, thank you Zuzka for teaching me to drink coffee and being a wonderful friend; Hany, Verča, Alča and Tom for filling me with energy when I lacked it; Jano for having fun solving difficult problems. I thank to all of the members of our group, especially Krishna, Anet, M arkéta, Noro, Radek and all others. I am thankful to Gabi Žoldák for his insights into the fluorescence experiments and for taking care of me, while I did my internship in Munich. I am deeply grateful for my beloved M íša. I thank you for everything you do, for supporting me all the time and especially for making me a better person. Na závěr chci poděkovat mým rodičům a celé rodině. M ami a tati, děkuji vám za všechno, co jste pro mě udělali a obětovali, abych mohl studovat. Díky za všechnu vaši podporu a lásku, bez vás bych to nikdy nedokázal. Věnováno rodičům Contents Abstract iv Abstrakt v Contents ix List of Figures xi List of Tables xiii List of Abbreviations xiv 1 Introduction 1 1.1 Two are better than one 1 1.2 14-3-3 protein family 3 1.2.1 Functions of 14-3-3 proteins 3 1.2.2 Interactions of 14-3-3 proteins 6 1.3 Tyrosine hydroxylase 7 1.3.1 Function and structure 7 1.3.2 Regulatory domain 8 1.3.3 Interaction with 14-3-3 proteins 10 1.4 Intrinsically disordered proteins 10 1.5 Interaction modelling 12 1.5.1 Homo- and heterodimers 13 2 Aims of the Thesis 17 3 Methods 18 3.1 Fluorescence assays 18 ix Contents 3.1.1 Sample preparation 18 3.1.2 Experimental setup 19 3.1.3 Model description 20 3.1.4 Error analysis 24 3.1.5 Detailed description of COPASI workflow . . 24 3.2 NMR experiments 33 3.2.1 Homodimerization of pS58_ 14-3-3£ 33 3.2.2 Regulatory domain of human tyrosine hydroxylase 1 34 4 Results and Discussion 36 4.1 14-3-3 dimerization 37 4.1.1 Equilibrium assays 38 4.1.2 Kinetic assays 41 4.1.3 Error analysis 42 4.1.4 Effect of external conditions 44 4.1.5 Effect of phosphorylation at S58 46 4.1.6 14-3-3a dimerization 49 4.1.7 Thermodynamic parameters 49 4.1.8 Heterodimerization remarks 55 4.2 Regulatory domain of hTHl 57 4.2.1 Phosphorylation of the unstructured part . . 57 4.2.2 Relaxation properties 60 5 Conclusions 62 Bibliography 64 Curriculum Vitae 78 List of Publications 80 List of Presentations 82 Paper 1 83 Paper 2 97 Paper 3 103 x List of Figures 1.1 Examples of functional consequences of oligomerization 2 1.2 Multiple sequence alignment of 14-3-3 isoforms . . . 4 1.3 The structure of 14-3-3^ protein dimer 5 1.4 A model of the dimeric regulatory domain of tyrosine hydroxylase 9 3.1 Screenshot of COPASI showing the Events dialog window 26 3.2 Screenshot of COPASI showing the Global Quantities dialog window 27 3.3 Screenshot of COPASI showing the Experimental Data dialog from the Parameter Estimation window. . . . 28 3.4 Screenshot of COPASI showing the Parameter Estimation window 30 3.5 Screenshot of COPASI showing the Parameter Scan window with the Report dialog open 32 4.1 The 14-3-3C homodimerization equilibrium assays . . 40 4.2 The 14-3-3C dimerization kinetic assays 43 4.3 The effect of external conditions on the Q/Q dimerization 45 4.4 Establishing the homodimerization K<± of p£ 47 4.5 Determination of pC/C heterodimerization dissociation constant 48 4.6 Characterization oia/a homo- and a/( heterodimerization 50 4.7 The effect of temperature on the £/£ and a/a dimerization 53 xi List of Figures 4.8 Phosphorylation of RD-hTHl as monitored by timeresolved NMR 59 xii List of Tables 4.1 Thermodynamic equilibrium parameters of 14-3-3 dimerization 54 4.2 Thermodynamic kinetic parameters of 14-3-3 dimerization 54 4.3 Example concentrations of 14-3-3^ forms after phosphorylation 56 4.4 Example concentrations of 14-3-3^ and a 56 xiii List of Abbreviations a 14-3-3cr protein c 14-3-3^ protein C _ A 14-3-3^ protein selectively labelled with AF647 C_D 14-3-3C protein selectively labelled with AF488 C Q 14-3-3C protein selectively labelled with AF488 PC 14-3-3^ protein phosphorylated at pS58 AF488 AlexaFluor488-C2-maleimide AF647 AlexaFluor647-C5-maleimide FRET Förster resonance energy transfer bTHl human tyrosine hydroxylase 1 IDP intrinsically disordered protein L-DOPA L-3,4-dihydroxyphenylalanine MST microscale thermophoresis NMR nuclear magnetic resonance NUS non-uniform sampling PAGE Polyacrylamide gel electrophoresis PKA cAMP-dependent protein kinase PRAK p38 regulated/activated protein kinase RD-hTHl regulatory domain of human tyrosine hydroxylase 1 RDC residual dipolar coupling SAXS small angle X-ray scattering SDS-PAGE sodium dodecyl sulfate-polyacrylamide gel electro- phoresis SQ self-quenching SSE sum of squared errors ssNOE steady-state 1 H - 1 5 N heteronuclear Overhauser en- hancement TCEP tris(2-carboxyethyl)phosphine TMR 6-tetramethyl-rhodamine-C6-maleimide xiv Chapter 1 Introduction 1.1 Two are better than one The majority of proteins is found forming some quaternary structure. Monomeric proteins account for less than 20 % of all proteins in Escherichia coli, while the ratio is even lower in higher organisms [1, 2]. The most common pattern are homodimers, which account for almost 40 % of all proteins in E. coli [1]. Even the proteins, that are monomeric by themselves, very often aggregate with other proteins to form complexes or engage in signalling cascades where the signal is transmitted by sequential interactions. Thus, the importance of protein-protein interactions is immense and ubiquitous. The fact of dimer formation might be interesting by itself, however, the process of association and dissociation of the dimer and its regulation is often important. As the proteins tend to have different properties and functions in monomeric and dimeric state, the regulation of the dimers creation and their dissociation back into monomers has to reflect this complexity. The difference between the states varies from mere increased stability of dimers due to reduced surface area; through changes in binding properties of substrates, including cooperativity as observed e.g. in haemoglobin; creation of new or shielding the existing active sites; to the increase of diversity in binding other protein partners into regulatory complexes [3]. Examples of particular functions of oligomerization are illustrated in Fig. 1.1. 1 Introduction 1.1 Two are better than one Figure 1.1: Examples of functional consequences of oligomerization. (A) After dimerization, the solvent accessible area is reduced and the stability increases, as indicated by the darker colour. (B) There are two binding sites on the orange protein. After dimerization, both can be utilized and involved in formation of a tertiary complex. (C) In the tetramer, binding of a substrate or a cofactor leads to a change in conformation of one subunit, indicated by the darker colour. Due to allostery, all other subunits are influenced structurally. (D) The dimerization alters the structure and creates a new active site. (E) The active site gets inaccessible after formation of the dimer. Adapted from [3]. 2 Introduction 1.2 14-3-3 protein family 1.2 14-3-3 protein family The 14-3-3 protein family is a good example of a complex dimeric behaviour in a network of interacting molecules. The family consists in human organism of seven isoforms labelled by Greek letters - j3, 7, £, £ i Q [4]- The isoforms are highly homologous (see Fig. 1.2) and form both homo- and heterodimers in native conditions. The dimerization tendencies differ among the isoforms - for instance, the a isoform is found almost exclusively in homodimeric state, the 7 forms homodimers and heterodimers with s and the e isoform itself is promiscuous forming mostly heterodimers with other isoforms. These tendencies were characterized by immunoprecipitation [5, 6] and mass spectrometry [7], however no quantitative dissociation constants were provided. The 14-3-3 proteins are very important regulatory proteins found throughout the whole animal and plant kingdoms [8]. In human, they are present in all tissues [9, 10]. The highest expression levels are found in human brain, where about 1 % of soluble proteins are members of the 14-3-3 family [11, 12]. The 3D structure, solved crystallograficaHy for each of the 7 human isoforms, contains nine a-helices arranged into a wide „u/' shape (Figure 1.3). 1.2.1 Functions of 14-3-3 proteins The 14-3-3 proteins form a central hub within the human interactome. There is more than 1200 reported protein partners and the number is expected to grow further [14, 15]. The 14-3-3 proteins are involved in most cellular processes including signal transduction [16, 17], cell cycle and apoptosis [18, 19], cytoskeletal organisation, stress response, DNA replication, viral infection and many others [20-22]. The functions of 14-3-3 in human brain involve among others neuronal signalling [23, 24], neurotransmitter synthesis [25, 26], or cytoskeletal stability and microtubule associated proteins (e.g. tau) regulation [27, 28]. The 14-3-3 involvement in various neurological diseases is studied extensively. Co-localization with parkin and a-synuclein was observed in brains of patients with Parkinson's disease [29, 30]. Significantly lowered 14-3-3 levels were found in patients suffering with schizophrenia and multiple sclerosis [31, 32]. The association of 14-3-3 with the neurofibrillary tangles, formed mostly by hy- 3 Introduction 1.2 14-3-3 protein family HI H2 H3 * mwuuww nwuuwi nMmxmwmviD e p s i i o n -MDDREDLVYQAKLAEQAERYDEMVESMKKVAGMDVELTVEERNLLSVAYKNVIGARRAIŠ Sigma —MERASLIQKAKLAEQAERYEDMAAFMKGAVEKGEELSCEERNLLSVAYKNWGGQRAA gamma -MVDREQLVQKARLAEQAERYDDMAAAMKNVTELNEPLSNEERNLLSVAYKNWGARRSS eta -MGDREQLLQRARLAEQAERYDDMASAMKAVTELNEPLSNEDRNLLSVAYKNWGARRSS t h e t a —MEKTELIQKAKLAEQAERYDDMATCMKAVTEQGAELSNEERNLLSVAYKNVVGGRRSA beta MTMDKSELVQKAKLAEQAERYDDMAAAMKÄVTEQGHELSNEERNLLSVAYKNWGARRSS zeta —MDKNELVQKAKLAEQAERYDDMAACMKSVTEQGAELSNEERNLLSVAYKNWGARRSS * . .*.********..* * * * . *.***********.* . * . . 59 58 59 59 58 60 58 H3 H4 H5 HUMM« ^mtmmm^utmvvnv mi e p s i i o n WRIISSIEQKEENKGGEDKLKMIREYRQMVETELKLICCDILDVLDKHLIPAAN—TGES 117 Sigma WRVLSSIEQKSNEEGSEEKGPEVREYREKVETELQGVCDTVLGLLDSHLIKEAG—DAES 116 gamma WRVISSIEQKTSADGNEKKIEMVRAYREKIEKELEAVCQDVLSLLDNYLIKNCSETQYES 119 eta WRVISSIEQKTMADGNEKKLEKVKAYREKIEKELETVCNDVLSLLDKFLIKNCNDFQYES 119 t h e t a WRVISSIEQKT—DTSDKKLQLIKDYREKVESELRSICTTVLELLDKYLIANAT—NPES 114 beta WRVISSIEQKT—ERNEKKQQMGKEYREKIEAELQDICNDVLELLDKYLIPNAT—QPES 116 zeta WRWSSIEQKT—EGAEKKQQMAREYREKIETELRDICNDVLSLLEKFLIPNAS—QAES 114 .****** * * . . * * * H5 H6 H7 mmwwwmT um»nmr\\\i\\vi(i9 vmww e p s i i o n KVFYYKMKGDYHRYLAEFATGNDRKEAAENSLVAYKAASDIAMTELPPTHPIRLGLALNF 177 Sigma RVFYLKMKGDYYRYLAEVATGDDKKRIIDSARSAYQEAMDISKKEMPPTNPIRLGLALNF 17 6 gamma KVFYLKMKGDYYRYLAEVATGEKRATWESSEKAYSEAHEISKEHMQPTHPIRLGLALNY 17 9 eta KVFYLKMKGDYYRYLAEVASGEKKNSWEASEAAYKEAFEISKEQMQPTHPIRLGLALNF 17 9 t h e t a KVFYLKMKGDYFRYLAEVACGDDRKQTIDNSQGAYQEAFDISKKEMQPTHPIRLGLALNF 174 beta KVFYLKMKGDYFRYLSEVASGDNKQTTVSNSQQAYQEAFEISKKEMQPTHPIRLGLALNF 17 6 zeta KVFYLKMKGDYYRYLAEVAAGDDKKGIVDQSQQAYQEAFEISKKEMQPTHPIRLGLALNF 174 .*** ****** ***.* * *. . ** * .*. . **.*********. H7 H8 H9 e p s i i o n SVFYYEILNSPDRACRLAKAAFDDAIAELDTLSEESYKDSTLIMQLLRDNLTLWTSDMQG 237 Sigma SVFHYEIANSPEEAISLAKTTFDEAMADLHTLSEDSYKDSTLIMQLLRDNLTLWTADNAG 236 gamma SVFYYEIQNAPEQACHLAKTAFDDAIAELDTLNEDSYKDSTLIMQLLRDNLTLWTSDQQD 239 eta SVFYYEIQNAPEQACLLAKQAFDDAIAELDTLNEDSYKDSTLIMQLLRDNLTLWTSDQQD 239 t h e t a SVFYYEILNNPELACTLAKTAFDEAIAELDTLNEDSYKDSTLIMQLLRDNLTLWTSDSAG 234 beta SVFYYEILNSPEKACSLAKTAFDEAIAELDTLNEESYKDSTLIMQLLRDNLTLWTSENQG zeta SVFYYEILNSPEKACSLAKTAFDEAIAELDTLSEESYKDSTLIMQLLRDNLTLWTSDTQG .* ** * . * * * * * * * * * * * * * * * * * * * * . .* * * . * * * * *. * *** . * * e p s i i o n DGEEQNKEALQDVEDENQ 255 sigma EEGGEAPQEPQS 248 gamma DDGGEGNN 247 eta EEAGEGN 246 theta EECDAAEGAEN 245 beta DEGDAGEGEN 246 zeta DEAEAGEGGEN 245 236 234 Figure 1.2: Multiple sequence alignment of seven human 14-3-3 isoforms showing high sequence homology. Conserved residues are highlighted in cyan, strongly similar residues are highlighted in orange. The a-helices H1-H9 are illustrated above sequences. The phosphorylatable serine S58 is indicated by a star and a bounding box. 4 Introduction 1.2 14-3-3 protein family Figure 1.3: The structure of 14-3-3^ protein dimer containing two noncovalently bound monomeric units. The monomeric units are distinguished by colour, the helices are labelled H1-H9, the ends of the chain are indicated by letters N (only upper view) and C (only lower view). The lower view contains the phosphorylatable serines S58 (highlighted by a dashed ellipse) including the phosphates for illustration purposes. Based on structure with PDB code 1A40 [13] with truncated C-termini. 5 Introduction 1.2 14-3-3 protein family perphosphorylated tau protein [33], is found in Alzheimer disease patients' brains [34-36]. The 14-3-3 proteins are suspected to be involved in the phosphorylation of tau [37-39]. However, the exact mechanisms behind the role of 14-3-3 in these diseases are yet to be revealed. While the 14-3-3 dimer-monomer equilibrium is dynamic and strongly dependent on the isoform and its total concentration [7], the majority of the protein is in the dimeric form at physiological conditions. However, the monomers do have different properties than the dimers, for instance it was reported that they exhibit higher chaperone-like activity than the dimers [40, 41], are unable to regulate Raf kinase activity [42], or modulate the activity of ion channels [43]. In vivo, the monomerization can be achieved by phosphorylation of serine S58 located at the dimeric interface [44]. This serine is conserved in five out of the seven human isoforms (/?, 7, e, ?7, C, see Fig. 1.2) [45]. Thus, phosphorylation plays a principal role in 14-3-3 functioning - it enables the binding of phosphorylated protein partners and also modifies the functions of 14-3-3 itself. 1.2.2 Interactions of 14-3-3 proteins The 14-3-3 proteins contain in each monomer a conserved amphipathic groove that binds phosphorylated proteins featuring one of the three consensus motifs: RSXpSXP, RX(Y/F)XpSXP, and a carboxyl-terminal pSX-COOH, where X stands for a non-proline residue and pS for a phosphoserine or phosphothreonine residue [46-48]. However, there are many protein partners that do not follow this consensus. There are several proposed regulation mechanisms that rely on the dimeric nature of 14-3-3 featuring two interaction sites, one in each subunit. The „clamp" [49], or „molecular anvil" [45] hypothesis suggests grasping a bivalent client protein into the two grooves and physically deforming it while the 14-3-3 itself changes only insignificantly. In the „adaptor" mechanism, two different clients bind into the 14-3-3 dimers to enhance their interaction. The binding sites are very often found within intrinsically disordered regions of the client proteins [50]. The disordered state is important for several reasons. Firstly, the flexible region is easily accessible by kinases that provide the phosphorylation needed for interaction with 14-3-3. Secondly and possibly more importantly, 6 Introduction 1.3 Tyrosine hydroxylase the disorder allows the vast number of diverse proteins to interact with structurally well defined 14-3-3 [51]. Following the binding to 14-3-3 protein, the previously disordered interacting region may undergo a transition to an ordered state [50]. Despite the loss of entropy due to this transition, the overall free energy of binding is negative as a result of a large favourite enthalpy change. This balance effectively uncouples the binding affinity from the specificity, leading to protein-protein interactions that are simultaneously reversible and highly specific, which is necessary for cell signalling and regulation. As the binding partners almost always feature a phosphorylated serine or threonine, it is very convenient to employ the 3 1 P NMR spectroscopy to assess the binding properties. The high sensitivity of phosphorus NMR, combined with the simplicity of ID spectra that usually contain only a handful of well resolved peaks (contrary to more common proton spectra containing thousands of overlapping signals), allows to investigate even complicated interaction scenarios of multimeric proteins phosphorylated at multiple sites, such as a dimeric regulatory domain of human tyrosine hydroxylase [52]. Moreover, the monomerization effect of phosphorylation at S58 of 14-3-3 can also be conveniently studied using the phosphorus spectra. 1.3 Tyrosine hydroxylase The following section describes a 14-3-3 interaction partner, human tyrosine hydroxylase 1 (hTHl), with a special focus on its regulatory domain (RD-hTHl) through which it interacts with 14-3-3 proteins. We studied this protein on its own (Paper 2, Paper 3) as well as its interaction with 14-3-3 (Paper 1). 1.3.1 Function and structure Tyrosine hydroxylase is an enzyme catalysing hydroxylation of amino acid tyrosine into L-3, 4-dihydroxyphenylalanine (L-DOPA). This conversion constitutes the first and rate-limiting step in the biosynthetic pathway leading to catecholamine neurotransmitters dopamine, noradrenaline, and adrenaline [53, 54]. Within human body, tyrosine hydroxylase is located in the brain and adrenal glands. The lack or deficiency of T H activity leads 7 Introduction 1.3 Tyrosine hydroxylase to severe diseases including progressive encephalopathy, infancy parkinsonism, or L-DOPA responsive dystonia [55-57]. Links of TH activity to neurological diseases like schizophrenia and bipolar affective disorder were also suggested [58, 59]. There are four isoforms of human tyrosine hydroxylase, originating from a single-copy gene due to alternative splicing of the corresponding mRNA [60, 61]. The difference lies mostly within the regulatory domain. While the isoform 1 is the shortest, isoforms 2, 3, and 4 contain inserts of 4, 27, and 31 residues, respectively. The inserts follow residue M30, directly preceding S31 of isoform 1, which can be phosphorylated [62]. This suggests that the isoforms differ in their regulation properties and might be recognized by different kinases. The human tyrosine hydroxylase forms homotetrameric quaternary structure [63, 64]. Each subunit comprises three domains: an N-terminal regulatory domain (residues 1-169 in isoform 1), a catalytic domain (170-450), and a short C-terminal tetramerization domain (451-497) [65, 66]. The tetramer is arranged as a dimer of dimers. 1.3.2 Regulatory domain The regulatory domain of human tyrosine hydroxylase 1 (RD-hTHl) contains 169 residues. The N-terminal 65 residues form an intrinsically disordered region, which plays an important role in regulation. The C-terminal part, containing about 100 residues, adopts ACT domain fold and forms a stable homodimer [67]. The structure is shown in Fig. 1.4. The regulation of tyrosine hydroxylase catalytic activity is achieved by the phosphorylation of the disordered region of the regulatory domain [69, 70]. There are four phosphorylation sites within the RD-hTHl - threonine T8 and serines S19, S31, and S40. While the serine sites are phosphorylated in vivo by several kinases and do serve important regulatory purposes, the threonine T8 appears to have no effect. In the inactive state of hTHl, the flexible N-tail is presumably covering the active site and physically blocking it [71, 72]. The phosphorylation (on S40) leads to release of the N-tail from the active site uncovering it and thus increasing the enzymatic activity. This effect can be enhanced and stabilised by interaction with 14- 8 Introduction 1.3 Tyrosine hydroxylase Figure 1.4: A model of the dimeric regulatory domain of tyrosine hydroxylase. The structured part was solved by NMR (PDB code 2MDA) [67]. The unstructured part was generated in an extended state using Modeller software [68] containing a transient helix (residues 42-55) and two phosphorylated serines (labelled and highlighted by ellipses) at positions S19 and S40. 9 Introduction 1.4 Intrinsically disordered proteins 3-3 protein that also protects the phosphoserines from phosphate removal by phosphatases [25]. 1.3.3 Interaction with 14-3-3 proteins Tyrosine hydroxylase was the the first recognized binding partner of 14-3-3 proteins [73]. The 14-3-3 interacts with RD-hTHl phosphorylated at either or both serines S19 and S40, while no interaction is observed with pS31 [74]. As there are two binding sites within the 14-3-3 dimer and up to four phosphoserines within the RD-hTHl dimer, the interaction behaviour is quite complex. A simplified version of this interaction, avoiding the RD-hTHl dimer by using an N-terminal 50-residue peptide, was studied using3 1 P NMR [52]. The results show pS19 having higher affinity towards 14-3-3£ than pS40 by a factor of 10, while three main complexes were observed depending on conditions - one 14-3-3 subunit was always occupied by pS19, while the other was binding either pS40 from the same RD-hTHl, or pS19 from a second peptide, or was free. 1.4 Intrinsically disordered proteins Many binding partners of 14-3-3 are either fully disordered, or at least their binding sites are located within disordered regions, including hTHl. Intrinsically disordered proteins (IDPs) in general are characterized by a lack of a stable and well defined tertiary structure. It is suggested that a third up to a half of all proteins in eukaryotic genomes contain long disordered regions [75-77], while the percentage rises to more than 70 % for signalling proteins [78]. The IDPs are playing key role in many cellular processes such as signal transduction, molecular recognition, protein phosphorylation, or regulation of transcription and translation [79]. Serious neurodegenerative diseases are directly linked with IDPs, including proteins tau and a-synuclein involved in Alzheimer's disease [80]. The main distinctive feature of IDPs is the disorder that can be understood in terms of potential energy surface. Unlike rigid ordered proteins that inhabit a distinct deep global minimum on their energy surface, the surface of IDPs resembles more a landscape full of shallow valleys and depressions along low ridges and hills [81], allowing many different conformations with similar stability. The resulting „structure" can be best studied in terms of an ensemble of 10 Introduction 1.4 Intrinsically disordered proteins dynamically changing individual structures. This flexibility allows the IDPs to easily interact with multiple binding partners with simultaneously high specificity and low affinity [75, 82, 83] and thus serve as important interaction hubs or signal transmitters. There are significant experimental and computational difficulties when it comes to IDPs characterization. The situation is even worse when the object of study contains both ordered and disordered regions, as exemplified by the regulatory domain of tyrosine hydroxylase (see above). Many traditional approaches of structural biology fail in such a scenario - X-ray crystallography suffers from the presence of disordered part hindering the crystallization or forcing the disordered part to fold into non-native state; cryoelectron microscopy averaging images is doomed to observe only the ordered part as the IDP part averages to noise; computational atomic simulations become prohibitively slow and expensive due to the need of a large simulation box with a huge number of solvent molecules. The difficulties of describing IDPs can be overcome by the usage of nuclear magnetic resonance (NMR) employing high-resolution high-dimensional (usually 5-7D) experimental procedures followed by automated data processing for assignment and evaluation of relaxation properties and residual dipolar couplings (RDCs) [84- 87]. The RDCs provide information regarding the compactness of the structure, while the relaxation properties indicate the amount of internal movement and flexibility. While these approaches work well for the disordered regions, the ordered domains exhibit much faster relaxation due to their rigidity, leading to lower signals and difficulty to employ high-dimensional (>3D) spectroscopy. These effects make it very difficult to design and run experiments that would be able to observe all residues in these hybrid proteins. 11 Introduction 1.5 Interaction modelling 1.5 Interaction modelling This section aims to describe the approach to model general molecular interactions from the viewpoint of thermodynamics and kinet- ics. In the simplest case of two molecules A, B interaction leading to AB complex, the system composition in equilibrium can be described by the following equations. A + B <==± AB (1.1) k0ff „ _ fcofl _ [A] [B] K i ~ K» ~ W ( ' where [A], [B], and [AB] denote the equilibrium concentrations of the separate molecules and their complex, kon and feQff denote the rate constants of association and dissociation of the complex, and denotes the dissociation equilibrium constant of the complex. The temporal development of the system is described by a set of differential equations ^ [A] = -kon [A] [B] + koS [AB] (1.3) ^ [B] = - A w [A] [B] + koS [AB] (1.4) ^ [ A B ] = fcon [A] [B] - fcoff [AB], (1.5) where [A], [B], and [AB] denote the concentrations of the separate molecules and their complex, while the left-hand site terms like 4T [A] denote the temporal change of individual concentrations. The boundary conditions are given by the laws of mass conser- vation. [A] + 2 [AA] + [AB] = Atot (1.6) [B] +2[BB] + [AB]=Bt o t, (1-7) where At o t and Bt o t represent the total concentration of monomeric units of A and B, that can be readily measured using classical biochemical equipment, for instance by UV-Vis absorption spectropho- tometry. 12 Introduction 1.5 Interaction modelling Both rate constants kon,k0s and equilibrium dissociation constant Kd are temperature dependent. This dependency can be described by the following formulas that allow to determine the values of transient (denoted by and equilibrium (denoted by •&) Gibbs free energy A C , enthalpy AH, and entropy A S . /coff = — e x p j ^ - — J (1.8) KkBT f AH*-TAS*\ k o S = — 6 X P V BT ) ( L 9 ) , , AH* 1 AS* , nkBT * * « = — -T + - R - + I a - t ( L 1 0 ) KA AG* = -RT In — =AH*-TAS* (1.11) 1M v ' Kd AH* 1 A S * h , m = ~ B - f+ - B - - ( L 1 2 ) where K in the Eyring equation (Eq. 1.8) [88] denotes the transmission coefficient that is usually assumed to be equal to one, kB is the Boltzmann constant and h is the Planck constant. All thermodynamic parameters here and in the whole thesis are considered in the direction of complex dissociation, i.e. reversed from the traditional convention. 1.5.1 Homo- and heterodimers Suppose we study the a dimeric protein A forming dimers A A with dissociation constant and dissociation rate constant £;0ff (Eq. 1.13). To track the dimerization, we typically label some amount of the monomeric units A by a label that does not interfere with the dimeric interface. The labelled protein will be denoted B for distinguishing it from the unlabelled protein A. As the label does not affect the dimerization, the dissociation properties of 13 Introduction 1.5 Interaction modelling dimer B B is identical to A A . 2 A ^ A A , ^ = ^ = g (1.13) 2B <==± B B , K D ) B B = ^ = |L= K d A A (L14) k0ff Kon L-D-Dj A . o ^ ^ . o ~ _ fcoff _ [A] [B] _ K d , A A A + B ^ T A B ' ' A B ~~ 2Aw ~~ ~[ABT ~~ 2 ( L 1 5 ) The association leading to complex A B happens with a rate constant that is twice larger than the rate constants feon for the other reactions. It follows from this that also the dissociation constant for the A B complex is half the dissociation for other complexes. One of the consequences is the composition of the final solution that contains (excluding monomers) dimers A A , A B , and B B in ratio 1 : 2 : 1 . The reason for the factor of 2 stems from the fact that the label allows to distinguish the subunits in the dimer, therefore there are two distinguishable ways to form the dimer (either A B or B A ) , while there is only one way to form the dimer A A . Note that if the label would interfere with the interaction, each of the three reactions (Eqs. 1.13-1.15) would have a separate, unrelated rate constants and dissociation constant and the system would get much more difficult to analyse due to having more degrees of free- dom. The change of concentrations of individual molecules in time can be described by the following set of differential equations. d dt d [A] = -2feo n ( [ A ] 2 + [A] [B]) + feoff (2 [AA] + [AB]) (1.16) ^ [B] = -2feo n ( [ B ] 2 + [A] [B]) + feoff (2 [BB] + [AB]) (1.17) ^ [AA] = feon [ A ] 2 - feoff [AA] (1.18) ^ [BB] = feon [B]2 - feoff [BB] (1.19) ^ [AB] = 2feon [A] [B] - feoff [AB] (1.20) The laws of mass conservation are identical to those described by Eqs. 1.6 and 1.7. 14 Introduction 1.5 Interaction modelling The initial conditions are an interesting topic to study closer. Usually, when performing an interaction experiment, one starts with pure protein samples with known concentrations determined using e.g. UV-Vis absorption spectrophotometry. The values of concentrations then represent the amount of monomeric units in a given volume (At o t , B t o t ) . The initial concentrations of monomers and dimers can be derived from the equilibrium equation (Eq. 1.13) and law of mass conservation ([A] +2 [AA] = A t o t ) in the following form. = yjK'j + 8 K d A t o t - Kd [ J 4 rAAl = A t o t v ^ d + S^dAtot - Kd [ J 2 8 The other protein is usually added to the solution of the first protein in a way that minimizes the change in volume in the experiment chamber. Therefore, a small amount (Vjnit) of a sample with high concentration Btot,init is added to achieve the final volume of V and final concentration Btot- However, the monomerdimer equilibrium at the high concentration differs (there are more dimers) from the equilibrium that would be adequate at the lower concentration in the final volume. If a kinetic measurement is precise enough, there is a noticeable lag in the interaction due to the time necessary for equilibration of protein B after its concentration jump. A detailed numerical simulation of the kinetic process must take this effect in consideration. The concentration of monomers after the dilution are described by the following equations. The dimer concentration can be easily calculated using the law of mass conservation. (1.21) (1.22) 15 Introduction 1.5 Interaction modelling V Btot,init = Btot • T~R— (1-23) "init [B]i n i t \JK\ + 8KdBtot,init - Kd = \ i^/KJ + 8Kd • B t o t • ^ - Kd ] (1.24) [B] = [B]ir V h l i t Jinit ( K d ^ ) + 8 ^ d - B t o f ^ - ^ d ^ v (1.25) where [B]i n i t denotes the concentration of monomers in the original high-concentration sample and [B] denotes the concentration of B monomers directly after the mixing. 16 Chapter 2 Aims of the Thesis The main focus of the thesis is to characterize the dimerization process of 14-3-3^ protein and the impact of serine 58 phosphorylation upon it. The particular aims toward this goal include developing a robust design and analytical methodology of several complementary fluorescence assays and using this methodology to determine the dissociation equilibrium and rate constants of the protein. Second, using the results to calculate the thermodynamic parameters of the dimerization. Third, measure the effect of S58 phosphorylation on the dimerization, both in case of homodimerization of phosphorylated proteins and heterodimerization between phosphorylated and wildtype variant. The aims of the other projects were to characterize the phosphorylation of the regulatory domain of tyrosine hydroxylase, and to measure its NMR relaxation properties to benchmark force fields used in molecular dynamics. • Biophysical characterization of dimerization properties of 14- 3-3C• Quantitative description of the effect of S58 phosphorylation upon dimerization of 14-3-3C• Description of disordered region of regulatory domain of tyrosine hydroxylase phosphorylation and relaxation properties. 17 Chapter 3 Methods This chapter describes the methods and analyses employed in the projects which results are presented in this thesis. The chapter focuses mainly on the experiments and methodologies I ran and used personally, while briefly describing methods used for sample preparation and others. 3.1 Fluorescence assays 3.1.1 Sample preparation All used 14-3-3 constructs were purified as published previously [52]. Here, I briefly describe the most important steps in expression, purification and labelling of the proteins. The sample preparation was performed mostly by my colleague Zuzana Trosanova. All 14-3-3 variants, containing at the N-terminus His-tag followed by the TEV cleavage site, were produced in BL21(DE3)RIL cells, which were subsequently homogenized and centrifuged. The supernatant was purified using a Ni affinity column followed by a gel filtration column. After that, the His-tag was cleaved using TEV protease and subjected to a second nickel affinity column, the flow-through was collected in this case. After reduction by TCEP, the next purification step comprised an anion exchange column, where the pure protein eluted in a single peak during elution buffer gradient (changing the concentration of NaCl). Finally, the sample was dialyzed into 20 mM phosphate buffer, pH 8.0. The purity of 18 Methods 3.1 Fluorescence assays the final 14-3-3 protein sample was verified by MALDI-TOF mass spectrometry. In order to specifically label 14-3-3 protein using fluorescent dyes, we designed a construct with inserted SVDACKGSSGG linker sequence, containing a single cysteine, preceding the N-terminus of the wildtype. To ensure specific labelling, two surface accessible cysteines C25 and C189 were mutated out to alanines. Three dyes were used in total - 6-tetramethyl-rhodamine-C6-maleimide (TMR) for self-quenching (SQ) experiments, AlexaFluor488-C2maleimide (AF488) as a Förster resonance energy transfer (FRET) donor, AlexaFluor647-C5-maleimide (AF647) as a FRET acceptor and for microscale thermophoresis (MST) detection. The optimized labelling protocol contains the following steps. Firstly, the protein in phosphate buffer (pH 6.8) is treated with TCEP to reduce all disulfide bonds, while not preventing the reaction between cysteine in the linker and the maleimide groups of the dyes. Then, the fluorescent dye in dimethyl sulfoxide is added to the protein at 30:1 molar ratio. After 45 minutes incubation in dark, the unreacted dye was removed using a 30kDa molecular sieve. The final sample quality was confirmed by MALDI-TOF mass spectrometry (to check the extent of labelling) and SDS-PAGE (to check the excess dye removal). 3.1.2 Experimental setup Fluorescence measurements were conducted mainly by Zuzana Trosanovä and Aneta Kozelekovä employing a FluoroMax-4 Spectrofluorometer (HORIBA Jobin Yvon). Before each measurement, the quartz cuvette was treated with 10mg/mL bovine serum albumine (BSA) for 30 min to avoid adhesion of the fluorescent dyes to the walls. For the SQ assay, A e x — 553 nm (bandwidth of 0.8 nm) and Aem = 575 nm (bandwidth 2.5 nm) were used. For the FRET assay, Ae x = 470 nm (bandwidth 4.2 nm) and A e m = 666 nm (bandwidth 10.5 nm) were used. Unless specified otherwise, the measurement was performed at 37 °C in 20 mM NaPi buffer, pH = 6.8. Microscale thermophoresis (MST) assay was performed by Zuzana Trosanovä at constant concentration of labelled 14-3-3(_AF647 (c = 0.5 nM) and varying concentration of unlabelled 14-3-3^ (c = 6pM-100nM) in a buffer containing 20 mM phosphate, 0.5mg/mL BSA and 0.05% TWEEN-20 at 30 °C. Binding studies were per- 19 Methods 3.1 Fluorescence assays formed in triplicates using 50% LED and 80% MST power with a Monolith NT.115Pico device (NanoTemper Technologies) using standard capillaries. 3.1.3 Model description The experimental data from fluorescence and MST assays were analysed by the author using software COPASI v4.24 [89]. Several different models were constructed, one for each assay. Each model comprises up to 21 chemical reactions defining relations between the various species present in the sample. The chemical equations, as well as differential equations for time dependency of concentrations, follow the general examples shown in section 1.5. Microscale Thermophoresis model The MST model, used for analysis of microscale thermophoresis experiments, describes the equilibria between nuorescently (by AF647) labelled (L) and non-labelled (N) 14-3-3 proteins and calculates the resulting normalized fluorescence intensity. The MST measures the system in an equilibrated state, therefore the full description of the system composition is achieved using the following static equilibrium equations. The differential equations, describing the changes of concentrations in time, required in other models below are not necessary here. 2 L ^ L L , Kd= [iM <3 '2) N + L ^ N L , f . M H ( , 3 ) F -F H +F [ N L ] + 2 [ L L ] f34l -Tnorm — *M " J r -TD - 7 j where FM and F D are the normalized fluorescence -Fnorm values corresponding to pure monomeric and dimeric protein. 20 Methods 3.1 Fluorescence assays Self-quenching model The following equations describe the model used for analysis of the self-quenching fluorescence assays employing a system of TMRlabelled (Q) and unlabelled (N) proteins. 2 N ^ N N , Kd = LL (3.5) 2 Q ^ Q Q , ifd = 4 (3-6) k0 ff ^ [N] = -2/co n ([N]2 + [N] [Q]) + koS (2 [NN] + [NQ]) (3.8) ^ [Q] = -2kon ([Q]2 + [N] [Q]) + koS (2 [QQ] + [NQ]) (3.9) [NN] = kon [N]2 - k0Q [NN] (3.10) ^ [QQ] = kon [Q]2 - koS [QQ] (3.11) ^ [NQ] = 2kon [N] [Q] - koS [NQ] (3.12) The measured fluorescence intensity is then calculated as follows. The raw value is used directly in the analysis of kinetic assay, while the intensity normalized by the total concentration of labelled protein is utilized for the equilibrium assay. hot = Ion • ([Q] + [NQ]) + /0 ff • 2 [QQ] (3.13) j _ T _[Q]_ , j 2[QQ] -'SO.norm 'on ^ i 'off ^ 'SQ,norm — 'on - 1" 'off ~7\ j (3-14) where the values of Ion and Iqq denote the normalized fluorescence intensities of non-quenched labelled protein and quenched labelled protein. FRET model The following equations describe the model used for analysis of the FRET fluorescence assays employing a system of 14-3-3 proteins 21 Methods 3.1 Fluorescence assays which can be unlabelled (N), labelled by the FRET donor AF488 (D), or labelled by the FRET acceptor AF647 (A). (1 d7 koff at ^ N N , koff dt ND, koff dt 2kon d z=± NA, — k0ff dt 2k0„( _ d = DA. koff dt D + D <=> DD A + A f i ' ' A A N + N "± N :\ N + D «= A XI) N + A<= A D + A «= A I> A [DD] = fcon[D]2 -feoff[DD] (3.15) [AA] = fe0n [A] -feoff[AA] (3.16) [NN] = /Con [N]2 -feoff[NN] (3.17) [ND] = 2kon [N] [D] - feoff [ND] (3.18) [NA] = 2feon [N] [A] -feoff[NA] (3.19) [DA] = 2feon [D] [A] -feoff[DA], (3.20) ^ [ A ] = -2feo n ([A]2 + [D] [A] + [N] [Af +feoff(2 [AA] + [DA] + [NA]) (3.21) ^ [D] = -2feo n ([D]2 + [D] [A] + [N] [D]) +feoff(2 [DD] + [DA] + [ND]) (3.22) ^[N] = -2feo n ([N]2 + [N] [A] + [N] [D]) +feoff(2 [NN] + [NA] + [ND]) (3.23) The observed signal J0bs is a linear combination of concentrations of fluorescently labelled proteins. -^obs = 7f r c t • [DA] + IA • ([A] + 2 [AA] + [NA]) + / D • ([D] + 2 [DD] + [ND]) (3.24) [DA] [A] + 2 [AA] JFRET.norm — Jfret "T T~p; r i a "T —p; H A tot + -Utot At o t + -Utot [D] + 2 [DD] + i d ' I I n ' (3 -2 5 ) A tot "T -Utot where /FRET,norm is the normalized fluorescence intensity of FRET pair (DA), IA and are the normalizedfluorescenceintensities of acceptor and donor labelled monomer.22 Methods 3.1 Fluorescence assays F R E T heterodimerization model As the heterodimerization assay was performed mostly in the equilibrium mode, the equations are given in the equilibrium form using dissociation constants for simplicity. However, the implementation of the model in COPASI is analogical to the former models, using rate constants as well. The presented equations describe such an assay setup, where acceptor (A) and donor (D) labelled 14-3-3£ and fluorescently nonlabelled S58-phosphorylated (p, or p£) 14-3-3£ proteins interact. The model is readily usable also for combinations of other isoforms, where the only difference is renaming of the entities. In order to model other possible combinations, the general model contains 21 analogical interactions of every possible pair of proteins with either label (or none) and possibly phosphorylated, described by three sets of equilibrium and rate constants. The full model is not shown here for the sake of brevity. The presented equations describe an experiment, where firstly acceptor- and donor-labelled non-phosphorylated proteins interact, followed by addition of nonlabelled phosphorylated protein, exemplified at Fig. 4.5. i ' 2 [AD] K d,e/c [Al2 A + A A A ' KdX/c = Kd,c/c [D]2 D + D ( >DD, KdX/c 1 1 [DD] Kd ,P /c t^ [p] [A] P + A \ ' PA > K t,P/c = ^^X] Kd,P/c [pi [Dl P + D ( >pD, Kd,p/C - U l [ 1 [PD] K d ' P / P , „„ [P]2 P + P < > PP, K d,p/p [PP] [A] + 2 [AA] + [AD] + [pA] = A t o t [D] + 2 [DD] + [AD] + [pD] = D t o t [p] + 2 [pp] + [pA] + [pD] = p t o t , 3.26) 3.27) 3.28) 3.29) 3.30) 3.31) 3.32) 3.33) 3.34) 23 Methods 3.1 Fluorescence assays where the symbols AD, pA and other correspond to dimers (_ A / ( _ A, p£/£_A, and analogous. The observed fluorescence intensity /FRET, norm normalized to the total concentration of nuorescently active molecules, is calculated as follows: ^FRET.norm = -^fret + / d 3.1.4 Error analysis The confidence intervals were estimated using a method described by Johnson et al [90], by evaluation of sum of squared errors (SSE) in dependency on Kd and/or fc0ffTo obtain the SSE contour plot, the parameter(s) were systematically varied over a range close to the best fit, while adjusting the other parameters (i.e. fluorescence scaling factors). The variation is achieved using Parameter Scan task within COPASI. For each combination of parameters rc, y, after the best possible adjustment was achieved, the SSE^y was calculated. The upper and lower limits of parameters were obtained at a threshold of 0.8 at the reciprocal normalized SSE plot (calculated as SSEmin/SSE^y). 3.1.5 Detailed description of COPASI workflow For the analysis of experimental data from the fluorescence and MST assays, the COPASI (COmplex PAthway Simulator) software [89] was used. This section describes in detail how to perform the curve fitting and error analysis for the experiments. The most general and complex model prepared is named f ret_kinetics.eps. This model allows to simulate and/or evaluate any experiment based on FRET assay, which can contain two different protein species (either different isoforms or a single isoform in phosphorylated and non-phosphorylated state), arbitrarily labelled by fluorescent dyes. From this general model, three specialized models are derived - f ret_homo_kinetics . cps for analysis of homodimerization kinetic assays, fret_dilutions.eps for analysis of general equilibrium assays, and fret_homo_dilutions.eps for analysis of homodimerization equilibrium assays. [AD] A-tot + Dtot [D] + 2 [DP] + [pD] Atot + D _ [A] + 2 [AA] + [pA] Atot + Dt-'tot •'tot (3.35) 24 Methods 3.1 Fluorescence assays The models for analysis of self-quenching assays are named sq_kinetics . cps and sq_dilutions. cps. Both these models are prepared also for heterodimeric experiments. The system of equation and general workflow is analogous to the FRET models. The main difference lies in calculation of the resulting signal. The last two models, mst_homo.eps and mst_hetero.eps are designed for MST experiments either in homo-, or heterodimerization setup. Again, the equations describing the concentrations and equilibria are analogous to the FRET models, the difference lies in the evaluation of measured signal. Firstly, the workflow for the analysis of the most complex heterodimerization FRET kinetic assay will be presented. Other kinetic assays are similar and can be easily derived. Next, the workflow for homodimerization SQ equilibrium assay will be described as an example of the equilibrium model type. The MST assay analysis is similar to this one and will be mentioned only briefly F R E T kinetic assay The workflow from raw data from fluorimeter to fitted parameters and errors consists of several steps - experimental details definition, best parameters fitting, and parameter scanning to get the error analysis. The course of experiment is divided into several stages separated by addition of various samples, for an example of an experiment see Figure 4.2C on page 43. The stages are defined using Events located within COPASI - Model - Biochemical, see Figure 3.1. The default events are addition of acceptor-labelled protein aX, addition of donor-labelled protein dY, and addition of non-labelled protein Y. There are two more stages defined (additions of another aY and Y) that can be used in specific experimental setups. The definition of events is fully extensible and can support arbitrarily complex experiments. The Events are all similar to each other and their definition differs only in the species involved. The specific details for Event N is governed by Global Quantities in the form of add_N_*. The Event 1 is triggered at time add_l_time (which might be equal to 0). At this time, the sample of aY is added to the measuring cuvette with original volume of V^uvette- The sample is usually added as a small volume add_l_volume of highly concentrated stock with 25 Methods 3.1 Fluorescence assays (a fret kinetics COPASI 4.30 (Build 240) C:/'Us. File Edit lools Window Help u a u g * v ss s s en a s/../fluomer/fret kinetics.eps COPASI ^ Model v Biochemical • Compartments [II Species [33] Reactions [21] Global Quantities [36] Events [5] Parameter Overview Parameter Sets [29] Mathematical Diagrams Tasks Output Specifications Functions [40] Units [35] Concentrations Search: [ # 1 Name 01 add phospho-acceptor Trigger Expression Time >= Values[add_l_time] Delayed No Delay Expression Assignment Target Compartments[cuvett e) c_dir c m i Com 2 02 add non-phospho donor Time >Values[add_2_tlme] NO dvdv dY Compartments[cuvette] c_dii c_m- Con 3 03 add non-labeled non-phospho Time >Values[add_3_time] No Y YY Compartments[cuvett e) c~dii Com • Madd non-phospho acceptor Time >=values[add_4_timel No Compartments[cuvetteI aYaY •Y Con c mi 5 05 add non-labeled non-phospho Time >Values[add_5_time] No Y YY Co mpartments[ cuvette) c mi c_dii5 New Event I Mew I ßelete Figure 3.1: Screenshot of COPASI showing the Events dialog win- dow. concentration add_l_conc_stock to reach the final concentration of add_l_conc_f inal. All these concentrations are specified as „total" concentrations of monomeric units. The actual concentrations of monomers and dimers are calculated so that their ratio corresponds to the ratio in equilibrated stock solution, i.e. with higher concentration of dimers than what should be present in the final (more diluted) solution. This concentration Jump" after which the sample equilibrates in the new conditions explains a lag in signal increase, sometimes observable after the addition of donor-labelled sample to the acceptor-labelled one already equilibrated in the cu- vette. The details of each analysed experiment are put to the Global Quantities window, see Figure 3.2. The values of add_N_time, add_N_conc_f i n a l and add_N_conc_stock must be filled with time of addition, final cuvette concentration and pipetted stock concentration of each sample; the added volumes (in litres) are calculated automatically, however it is recommended to double check those as well. Also recommendable is to check the definitions of Events, es- 26 Methods 3.1 Fluorescence assays & fre: kinetics - COPASI 4.30 (Build 24C) C:/Users/.../fluor .File Edit Tools window Help ner/fret kinetics.eps u ; Ü JH y 9 * * s i i CB .*J Concentrations " v COPASI v Model v Biochemical Compartments [1) > Species [33] Reactions [21] Global Quantities [36] v COPASI v Model v Biochemical Compartments [1) > Species [33] Reactions [21] Global Quantities [36] Search: v COPASI v Model v Biochemical Compartments [1) > Species [33] Reactions [21] Global Quantities [36] 1 Name Kd_XX Type fixed Unit nmol/l Initial Value [Unit] 5e+06 Transient Value Rate [Unit] [Unit/min] nan 0 Initial Expression [Unit] v COPASI v Model v Biochemical Compartments [1) > Species [33] Reactions [21] Global Quantities [36] 2 KdYY fixed nmol/l 5 nan 0 Events [5] Parameter Overview Parameter Sets [29] Mathematical v Tasks > Steady-State > Stoichiometric Analysis > Time Course Metabolic Control Analysis > Lyapunov Exponents Time Scale Separation Analysis • Cross Section Parameter Scan > Optimization v Parameter Estimation Result > Sensitivities 3 KdXY fixed nmol/l 2609.95 nan 0Events [5] Parameter Overview Parameter Sets [29] Mathematical v Tasks > Steady-State > Stoichiometric Analysis > Time Course Metabolic Control Analysis > Lyapunov Exponents Time Scale Separation Analysis • Cross Section Parameter Scan > Optimization v Parameter Estimation Result > Sensitivities 4 k_off_XX fixed l/min 0.1 nan 0 Events [5] Parameter Overview Parameter Sets [29] Mathematical v Tasks > Steady-State > Stoichiometric Analysis > Time Course Metabolic Control Analysis > Lyapunov Exponents Time Scale Separation Analysis • Cross Section Parameter Scan > Optimization v Parameter Estimation Result > Sensitivities 5 k_off_yv fixed l/min 0.156 nan 0 Events [5] Parameter Overview Parameter Sets [29] Mathematical v Tasks > Steady-State > Stoichiometric Analysis > Time Course Metabolic Control Analysis > Lyapunov Exponents Time Scale Separation Analysis • Cross Section Parameter Scan > Optimization v Parameter Estimation Result > Sensitivities 6 k_off_XY fixed l/min Steady-State > Stoichiometric Analysis > Time Course Metabolic Control Analysis > Lyapunov Exponents Time Scale Separation Analysis • Cross Section Parameter Scan > Optimization v Parameter Estimation Result > Sensitivities 7 k_on_XX assignment l/hmol*mm) 2e-08 nan nan Values[k_off_ Events [5] Parameter Overview Parameter Sets [29] Mathematical v Tasks > Steady-State > Stoichiometric Analysis > Time Course Metabolic Control Analysis > Lyapunov Exponents Time Scale Separation Analysis • Cross Section Parameter Scan > Optimization v Parameter Estimation Result > Sensitivities 8 k_on_W assignment l/(nmol*min) 0.0312 nan nan Values[k_off_ Events [5] Parameter Overview Parameter Sets [29] Mathematical v Tasks > Steady-State > Stoichiometric Analysis > Time Course Metabolic Control Analysis > Lyapunov Exponents Time Scale Separation Analysis • Cross Section Parameter Scan > Optimization v Parameter Estimation Result > Sensitivities 9 k_on_XY assignment l/(nmol*min) 7.66297e-05 nan nan Values[k_off_ Events [5] Parameter Overview Parameter Sets [29] Mathematical v Tasks > Steady-State > Stoichiometric Analysis > Time Course Metabolic Control Analysis > Lyapunov Exponents Time Scale Separation Analysis • Cross Section Parameter Scan > Optimization v Parameter Estimation Result > Sensitivities 10 double_k_on_XX assignment l/(nmol*min) 4e-08 nan nan 2 'Values [k_oi Events [5] Parameter Overview Parameter Sets [29] Mathematical v Tasks > Steady-State > Stoichiometric Analysis > Time Course Metabolic Control Analysis > Lyapunov Exponents Time Scale Separation Analysis • Cross Section Parameter Scan > Optimization v Parameter Estimation Result > Sensitivities 11 double_k_on_W assignment l/(nmol*min) 0.0624 nan nan 2*Values[k_oi Events [5] Parameter Overview Parameter Sets [29] Mathematical v Tasks > Steady-State > Stoichiometric Analysis > Time Course Metabolic Control Analysis > Lyapunov Exponents Time Scale Separation Analysis • Cross Section Parameter Scan > Optimization v Parameter Estimation Result > Sensitivities 12 lnten5ity_donor fixed l/nmol 2300 nan 0 Events [5] Parameter Overview Parameter Sets [29] Mathematical v Tasks > Steady-State > Stoichiometric Analysis > Time Course Metabolic Control Analysis > Lyapunov Exponents Time Scale Separation Analysis • Cross Section Parameter Scan > Optimization v Parameter Estimation Result > Sensitivities 13 mtensity_acceptor fixed l/nmol 70 nan 0 Events [5] Parameter Overview Parameter Sets [29] Mathematical v Tasks > Steady-State > Stoichiometric Analysis > Time Course Metabolic Control Analysis > Lyapunov Exponents Time Scale Separation Analysis • Cross Section Parameter Scan > Optimization v Parameter Estimation Result > Sensitivities 14 IntensityFRET fixed l/nmol 4000 nan 0 Events [5] Parameter Overview Parameter Sets [29] Mathematical v Tasks > Steady-State > Stoichiometric Analysis > Time Course Metabolic Control Analysis > Lyapunov Exponents Time Scale Separation Analysis • Cross Section Parameter Scan > Optimization v Parameter Estimation Result > Sensitivities IS lntensity_total assignment 0 nan nan Values[lntens Events [5] Parameter Overview Parameter Sets [29] Mathematical v Tasks > Steady-State > Stoichiometric Analysis > Time Course Metabolic Control Analysis > Lyapunov Exponents Time Scale Separation Analysis • Cross Section Parameter Scan > Optimization v Parameter Estimation Result > Sensitivities 16 add_l_time fixed min 14.2 nan O Linear Noise Approximation add_2_time fixed min 83.7 nan 0 Time Course Sensitivities > Output Specifications > Functions [40| > Units [35] 18 add_3_time fixed min 129.2 nan 0Time Course Sensitivities > Output Specifications > Functions [40| > Units [35] 19 add_4_time fixed min 1000 nan 0 Time Course Sensitivities > Output Specifications > Functions [40| > Units [35] 20 add_5_time fixed min 1000 nan 0 Time Course Sensitivities > Output Specifications > Functions [40| > Units [35] 21 add_l_CQrc_finsl fixed nmol/l 2000 nan 0 22 < add 2 cone final fixed nmol/l 20 nan 0 > Mew Delete Delete All Figure 3.2: Screenshot of COPASI showing the Global Quantities dialog window. pecially the order of samples. The numbers in the Quantity names denote the affiliation to the Events, not necessarily the order of addition, although it is less prone to error to follow the order. Note, that if the first addition time is larger than zero, the experiment is started without any labelled sample and the measured signal consists of blank only - this is highly recommendable setup, if possible as it reduces the fitted parameter space (because the value of blank is well established). The parameter fitting is performed in the Parameter Estimation window within COPASI - Tasks. The input is the data from fluorimeter in the form of two-column text file (might be in CSV format) with time (in minutes) in the first column and measured fluorescence signal (in raw counts) in the second column. The data should not contain any blank lines. This file is loaded through the Experimental Data dialog, where the Experiment Type is set to Time Course and the file columns are defined as Time and dependent with Model Object specified as Values [Intensity_total]. The Separator is usually the default tabulator, however other char- 27 Methods 3.1 Fluorescence assays Experimental Data File H Q 0 Experiment - 2000nM-pS58ZETA-647_20nM-ZETA-4S8.txt Experiment Experiment Experiment Header 1 0 Copy Settings D from previous CI to next Experiment Type O Steady State ® Time Course Weight Method Mean Square First Row 1 Separator 0 to all following Last Row 321 0 » 0 Normalize Weights per Experiment Column Name time Type Time Model Object Weight signal intensity dependent Values[lntensity_total] (1) OK Reven Cancel Figure 3.3: Screenshot of COPASI showing the Experimental Data dialog from the Parameter Estimation window. acters (e.g. comma or semicolon) might be used. The Header field might be filled if the file contains column names. After input data loading, the parameters to fit are specified in the main Parameter Estimation dialog as Objects, see Figure 3.4. The parameters are selected from the COPASI dialog tree (accessible by clicking the icon next to Object field); it is possible to select more parameters at once. In the analysis of FRET heterodimeric kinetic assay, the usual fitted parameters comprise of Kd_XY, k_off_XY, Intensity_acceptor (signal intensity of a monomeric acceptor-labelled protein at concentration of 1 nM), Intensity_donor, Intensity_FRET (signal intensity of a dimer containing both acceptor and donor dye, at concentration of InM), and blank (signal intensity measured without any labelled protein). After selecting the parameters, it is advisable to fill the Start Value for each parameter as precisely as possible; if this value is too far from the real value, the fitting might get stuck or find an 28 Methods 3.1 Fluorescence assays unreasonable set of parameter values resulting in a poor fit. The Bounds are usually less important parameters to set, however it is wise to set especially the Lower Bound to at least a very small positive number (such as le-6) to prevent numerical instabilities. The full heterodimerization model is described by a number of independent parameters, including three rate constants and three dissociation constants, one for each type of interaction - homodimerization to form X X , homodimerization to form Y Y and heterodimerization to form X Y dimer, see Eqs. 3.35. However, the experimental data do not contain sufficient information to estimate all these parameters simultaneously. Therefore, it is highly advisable to fill the known values of homodimerization parameters into Global Quantities before the estimation of heterodimeric parameters. The homodimerization parameters should be known from previous experiments using simpler homodimeric setup. The fitting is started by hitting the Run button. The recommended Method is NL2SOL [91] - an algorithm for non-linear leastsquares optimization. This method is very fast and able to find the locally optimal solution. However, it tends to fail when the Start Values of some parameters are too far from the true values, then the algorithm finds some non-optimal solution which manifests as a poor fit, often characterized by parameter values hitting their bounds. In such a case, either the Start Values must be set better, or (when the true values are not even estimated yet) some more robust method must be chosen. The experiences of the thesis author show that the methods most likely to find good results and globally optimal parameters are Evolutionary Programming, Genetic Algorithm, and Particle Swarm. The results of fitting are accessible in the Result window at the Parameters tab. When a good fit is obtained, it might be a good idea to Create Parameter Sets to save the resulting values and (maybe even more importantly) the experiment definition (times, concentrations, etc). For that purpose, run method Current Solution Statistics. The prepared parameter set can be found under COPASI - Model - Biochemical - Parameter Sets with a timestamp and experiment name. The final step in the workflow is the error analysis. The goal is to fit the data with selected parameters (usually and fc0ff) fixed at several values, find the best objective function (sum of squared errors, SSE) value by varying other parameters and plot the results to a 2D chart (referred to as a heatmap, see Fig. 4.2B on 29 Methods 3.1 Fluorescence assays <& rret_kiiielits - COPASI 4.30 (Build 240) C:AJiers/../riuurrier/rrfcl_kJri File Edit Tools Window Help i J U O ^ ^ S ?s —•' -4 Concentrations COPASI v Model v Biochemical Compartments [1| Species [33] Reactions [21] Global Quantities [36] : Events [5] Parameter Overview ParameterSets [29] Mathematics Diagrams " Tasks Steady-State Stoichiometric Analysis T ° ICom Metabolic Control Analysis Lyapunov Exponents TiniF Sralp Separation Analysis Cross Section Parameter Scan Optimization Parameter Estimation Result Sensitivities Linear Noise Approximation Time Course Sensitivities Output Specifications Functions [40] Units [35] Parameter Estimation • Randomize Start Values Ö Create Parameter Sets Parameters (3) Constraints (0) 0 Calculate Statistics D update model D executable Experimental Data Validation Data • Use Time Sens 1 0 < Values[lntensity_donor] .InitialVal je < 50000; Start Value = 2300 2 0 s values[intensity_FRET] .Initlalvalue £50000; Start value =15000 3 le-06 Species [33] Reactions [21] Global Quantities [36] Events [5] Parameter Overview » Parameter Sets [39] Mathematical Diagrams v Tasks Steady-State Stoichiometric Analysis Time Course Metabolic Control Analysis ( Lyapunov Exponents lime Scale Separation Analysis Cross Section Parameter Scan Optimization v Parameter Estimation Resuh Sensitivities • Linear Noise Approximation Time Course Sensitivities Output Specifications Functions [d0\ Units [35] Concentrations » Parameter Scan • • I! update model • executable New scan item: Scan Scan Intervals O Values Object Values[Kd_XY]. InitialValue Intervals min max 30 0.01 11000000 1 0 logarithmic Scan <9I Intervals O Values Object Values[k_off_XY].lnitialValue Intervals min max 30 0.01 1 1 • logarithmic T 3 S k Parameter Estimation * • Continue from Current State • output during subtask execution Q Continue an Error (&CQRepurlDe[iriiliuriSelei Report Template HEATMAP - Scan Parameters and Target function of parameter estimation Filename ^ho_kinetics/ou:puts/l-e3fn3pE/2000nM-pS58ZETA-647_20nM-ZETA-4SS_heatmap_large.txt | , j • Append 0 Confirm Overwrite I OK I Cancel I Bun Renort Qutput Assistant Figure 3.5: Screenshot of COPASI showing the Parameter Scan window with the Report dialog open. 32 Methods 3.2 NMR experiments lower than noise. The input data file is then constructed so that it contains two numerical columns - the protein concentration, and the mean signal observed for the corresponding concentration, calculated as an average from a plateau reached after the equilibration. Unlike in the kinetic assay, there is no setup to make, the Parameter Estimation can be launched directly. The input file is loaded through the Experimental Data window. The Experiment Type should be set to Steady State. The first column should be set as Independent with Object of [qX_tot] _0 (found in Species Initial Concentrations), in case of self-quenching. The second column should be set as Dependent Object of Values [I_total]. The fitted parameters are usually the dissociation constant Kd, the normalized fluorescence intensity of quenched and non-quenched labels I_off, I_on and possibly the signal value of pure buffer I_blank. To perform the error analysis, the Values [Kd] . InitialValue should be scanned, while removed from the Parameter Estimation. As only one parameter is scanned, the result can be plotted as a ID curve. The Microscale Thermophoresis assay is similar to the equilibrium assays described above. The input data are expected in the form of two-column text file with columns of unlabelled protein concentration and normalized fluorescence intensity -Fnorm expressed in thousandths (%o), for a reference data see Fig. 4.IE on page 40. The experiment analysis needs to start by setting the concentration of labelled protein [WT_label_tot], which is kept constant throughout the experiment. During data loading, the columns should be set to [WT_tot] _0 as the independent and MST_signal as the dependent variable. The fitted parameters are then the dissociation constant Kd and the signal intensities that would be observed for pure monomeric and dimeric proteins (l_m and I_d). The error analysis is performed with Kd as the scanned independent variable, resulting into a ID curve. 3.2 N M R experiments 3.2.1 Homodimerization of pS58_14-3-3< The 3 1 P spectra were measured over the course of 18 hours at 202.49 MHz using Bruker Avance NEO 500 MHz spectrometer equipped with a dual-band Prodigy cryoprobe. The pS58_14-3-3C sam- 33 Methods 3.2 NMR experiments pie concentration (i.e. concentration of monomeric subunits) was 1.1 mM. The peaks were fitted by Lorentz curves using a built-in deconvolution function in TopSpin software v4.0.8 (Bruker). The obtained peak areas (including confidence intervals) converted to concentrations of monomers and dimers based on the known total protein concentration. The concentrations then allowed to directly calculate the resulting together with its confidence intervals. 3.2.2 Regulatory domain of human tyrosine hydroxylase 1 The assignment of non-phosphorylated regulatory domain of hTHl (RD-hTHl) was performed using a Bruker 850 MHz US2 spectrometer equipped with cryogenic triple-resonance probe head. The kinetic measurements as well as assignment of the phosphorylated RD-hTHl were performed using a Bruker 600 MHz spectrometer equipped with a cryogenic triple-resonance probe head. Both probe heads were equipped with 2-axis gradient coils. All measurements were done at a temperature of 293.2 K. For the assignment of non-phosphorylated and doubly S19/S40 phosphorylated RD-hTHl, 3D HNCO [92], 5D HN(CA)CONH and 5D HabCabCONH [84, 93] were measured. All experiments were carried out with non-uniform sampling of the indirectly detected domains. The time schedule was generated using Poisson disk sampling on a grid, introducing distance constraints between points. The density of points was set according to a Gaussian distribution (a = 0.5). The 3D spectra were processed using a multidimensional Fourier transform [94], the 5D spectra were processed by a Sparse Multidimensional Fourier Transform [95]. The phosphorylation was monitored using time-resolved 2D HSQC experiments [96]. In both cases, the maximal evolution times were set to 102 and 64 ms, respectively. The time resolution was achieved using non-uniform oversampling in the indirect domain. The sampling schedule comprised 16,000 points in total for both phosphorylations. The size of the schedules exceeded the size of regular Nyquist grid 125 times. The total measuring time was 40 hours. The non-uniformly sampled (NUS) kinetic spectra were processed using coprocessed Multidimensional Decomposition approach [96]. The method slices the long 2D HSQC experiment into a se- 34 Methods 3.2 NMR experiments ries of time windows that are reorganized into a pseudo-3D spectrum with reaction time as the third dimension. The planes of this pseudo-3D spectrum differ only in the intensities of the peaks affected by the reaction, not in their locations. This invariance allows effective processing of the whole pseudo-3D spectrum at once, providing better resolution than achievable by a series of traditional 2D HSQC experiments. The method outputs a series of 2D spectra, where the intensities of the peaks can be analysed and rate constants can be calculated. The spectra were processed using nmrPipe [97], if not stated otherwise, and analysed by NMRFAM-Sparky [98]. 35 Chapter 4 Results and Discussion This chapter comprises the results of two main projects accompanied with short commentary on two additional projects I collaborated on. The first and largest project deals with the 14-3-3£ protein dimerization. The goal was to determine the kinetic and thermodynamic parameters of the dimerization. On top of this goal, we investigated the dimerization properties of 14-3-3C phosphorylated at serine S58. This work is summarized in Paper 1 with me as a shared first author. The methodology was successfully used also for isoform cr, the unpublished results are presented within this chapter. Our research group continues the efforts in this area. Using the models and analyses implemented by myself and described in this thesis, they have recently demonstrated that the often used approaches (e.g. phosphomimetics) to mimic the behaviour of phosphorylated S58 are not representative of the real dimerization properties, see [44]. In the second project, we worked with the regulatory domain of human tyrosine hydroxylase 1 (RD-hTHl) and its doubly phosphorylated variant. The goal was to perform the resonance assignment for all variants, investigate the kinetics of phosphorylation using non-uniform sampling methods and to analyse the effect of phosphorylation on the structural properties. The results of this work are presented in Paper 2, where I am the first author. The third project was focused on proteins containing both structured and intrinsically disordered regions, including RD-hTHl. This 36 Results and Discussion 4.1 14-3-3 dimerization work [99] was led by Vojtech Zapletal. It was focused on the validation of parameters of several different force-field for the prediction of structural ensembles of IDPs and disordered regions within hybrid proteins containing both disordered and rigid parts. While Vojtech performed the computational part, I contributed by measuring and analysing the relaxation NMR data serving as a reference for the atomistic simulations. 4.1 14-3-3 dimerization Although the dimeric nature of 14-3-3 proteins is well known, the important biophysical characteristics of the dimerization process remained unknown. Our study focused on the equilibrium dissociation constant and rate constants fc0ff and kon of the human isoform £. We studied extensively also the impact of phosphorylation at S58 on the dimerization properties. In less detail, the heterodimerization of a/C, pair was investigated. The 14-3-3^ homodimers are too stable for the monomerization to be observable using standard, so called traditional, biophysical techniques like native PAGE or analytical ultracentrifugation, that are sensitive enough for dimeric in micromolar range. This obstacle has been overcome by using highly sensitive fluorescence techniques employing Förster resonance energy transfer (FRET) and self-quenching (SQ) phenomena. These techniques, performed at a standard benchtop fluorimeter, were complemented with microscale thermophoresis (MST) measurements. The phosphorylation at S58 leads to a significant shift in the equilibrium from the dimeric to the monomeric form as the dissociation constant increases by many orders of magnitude. This very weak interaction between the phosphorylated monomers is prohibitively difficult to study using classical biophysical methods as the phosphorylated 14-3-3 exists almost exclusively in a monomeric state at the micromolar concentrations used by these methods. The problem was overcome by usage of P NMR spectroscopy that allowed to observe simultaneously both monomeric and dimeric form in the solution. My task was to design the experiments, develop a model for their description and use it to analyse the obtained data. Two sets of assays were designed to allow the determination 37 Results and Discussion 4.1 14-3-3 dimerization of thermodynamic (Kd) and kinetic (fc0ff) parameters of the 14-3- 3£ dimerization. All the presented assays use proteins labelled by fluorescent labels, namely 6-tetramethylrhodamine-C6-maleimide (TMR) for self-quenching assays, AlexaFluor647-C5-maleimide (AF647) serving as an acceptor in FRET assays and source of signal in MST assay, and AlexaFluor488-C2-maleimide (AF488) serving as a donor in FRET assays. 4.1.1 Equilibrium assays In the three equilibrium assays designed to determine the value of .Kd, the signal is measured in static conditions, when the system reaches equilibrium. The observed signal intensity then corresponds to the composition of the sample, allowing to calculate the dimerization dissociation constant. The general model describing the sample within equilibrium assays comprises several monomeric species and all possible dimeric pairs. The model then contains all possible equilibrium reactions between monomers leading to dimers, accompanied with corresponding dissociation constants. There are a few details requiring caution. Firstly, we assume that the dissociation constant is not affected by the fluorescent labels. This assumption is backed by the fact, that the value of K& extracted from all three presented assays using different labels and different approaches is very consistent and by the fact, that the labels are attached by a long flexible linker providing with enough freedom to not interfere with the dimeric interface. Secondly, the dissociation constant for interaction between the same proteins (e.g. non-phosphorylated 14-3-3£) bearing differing labels, is equal to half the dissociation constant of unlabelled dimer. This stems from the fact that the labelling breaks the symmetry present in the unlabelled dimer, leading to higher association rate constant by factor 2 as the heterodimer can be formed in two distinguishable ways but dissociated in a single way compared to the homodimer that is symmetric and can be both formed and dissociated in a single way In the MST assay, the AF647 labelled protein sample at a very low concentration (0.5 nM) is titrated with unlabelled protein. As the total £ concentration increases, dimers tend to dominate the sample composition. Therefore, at the start of the experiment, we observe the labelled protein mostly in monomeric 38 Results and Discussion 4.1 14-3-3 dimerization state, while towards the highest concentrations, we mostly observe dimeric proteins with one label attached. As the thermophoretic properties (e.g. charge, hydration shell) differ for the monomer and dimer, we observe a change of the normalized fluorescence intensity (Fig. 4.1g). Using the above described model (Eq. 3.4), we obtained the value of dissociation constant, = (3.2 ± 1.9) nM, see Fig. 4.1E,F. The analytical model provided by the MST vendor (NanoTemper Technologies) was not used as it assumes standard 1:1 receptor-ligand binding process, which is not applicable in this case where there are multiple interconnected equilibria. The self-quenching equilibrium assay make use of the fact that the intensity of TMR fluorescence decreases if two of its molecules come to close contact. The experiment start with a concentrated sample of ( _ T M R having a low fluorescence signal as there is a high percentage of dimers bearing two TMR labels quenching each other. The sample gets diluted to lower concentrations at each step. As a result of the dilution, the absolute fluorescence signal intensity decreases, therefore we normalize the signal with the total concentration of the protein. The normalized signal intensity increases with the decrease of protein concentration as the relative concentration of non-quenched monomers increases. The fitting of the experimental data using Eq. 3.14 provided the resulting value of dissociation constant of = (7.7 ± 1.1) nM, see Fig. 4.1A,B. The FRET assay is based on FRET effect - the donor (AF488) is excited and the energy is passed to the acceptor (AF647) if the two molecules are in close proximity, i.e. bound to subunits of one dimer. The assay setup is similar to the self-quenching assay - we start with concentrated sample containing both donorand acceptor-labelled 14-3-3 proteins. We proceed with dilution of the sample increasing the percentage of monomers, that are incapable of the FRET effect due to the presence of only one fluorescent dye molecule. Therefore, the relative fluorescence intensity decreases with decrease in the protein concentration. Evaluation of the measured data, according to Eq. 3.25 provides value of Kd = (4.6 ± 0.7) nM, see Fig. 4.1C,D. In summary, all three assays, while using differing underlying principles and different fluorescent dyes, provided very similar values of the dissociation constant. The final value, averaged with weights corresponding to the uncertainties of the individual assays, equals to (5.5 ± 0.8) nM. 39 Results and Discussion 4.1 14-3-3 dimerization B O in 1 10 100 ICQ] (nM) 0.01 0.1 1 10 100 K„ (nM) I- LU 1 10 100 [CA] + KJ>] (nM) 0.01 0.1 1 10 100 K„ (nM) I- 730 4 K_A] = 0.5 nM K d = 3 . 2 ± 1 . 9 n M I I I I L 0.001 0.01 0.1 1 10 100 H (nM) 0.01 0.1 1 10 100 M n M ) Figure 4.1: Experimental data and the error analysis of fluorescence equilibrium assays determining the homodimerization of 14-3-3^. Concentration dependence of fluorescence intensity within the SQ (A) and FRET (C) assays normalized to total protein concentration. Normalized MST fluorescence (E) of 0.5 nM C _ A titrated by unlabelled (. Error analyses of the SQ (B), FRET (D), and MST (F) assays. The solid lines show normalized sum of squares of errors, while the dashed lines show a threshold of 0.8. Adapted from Paper 1. 40 Results and Discussion 4.1 14-3-3 dimerization 4.1.2 Kinetic assays The two kinetic assays described in this section are designed to capture the dynamic nature of dimerization/monomerization process. For this task, time dependent fluorescent data were measured. The model describing these assays consists again of several differently labelled monomeric protein units combining into all possible dimers. While for the equilibrium assays, only dissociation constants were sufficient, for the kinetic assays, all rate constants have to be considered. This leads to a complex set of differential equations describing the development of concentration of all present species in time, see Section 3.1.3. This set of equations is integrated numerically using software COPASI [89]. Both the kinetic assays are designed as multi-phase experiments. The experiment is started with a cuvette filled with a buffer only (this phase is not shown in the figures), after several minutes of measuring the blank, a small volume of concentrated labelled (by TMR or AF647) protein in equimolar ratio is added and let to equilibrate. In case of FRET assay, a sample of donor (AF488) labelled protein is added, forming FRET pairs (heterodimer containing both donor and acceptor label) over the course of several dozen minutes. The last part starts by addition of an unlabelled protein in a high excess (typical molar ratio 100:1). This leads to disruption of either FRET pairs or self-quenched dimers by formation of "hetero"dimers containing only one fluorescently labelled subunit. The resulting value from self-quenching kinetic assay is fc0fr = (3.3 ± 0.7) x 10"3 s-1 . The FRET kinetic assay provided values of koS = (2.77 ± 0.09) x 10"3 s"1 and Kd = (3.6 ± 2.4) nM. The averaged value is (2.8 ± 0.4) x 1 0 _ 3 s _ 1 , corresponding to a mean dimer life-time of (6.0 ± 0.4) min. All these values were measured at37°C. The error analysis (described in detail in the next section) shows that the FRET kinetic assays allows to determine not only the value of kinetic rate constant fc0ff5 but also a good estimate of the dissociation constant This feature stands the FRET kinetic assay slightly above all the other assays as it provides more information than the self-quenching assay, while being faster and less sample consuming than the equilibrium assays. The possible downside of the need of three different protein samples (two labelled by dyes and one unlabelled) can be actually an advantage as it allows to 41 Results and Discussion 4.1 14-3-3 dimerization easily perform experiments using multiple different proteins simultaneously. Therefore, it is possible (as demonstrated below) to measure for instance the dissociation constant of a heterodimer between wildtype 14-3-3£ and its S58-phosphorylated variant. Other case might be measurement of heterodimerization of two 14-3-3 isoforms. The heterodimerization experiment requires labelling one protein with acceptor label, while the other with donor label, the unlabeled sample might be of whichever protein is more suitable. A l though in principle it should be possible to determine all three equilibrium and three rate constants via fitting, the variables are too correlated and numerous, therefore the fitting procedure converges slowly and the uncertainties tend to be large. However, the problem can be overcome by using good estimates of the homodimeric parameters, therefore fitting only the heterodimerization constants. The model for these experiments is a generalized and more complex form of the models used for homodimerization assay, taking into account all possible monomers, homo- and heterodimers. 4.1.3 Error analysis The uncertainty of resulting parameters can not be taken directly from the least-squares optimization method as the underlying model is quite complex and the parameters are often correlated. The problem of error analysis was overcome using the approach described by Johnson et al. [90]. The parameters of interest (usually or feoff) are varied across a range of values, while all other parameters are adjusted to achieve the best possible fit. For each (pair of) parameter value, we calculate the sum square error. The resulting values are then normalized by taking the reciprocal value as SSEm in /SSEa;) 2 / . In this way, the best fit is represented by value of 1, while the low numbers (close to 0) represent rather poor fit to the data. The threshold for error estimation was set to an SSE 25 % higher the best found SSE (corresponding to SSEm in /SSEa;) 2 / = 0.8). This value is somewhat arbitrary in the sense that it can not be rigorously derived from statistics, but it is rather a reasonable and conservative value that works for a broad range of models [90]. The described approach not only yields values of the parameter estimation uncertainties, but it provides with a visual represen- 42 Results and Discussion 4.1 14-3-3 dimerization A B Figure 4.2: The 14-3-3£ dimerization kinetic assays. SQ (A) and FRET (C) kinetic profiles at 37 °C. The dissociation rate constant fcoff, dimer mean lifetime r0fr and Kd (in case of FRET assay) are given. At each stage, a schematic representation of the most abundant species is included. The fit (red) is overlaid with experimental data (black). Dependence of sum square error (SSE) on dissociation rate constant fc0ff and dissociation constant Kd from the SQ (B) and FRET (D) assay. Minimum SSE (indicated by the symbol '+') divided by the SSE obtained at each rc, y coordinate (SSEm in /SSEa;) 2 / ) is presented. Adapted from Paper 1. 43 Results and Discussion 4.1 14-3-3 dimerization tation of the parameter estimation. An example of this analysis can be seen in Fig. 4.2. The confidence contour map presented for self-quenching assay (panel B) shows the case when one of the parameters (Kd in this case) is not well defined by the data. For every considered value below 100 nM, it is possible to choose other parameters in such that the resulting sum of squared errors is at most 25 % higher than for the best possible fit. This result means that while we can not determine the value of Kd, its upper limit was determined to 100 nM. The other parameter of interest, fc0ff is much more constrained and we can be confident that the real value will be in close proximity of the best fit. On the other hand, the FRET assay, while being more complex, bears more information and we are able to estimate not only k0^, but also with high confidence, see Fig. 4.2D). Very similarly, the errors can be estimated also for the remaining assays, where the analysis is simpler as only is evaluated, resulting in ID SSE profiles, see Fig. 4.1. 4.1.4 Effect of external conditions Using the described FRET and SQ assays, we examined multiple external factors that could potentially affect the 14-3-3 dimermonomer equilibria. The temperature and binding partner presence were the most significant factor affecting the 14-3-3^ dimerization. The change of pH had negligible and unsystematic influence (Fig. 4.3C,F). The analysis of temperature dependence is described in the section 4.1.7. The effect of buffer ionic strength was small (Fig. 4.3B,D) suggesting that hydrophobic interactions play a substantial role in the dimer stabilization. In case the ionic interactions were the main contributor to the dimerization, much stronger effect of the salt concentration would be expected. However, as we observe only a mild decrease in the fc0ff and Kd, compared to e.g. temperature or binding partner presence, we deduce the contribution of the polar contacts to be rather minor. The effect of binding partner presence is very strong and potentially very important for the role of 14-3-3 dimerization in vivo. As a prototypical 14-3-3 binding partner, we used the S19,S40phosphorylated regulatory domain of human tyrosine hydroxylase 1 (RD-hTHl). This domain is dimeric by itself containing two 44 Results and Discussion 4.1 14-3-3 dimerization B - • I I 1 0 1 1 OB 0 2 4 6 8 10 [RD-hTHl] (ekv.) PH 0 2 4 6 [RD-hTHl] (ekv.) 200 400 [NaCI] (mM) PH Figure 4.3: The effect of external conditions on the Q/Q dimerization. The presence of RD-hTHl lowers the rate constant feQff (A) and the dissociation constant (D) significantly. The increase in NaCI concentration (increasing ionic strength) slightly lowers the feoff (B), while having no systematic effect on the (E). The change in pH does not affect the feQff (C) nor the (F) significantly. The data were acquired using the FRET assay (empty circles), while the effects of NaCI and pH were studied also using the SQ assay (filled circles). Adapted from Paper 1. phosphoserines in each subunit, see Fig. 1.4. It was shown previously [52, 74] that the pS19 is the main interaction target for 14-3-3 proteins. The FRET kinetic assay was used to assess the impact of RDhTHl binding (Fig. 4.3A,D), using only the first part of the assay mixing acceptor and donor labelled proteins. The observed slowdown of the signal intensity increase was very prominent (more than 20 times smaller feQff) for RD-hTHl concentrations above 100 nM. The observed behaviour of dimer stabilization by a binding partner is highly relevant for the description of 14-3-3 in crowded cellular environment. In such conditions, surrounded by other proteins, the 14-3-3 dimers would be probably stabilized when bound to part- 45 Results and Discussion 4.1 14-3-3 dimerization ner proteins. However, as the data on phosphorylated 14-3-3 show, the monomers might be a highly relevant species especially following phosphorylation on S58. This evidence shows that while 14-3-3 acts as a hub in protein-protein interactions often relying on its dimeric nature, the protein itself is subject to several mechanisms that affect its dimeric properties. 4.1.5 Effect of phosphorylation at S58 The 14-3-3^ protein can be phosphorylated at S58 located at the dimeric interface. This modification strongly affects the dimeric properties, but the exact extent was previously unknown. Within this study, we measured both the homodimerization of the phosphorylated variant as well as the heterodimerization constants for its interaction with wildtype 14-3-3CFirstly, the pS58_14-3-3£ (denoted as p£ further on) homodimerization Kd was measured. We attempted to measure this constant using native PAGE experiment. However, the results (see Fig. 4.4A) were quite unconvincing as the used concentrations (up to 200 uM) were well below the K& and analysis of each lane provided inconsistent values in millimolar range. To measure more reliable data, we turned to NMR experiments. The system was particularly suitable for P NMR as there was exactly one phosphate present in each subunit. Moreover, the millimolar protein concentration appropriate for NMR was close to the estimated dissociation constant giving possibility to observe both monomeric and dimeric species. The sample of p£ (c = 1.1 mM) was measured at three distinct temperatures of 10, 20, and 37 °C, see Fig. 4.4. The peak assignment was achieved by measuring the sample at lower concentrations (down to 0.13 mM). As the protein concentration lowers, the relative proportion of monomeric form increases, as well as the relative volume of its peak in NMR spectrum, thus allowing to identify the positions of dimeric (about —2ppm) and monomeric signals (about 2ppm). Following the assignment, the peaks were fitted by Lorentz curves and integrated to obtain the areas under them, including the uncertainty stemming from the curve fitting. The concentrations of each form were derived from the areas and then employed to calculate dissociation constants for each temperature (Fig. 4.4C). The values of thermodynamic 46 Results and Discussion 4.1 14-3-3 dimerization c (MM) 10 dimer monomer 200 B 4 3 2 1 0 -1 -2 -3 -4 3 1 P chemical shift (ppm) 1000/T (K"1 ) Figure 4.4: Establishing the homodimerization K& of pS58_£. (A) Native PAGE gel containing non-phosphorylated Q at 10 uM concentration as reference, and p£ at concentrations of 1, 10, 20, 50, 75, 100, 150 and 200 uM. All concentrations are expressed per monomeric units. (B) Overlay of 3 1 P NMR spectra of pC measured at temperatures 10°C (blue), 20°C (magenta) and 37°C (red) with schematic images of respective species above the corresponding peaks. (C) Dependence of values (from the 3 1 P NMR data) on temperature with linear regression line. Adapted from Paper 1. 47 Results and Discussion 4.1 14-3-3 dimerization A B 0.01 0.1 1 10 100 1000 1 2 3 4 5 M (MM) Kd (|JM) Figure 4.5: Determination of pC/C heterodimerization dissociation constant. (A) Experimental data and fitted curve of FRET heterodimerization equilibrium assay. The total concentrations of donor- and acceptor-labelled £ is kept constant at 20 nM, while the concentration of p£ is varied. (B) Error analysis of the experiment. The solid line shows normalized sum of squares of errors, while the dashed line shows a threshold of 0.8. Adapted from Paper 1. state functions were derived from the temperature dependence using Eq. 1.12. The value of dissociation constant at 37°C was calculated as K<\ = (4.3 ± 0.3) mM. The Gibbs free energy at 37 °C is equal to (13.8 ± 0.7) kJ mol"1 . To determine the heterodimeric dissociation constant between phosphorylated and non-phosphorylated £, we slightly modified the FRET equilibrium assay. The equilibrated mixture of C _ D / C _ A in equimolar ratio was titrated with nonlabelled p(. The observed fluorescence intensity decreases with increasing concentration of p£ leading to formation of the pC/C heterodimer and dissociation of the FRET pair C _ D / C _ A , see Fig. 4.5. The heterodimerization value was determined to (2.5 ± 0.3) uM. The very high dissociation constant Kd = 4.3 mM of phosphorylated 14-3-3^ does not automatically imply that it would be purely in monomeric form in vivo. If there is non-phosphorylated £ present, the heterodimers pC/C will occur to some extent given by the heterodimerization Kd- The concentrations of individual forms are given below in Tab. 4.3. The results show that under estimated 48 Results and Discussion 4.1 14-3-3 dimerization physiological conditions (10% phosphorylated), about 15% of all pC occurs in dimeric form, bound into pC/C dimer. 4.1.6 14-3-3cr dimerization The kinetic FRET assay were employed also for another 14-3-3 isoform, namely 14-3-3<7, which plays an important role in some human cancers acting as a tumor suppressor [100]. The presented results were not published elsewhere yet. The homodimerization assay was measured at five temperatures ranging from 25 °C to 40 °C and fitted using COPASI with the same model as 14-3-3£, for illustration see Fig. 4.6A,B documenting the results at 37 °C. Both Kd and fc0ff values were successfully obtained from all the measured temperatures, allowing to fit the temperature dependencies and calculate the values of Gibbs free energy, entropy and enthalpy, presented in Tables 4.1, 4.2, and Fig. 4.7C,D. The CT/C heterodimerization was characterized using a modified FRET assay. Acceptor labelled f ,1> I 1 1 1 1 41.4 0.8 0.3 3.15 3.2 3.25 3.3 3.35 3.4 3.45 1 0 0 0 ^ 0 0 Figure 4.7: The effect of temperature on the £/C and a/a dimerization. (A) Rate constant fc0ff of C/C interaction was measured using both FRET (empty circles, dashed line) and SQ assay (filled circles, solid line). The data were fitted using Eyring equation Eq. 1.8 to provide transient thermodynamic parameters. (B) Dissociation constant of £C dimerization measured by FRET assay. The data were fitted using Eq. 1.12 to get the equilibrium thermodynamic parameters. (C,D) Temperature dependence of k0s and Kd of a/a homodimerization measured by FRET. Data were analysed analogically to the £/£ case. The errorbars were omitted in case they were smaller than the markers. 53 Table 4.1: Thermodynamic equilibrium parameters of 14-3-3 dimerization. Standard state (c° = 1M and T = 310 K) is assumed. The thermodynamic values for pC/C heterodimer were not determined. interacting Kd AG* AH* AS* isoforms [mol dm- 3 ] [kJ mol"1 ] [kJ mol"1 ] [J m o l ^ K - 1 ] c/c (4.6 ±0.7) x 10"9 48.8 ±0.1 40 ± 7 -29 ± 2 2 P C / P C (4.3 ±0.2) x 10"3 13.8 ±0.7 -38.3 ±0.7 -168 ± 3 P C / C (2.5 ±0.3) x 10"6 33.2 ±0.3 n.d. n.d. a/a (2.4 ± 1.3) x 10"9 51 ± 2 95 ± 3 1 150 ± 100 (9 ± 5 ) x 10"9 48 ± 2 n.d. n.d. Table 4.2: Thermodynamic kinetic parameters of 14-3-3 dimerization. Values are given for temperature T = 310 K. interacting k0f[ AGi AH* AS* isoforms [lO^s"1 ] [kJmol"1 ] [kJmol"1 ] [J m o l " 1 ^ 1 ] C/C 2.8 ±0.4 91.2 ±0.4 110 ± 6 60 ± 2 0 « 1.0 g OJ o sz OJ o o =3 ( » Lj~J ka = (3.3±0.5)x103 s'1 r0(l = 5.1+0.7 min ] I I I L 0 10 20 30 40 50 60 Time (min) 100 E- 10 = 0.01 n 1 0.1 s 2.0 2.5 3.0 3.5 4.0 4.5 5.0 3 „-1, 0.95 0.9 0.85 0.8 koff (10J s CO o S c 0) ocz CD o o 400 300 200 100 0 - 4 r ) c - _ » ,ioox L - ^ J — / l x L — J \ = 3.6+2.4 nM * A;(„ = (2.77±0.09)x10-3 S', i TM = 6.0±0.2 min I I I I 100 m 0 20 40 60 80 100 120 Time (min) 0.01 n 1 - 0.95 0.9 0.85 0.8 2.5 2.6 2.7 2.8 2.9 3.0 koff (10"3 s-1 ) Figure 4. Design of 14-3-3C dimerization kinetic assays. SQ (a) and FRET (c) kinetic profiles at 37 °C. The dissociation rate constant /coff, dimer mean lifetime T0n and Ka (in case of FRET assay) are given. At each stage, a schematic representation of the most abundant species is included. The fit (red) is overlaid with experimental data (black). Dependence of sum square error (SSE) on dissociation rate constant /coff and dissociation constant Kd from the SQ (b) and FRET (d) assay. Minimum SSE (indicated by the symbol'+') divided by the SSE obtained at each x,y coordinate (SSEm i n /SSEx y ) is presented. The rate constant corresponds to a mean dimer lifetime of (6.0 ± 0.4) min. External factors modulating 14-3-3; dimerization To reveal the influence of various biophysical determinants on 14-3-3C dissociation, we measured the dependence of koii and Kd on temperature, presence of a binding partner (regulatory domain of human tyrosine hydroxylase phosphorylated at S19 and S40, dpRD-hTH1), salt concentration, and pH. While kD» values were determined by both SQ and FRET kinetic assays except for the binding partner presence, the presented Kd values were obtained from the FRET kinetic assay for the reasons listed in previous section. Temperature. The temperature change from 313 K to 293 K led to a decrease of kon value from 4.3x10~3 s-1 to 0.25 x 10"3 S"1 . The obtained data (Figure 5(a)) were fitted using the Eyring equation42 in its logarithmic form (Eq. (2)). This fitting provided the transient thermodynamic parameters (enthalpy AH* and entropy AS*) of the dissociation free energy barrier (Table 1). The Gibbs free energy of the barrier at temperature T was calculated as AG*(T) = AH* - TASK AH* 1 AS* , fkBT\ ln/cofl = - —T + — + ln(^—j, (2) where R stands for universal gas constant, T for thermodynamic temperature, kB for Boltzmann constant and h for Planck constant. In addition to a slower dissociation rate, the decrease in temperature also led to a small Z Trosanova, P. Lousa, A. Kozelekova, et al. Journal of Molecular Biology 434 (2022) 167479 o FRET - - © - - 5s SQ i 1 0.0032 0.0033 0.0034 1/T(K-1 ) 0.0032 0.0033 0.0034 1/T(K"1 ) - • I I I S I I 0 2 4 6 8 10 [hTH] (ekv.) t 15 10 200 400 600 [NaCI] (mil) JLA 200 400 [NaCI] (mM) Figure 5. Homodimeric 14-3-3C dissociation rate constant /coff (a-d) and dissociation constant Kd (e-h) dependence on external factors: temperature (a,e), binding partner (dpRD-hTH1) concentration (b,f), NaCI concentration (c,g), and pH (d,h). The temperature dependence is plotted as the logarithm of the rate constant (dissociation constant, respectively) on the reciprocal temperature. The error bars were omitted in cases where they are smaller than symbol size. Presented data were obtained from SQ (closed circles, solid line) and FRET kinetic assays (empty circles, dashed line). reduction in K"d values (Figure 5(e)). The thermodynamic equilibrium parameters were fitted using Eq. (1). The results are shown in Table 1 and free energy profiles are presented in SI (Figure S3.10). Binding partner presence. To evaluate the influence of a binding partner on the dissociation of 14-3-3C protein, we performed the FRET kinetic assay in the presence of a known protein client, the regulatory domain of tyrosine hydroxylase 1 phosphorylated at S19 and S40 (dpRD-hTH1 ) 4 3 ~ 4 5 The FRET kinetic assay was performed as described above, with dpRD-hTH1 (0 to 10-fold excess, c = 0, 20, 40, 100 and 200 nM) present from the start of each experiment. The assay was performed three times at each listed partner concentration. Figure 5(b) shows that increasing concentration of the binding partner led to an exponential decrease of the dissociation rate constant at 37 °C from 2.8x10~3 s~1 to 0.13 x 10~3 s~1 , i.e. more than 20-fold. Also, a decrease in Kd values (down to subnanomolar values) was observed (Figure 5 (f))Buffer ionic strength and pH. The 14-3-3C dimer dissociation rate constants were also determined in varying buffer ionic strength and pH conditions. The increase of NaCI concentration from 0 mM to 600 mM lowered the dissociation rate constant (/(off) by a factor of 2 (Figure 5(c)), while having an insignificant effect on Kd values (Figure 5(g)). The change in pH in range of 5.7 to 9 did not alter /coff (Figure 5(d)) nor Kd (Figure 5(h)) in any significant and/or consistent way. Discussion Members of the 14-3-3 protein family are capable of forming a variety of homo and heterodimers that exist in equilibrium with the corresponding monomers.46 The populations of individual oligomeric forms depend on the expression levels of individual 14-3-3 isoforms (including post-translational modifications) and the dissociation constants (K"d) of all relevant homodimerizations and heterodimerizations. Surprisingly, the dimerization dissociation constants are unknown despite the high abundance of 14-3-3 isoforms in the human brain and their regulatory role with hundreds of client proteins.47,48 To our knowledge, this is the first study providing quantitative Kd values for any 14-3-3 isoform (£ in particular) and its phosphorylated form. The presented assays accompanied with theoretical models for their analysis have general applicability not only for other 14-3-3 isoforms but also for other protein families capable to form homo- and heterodimers. Impact of Ser58 phosphorylation on 14-3-3 monomer-dimer equilibria Non-phosphorylated 14-3-3 proteins form stable dimers, however the impact of phosphorylation on Z. Trosanova, P. Lousa, A. Kozelekova, et al. Journal of Molecular Biology 434 (2022) 167479 the oligomerization state is still being debated. There are a few key issues in previously published studies. First, phosphomimicking or monomeric mutants were used, which have revealed themselves as an inadequate approximation to phosphorylated protein in terms of oligomeric behavior. 4 , 3 1 '3 3 '4 & _ 5 1 Second, the phosphorylated 14-3-3 was often contaminated by nonphosphorylated protein, due to insufficient cleanup after phosphorylation or measurements being conducted in vivo. In those studies, the presence of non-phosphorylated protein could significantly influence the oligomeric state, which was not considered.24 '26 '28 We revealed that 14-3-3£ forms tight homodimers with the Kd of (5.5 ± 0.8) nM as determined by fluorescence and MST assays. On the other hand, 14-3-3£ phosphorylated at Ser58 (without the presence of the non-phosphorylated form) has the homodimerization Kd equal to (4.3 ± 0.3) mM as determined by solution 3 1 P NMR spectroscopy. Given these dissociation constants, we could expect that at micromolar concentrations, £ exists dominantly in the dimeric, while p£ in the monomeric state. However, assuming an estimated upper limit of 14-3-3£ concentration in the human brain3 7 '3 8 , 5 2 , 5 3 being 100 uM, and a certain level of phosphorylation, one cannot neglect the formation of p£/£ heterodimers. Here, we identified the heterodimerization Kd (p£/£) being (2.5 ± 0.3) uM which is highly relevant at physiological 14-3-3 concentrations. This value also suggests that previous studies of p£ in the presence of £ observed p£/£ heterodimers that could be mistaken for p£/p£ homodimers. To illustrate the importance of heterodimerization, we apply the determined Kd values and assume a 10% level of £ phosphorylation (i.e. [£] = 90 LiM and Jp£] = 10LiM), being in physiological range. • Under these conditions, we expect the presence of the following entities at the following concentrations: £/£ (44 uM), £ (0.45 uM), p£ (8.5 uM), p£/£ (1.5 uM), and p£/p£ (0.02 uM). It follows that a significant portion (approx. 15%) of p£ is in the heterodimeric state (p£/£) in physiological concentration range. The high degree of homology between 14-3-3 isoforms and the fact that most of them have a phosphorylatable serine at the position corresponding to Ser58 of £, leads to another important point. We expect a similar impact of phosphorylation on a wide range of 14-3-3 isoforms as we described here for £, i.e. strong monomerization corresponding to the Kd in millimolar range, followed by heterodimerization with the non-phosphorylated isoforms, with Kd of the heterodimers in the low micromolar range. The obtained Kd,k0n and kon values allow the characterization of 14-3-3£ dimer and monomer populations as well as their mean life-times (^off, Ton) for any given total concentration. In Section 2.6 in SI, we provide a visualization of the oligomeric populations for any percentage of the phosphorylated protein (range 0-100%), based on the Kd values presented here. The equilibria are dynamic and, on average, it takes about 6.5 min (T0ff = 1 /kon) for the £/£ dimer to dissociate and about 4 s (T0n = [M]//c0n) for a £ monomer to associate with another 14-3-3£ monomer into a dimer. Note that in a complex environment (e.g. intracellularly) the £ and p£ will also exist as heterodimers with other isoforms. Effect of external conditions Figure 5 shows that the only considerable change of 14-3-3£ homodimerization Kd (by a factor of ca. 20) was observed in the presence of its binding partner (dpRD-hTH143 ) at ratios above 5:1 (corresponding to 30% of bound 14-3-3£ based on Kd = 400 nM54 ). It is well established that 14-3-3 proteins change the conformation of their binding partners.55-57 Here we show (to our knowledge for the first time) that such interaction has a significant impact on 14-3-3 itself by stabilizing it in the dimeric state. The novel observation that a binding partner has also a marked physical effect on the 14-3-3 protein might have additional downstream consequences considering 14-3-3 as a hub in a large signaling network, binding hundreds of phosphorylated proteins5 8 The effect of temperature on Kd and koii allows us to determine the values of both equilibrium and transient-state enthalpy, entropy and Gibbs free energy. In case of non-phosphorylated £ dimerization, the main contribution to free energy comes from enthalpy, which being negative leads to increase of Kd with higher temperature. Surprisingly, the dimerization process for phosphorylated p£ is different. In this case, the Kd decreases with increasing temperature, therefore the enthalpy is positive and the dimerization is driven by entropy. The sign change in enthalpy might be explained by the strong repulsion between negatively charged phosphoserines (at position 58) at both monomers (Figure 1(a)). The change in entropic part could be explained by destabilization of the dimeric interface due to the presence of pSer58. The phosphorylation probably also impacts the properties of the water shell around the interface. The minor influence of salt on the dissociation rate constants (both Kd and kDn) indicates that the overall dimer stability originates predominantly from hydrophobic interactions, being to some extent compensated by polar contacts. There are two complementary hydrophobic patches at the 14-3-3£ dimer interface, around residues L12 and M78. To support our hypothesis regarding their major role in dimer stabilization, we constructed a double-point mutant L12K/M78E.39 Such a construct may in principle form one additional intersubunit salt bridge between K12 and E78, while Z Trosanova, P. Lousa, A. Kozelekova, et al. Journal of Molecular Biology 434 (2022) 167479 retaining the overall 14-3-3C charge. However, the oligomerization analysis of this construct clearly displays the monomeric state at micromolar concentrations (Figure 1). Design of 14-3-3 monomer-dimer assays We present two types of assays addressing the dimerization of tight dimers: the equilibrium assays for steady-state experiments; and kinetic assays assessing dynamic properties. Both types utilized labeling by fluorescent dyes attached through a flexible linker. All three equilibrium assays (SQ, FRET, MST) provided Kd values in the low nanomolar range averaged to Kd = (5.5 ± 0.8) nM. The results are very consistent, even though the individual assays differ in their physical principles. In SQ and FRET assays, the oligomeric state of 14-3-3C influences the intensity of fluorescence. In the MST assay, the oligomeric state affects the thermophoretic mobility. The MST was evaluated using a general monomer-dimer model (as described in SI, Section 2.4), which differs significantly from the receptor-ligand model provided by the manufacturer (NanoTemper). The model applied to fitting the FRET kinetic assay data provided a Kd value of (3.6 ± 2.4) nM, as well as the kinetic parameters. While the Kd value is consistent with the values obtained from equilibrium assays, the uncertainty from the FRET kinetic assay (Figure 4(c)) is larger. On the other hand, this assay requires lower sample consumption (factor of 6) and can be performed much faster (ca. 2 vs. 10 h). For these reasons, the FRET kinetic assay was chosen to measure the dependence of Kd on temperature, binding partner presence, ionic strength, and pH. The same fluorescence assay design as presented for tight 14-3-3C homodimers, can be applied to any other 14-3-3 isoform. For heterodimerization purposes, the FRET assay is particularly useful as it allows the labeling of one 14-3-3 isoform with an acceptor, while the other with a donor molecule after which the two samples are mixed. We proved the applicability for the pC7( heterodimerization. Notably, the presented approaches do not require any special instrumentation - the fluorescence assays can be performed on regular cuvette fluorometers and 1D 3 1 P NMR spectra can be acquired using NMR spectrometers available at most chemistry departments. The determination of Kd for all 14-3-3 dimers, coupled with an increasing number of proteomic studies of 14-3-3 expression profiles in various tissues or diseases, will enable to characterise the absolute populations of all present 14-3-3 homoand heterodimers as well as monomers. For this purpose, the corresponding homo- and heterodimeric Kd values are desired, warranting further research. Quantitative description of this intricate network could be of particular interest to computational biologists as a basis for predictive analyses or advanced experimental design. Conclusions The oligomeric state of 14-3-3C, a major isoform in the human brain, was quantitatively examined by a variety of methods. For quantitative analysis of the equilibrium and kinetics of tight 14-3-3C dimerization, we designed fluorescence assays based on the FRET and self-quenching phenomena. Phosphorylation at Ser58 impacts the 14-3-3C dimerization most dramatically, shifting the homodimerization Kd by 6 orders of magnitude from the value of (5.5 ± 0.8) nM (£/£) to (4.3 ± 0.3) imM (pC/pQ. The heterodimerization Kd of pC/C was determined as (2.5 ± 0.3) uM implying that at typical expression levels of pC, about 85% exist as monomers and approximately 15% exist in the heterodimeric state (pC/Q. The kinetic studies provided the 14-3-3C homodimer dissociation rate constant of 92.8 ±0.4) x 10~3 s~\ corresponding to a mean dimer lifetime of (6.0 ± 0.4) min at 37 °C. The £/£ dissociation free-energy barrier of 90 kJ mol- 1 is dominated by the enthalpic contribution of 110 kJ mol- 1 . Assuming the 14-3-3C concentration in human brain being 100 u.M, the concentration of its monomeric form is expected to be 0.5 uM with a mean life-time of 4 s. The 14-3-3C dimer dissociation constant Kd is almost insensitive to temperature, ionic strength and pH. On the other hand, the dimerization affinity increases significantly after addition of a binding partner (dpRD-hTH1), which suggests possible regulation of the 14-3-3 protein by its client proteins. The presented self-quenching and FRET based assays can be employed to determine homo- and heterodimerization dissociation constants of other proteins, for instance of other 14-3-3 isoforms. Materials and Methods Design of constructs A codon-optimized cDNA for 14-3-3C (GenScript) with two mutations C25A and C189A (solvent accessible cysteines) was inserted into a pET15b plasmid including a TEV cleavable His-tag.59 This construct (called 14-3-3C in the text) was used as the parental construct for the monomeric mutant and Ntail construct. In the monomeric construct, mutations L12K and M78E were introduced using site-directed mutagenesis (QuikChange, Agilent). For specific fluorescent labeling, we designed a 14-3-3Ntail construct, where we inserted an SVDACKGSSGG sequence preceding the Nterminus. Sequences of the constructs were verified by Sanger sequencing. Proteins were expressed Z. Trosanova, P. Lousa, A. Kozelekova, et al. Journal of Molecular Biology 434 (2022) 167479 and purified as described previously, the most important steps are outlined in SI (Section 1). Preparation of phosphorylated 14-3-3C 14-3-3C protein (WT or Ntail) was phosphorylated using recombinant PKA (for PKA preparation protocol, see SI Section 1.2) in PKA buffer (20 mM Tris, 15 mM MgCI2, 3 mM NaN3; pH = 7.4) with 1 mM ATP. The phosphorylation mixture was incubated at 37 °C for 4 h. Afterwards, phosphorylated protein was isolated using anion exchange chromatography (HiTrap Q HP, GE Healthcare). Protein identity, purity and oligomeric state were verified by MALDI-TOF MS, SDS-PAGE and Native PAGE, respectively. According to LC-MS/MS analysis, the final sample was quantitatively phosphorylated at Ser58 (99.8%) and additional minor phosphorylation was detected at Ser28 (3-4%). For detailed protocol of pS58_14-3-3C preparation, see SI, Section 1.3. Maleimide-dye labeling procedure The Ntail constructs were specifically labeled using a 30-fold excess of the fluorescent dyes, AlexaFluor647-C5-maleimide (AF647), AlexaFluor488-C2-maleimide (AF488, both Thermo Fisher Scientific) and 6-tetramethylrhoda mine-C6-maleimide (TMR, AAT Bioquest) for 30 min at room temperature in the dark. Before labeling, dialyzed protein was treated with 10 mM TCEP for 30 min at room temperature to remove potential disulfide bonds. After labeling, unreacted free dye was removed using a molecular sieve Vivaspin6 MWCO 10 or 30 kDa (GE Healthcare), depending on protein oligomeric state. The completion of labeling was confirmed by MALDITOF analysis (SI, Section 1.4). Native polyacrylamide gel electrophoresis (Native PAGE) Protein samples were mixed with loading buffer (final composition: 1% glycerol, 1% Bromophenol Blue, 80 mM TRIS-HCI; pH = 6.8) and were loaded onto a freshly prepared 12.5% native minigel. The gel was run for 200 min at a constant voltage of 95 V in a native electrophoretic buffer (2.5 mM TRIS-HCI, 19.2 mM glycine; pH = 8.3). In order to prevent the thermal denaturation of studied proteins, the apparatus was cooled down on ice during the whole procedure. Subsequently, the gel was stained by Coomassie Brilliant Blue R- 250 (AppliChem, Darmstadt, Germany). Fluorescence measurements setup Fluorescence measurements were conducted employing a FluoroMax-4 Spectrofluorometer (HORIBA Jobin Yvon). Before each measurement, the quartz cuvette was treated with 10 mg/mL BSA for 30 min to avoid adhesion of the fluorescent dyes to the walls. For the SQ assay, lex = 553 nm (bandwidth of 0.8 nm) and l e m = 575 nm (bandwidth 2.5 nm) were used. For the FRET assay, l e x = 470 nm (bandwidth 4.2 nm) and lem = 666 nm (bandwidth 10.5 nm) were used. Unless specified otherwise, the measurement was performed at 37 °C in 20 mM NaPi buffer, pH = 6.8. Microscale Thermophoresis Microscale Thermophoresis (MST) assay was performed at constant concentration of labeled 14- 3-3C_AF647 (c = 0.5 nM) and varying concentration of unlabeled 14-3-3C (c = 6 pM- 100 nM) in a buffer containing 20 mM phosphate, 0.5 mg/mL BSA and 0.05% TWEEN-20 at 30 °C. Binding studies were performed in triplicates using 50% LED and 80% MST power with a Monolith NT.115Pico device (NanoTemper Technologies) using standard capillaries. Fluorescence and MST data analysis All data gathered from fluorescence assays and MST were analysed within COPASI software v4.2461 using a set of chemical equations describing the equilibria among several labeled and unlabeled variants of £ and p£. All used models are described in detail in SI, Sections 2 and 3. The models consist of several phases described by kinetic differential equation, separated by additions of labeled or unlabeled protein samples. By numerical integration, time dependencies of concentrations of all chemical species are obtained. The signal is then calculated (using Eqs. S3.6 for SQ and S3.16 for FRET) and fitted to the experimental data. The confidence intervals were estimated using a method described by Johnson etal.,62 by evaluation of sum of squared errors (SSE) in dependency on Kd and/or kon. To obtain the SSE contour plot, the parameter(s) were systematically varied over a range close to the best fit, while adjusting the other parameters (i.e. fluorescence scaling factors). For each combination of parameters x, y, after the best possible adjustment was achieved, the S S E x y was calculated. The upper and lower limits of parameters were obtained at a threshold of 0.8 at the reciprocal normalized SSE plot (calculated as SSEmin/SSEXj,). The resulting plots and contour plots are shown in SI (Sections 2 and 3). NMR experiments and analysis The 3 1 P spectra were measured over the course of 18 h at 202.49 MHz using Bruker Avance NEO 500 MHz spectrometer equipped with a dual-band Prodigy cryoprobe. The pC sample concentration was 1.1 mM. The peaks were fitted by Lorentz curves using a built-in deconvolution function in TopSpin software v4.0.8 (Bruker). The obtained peak areas (including confidence intervals) were Z Trošanová, P. Louša, A. Kozeleková, et al. Journal of Molecular Biology 434 (2022) 167479 used directly to calculate the resulting Kd together with its confidence intervals. DECLARATION OF COMPETING INTEREST The authors declare that they have no k nown competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement This research was financed by MEYS-CR InterExcellence Inter-Action grant scheme (No. LTAUSA18168) and by the Czech Science Foundation (No. GF20-05789L). CIISB research infrastructure project LM2015043, funded by Ministry of Education, Youth and Sports of the Czech Republic (MEYS CR), is gratefully ack nowledged for partial financial support of the measurements at Biomolecular Interactions and Crystallization and Proteomics Core Facilities, at the Central European Institute of Technology (CEITEC), Masaryk University. We also ack nowledge the MSMTDAAD mobility grant (project No. 7AMB16DE005). The cDNA of the PKA catalytic domain was k indly provided by assoc. prof. Lumi'r Krejf (Faculty of Medicine, Masaryk University). We thank dr. Karel Drbal, Department of Cell Biology, Faculty of Science, Charles University in Prague, for access to Monolith NT.115 Pico device purchased from Grant No. CZ. 1.05/4.1.00/16.0347. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at https://doi. org/10.1016/j .jmb.2022.167479. Received 22 October 2021; Accepted 31 January 2022; Available online 5 February 2022 Keywords: 14-3-3; phosphorylation; dimerization; dissociation constant; FRET; NMR f Equal contribution. References 1. Steinacker, P., Aitken, A., Otto, M., (2011). 14­3­3 proteins in neurodegeneration. S eminars Cell Develop. Biol. 22, 696­704. 2. Mackintosh, C , (2004). Dynamic interactions between 14­ 3­3 proteins and phosphoproteins regulate diverse cellular processes. Biochem. J. 381, 329­342. 3. Berg, D., Holzmann, O , Riess, O., (2003). 14­3­3 proteins in the nervous system. Nature Rev. Neurosci. 4, 752­762. 4. Sluchanko, N.N., Gusev, N.B., (2011). Probable participation of 14­3­3 in tau protein oligomerization and aggregation. J. Alzheimers Dis. 27, 467­476. 5. Sluchanko, N.N., Gusev, N.B., (2010). 14­3­3 proteins and regulation of cytoskeleton. Biochemistry (Moscow) 75, 1528­1546. 6. Zhao, J., Meyerkord, C L , Du, Y., Khuri, F.R., Fu, H., (2011). 14­3­3 proteins as potential therapeutic targets. Seminars Cell Develop. Biol. 22, 705­712. 7. Freeman, A.K., Morrison, D.K., (2011). 14­3­3 Proteins: Diverse functions in cell proliferation and cancer progression. S eminars Cell Develop. Biol. 22, 681­687. 8. Gan, Y., Ye, F., He, X.X., (2020). The role of YWHAZ in cancer: A maze of opportunities and challenges. J. Cancer 11, 2252­2264. 9. Uhart, M., Bustos, D.M., (2013). Human 14­3­3 paralogs differences uncovered by cross­talk of phosphorylation and lysine acetylation. PLoS ONE 8, e55703. 10. Yaffe, M.B., Rittinger, K., Volinia, S., Caron, P.R., Aitken, A., L effers, H., Gamblin, S.J., Smerdon, S.J., etal., (1997). The structural basis for 14­3­3:phosphopeptide binding specificity. Ce//91, 961­971. 11. Xiao, B., Smerdon, S.J., Jones, D.H., Dodson, G.G., Soneji, Y., Aitken, A., Gamblin, S.J., (1995). Structure of a 14­3­3 protein and implications for coordination of multiple signalling pathways. Nature376, 188­191. 12. L iu, D., Bienkowska, J., Petosa, O , Collier, R.J., Fu, H., Liddington, R., (1995). Crystal structure of the zeta isoform of the 14­3­3 protein. Nature 376, 191­194. 13. Yang, X., L ee, W.H., Sobott, F., Papagrigoriou, E., Robinson, C.V., Grossmann, J.G., Sundstróm, M., Doyle, D.A., et al., (2006). Structural basis for protein­protein interactions in the 14­3­3 protein family. Proc. Natl. Acad. Sci. US A 103, 17237­17242. 14. L iu, J.Y., Li, Z., Li, H., Zhang, J.T., (2011). Critical residue that promotes protein dimerization: A story of partially exposed Phe25in 14­3­3CT. J. Chem. Inf. Model. 51, 2612­ 2625. 15. Yaffe, M.B., (2002). How do 14­3­3 proteins work? Gatekeeper phosphorylation and the molecular anvil hypothesis. FEBS Lett. 513, 53­57. 16. Tzivion, G., L uo, Z., Avruch, J., (1998). A dimeric 14­3­3 protein is an essential cofactor for Raf kinase activity. Nature 394, 88­92. 17. L andrieu, I., L acosse, L ., L eroy, A., Wieruszeski, J.­M., Trivelli, X., Sibille, N., Schwalbe, H., Saxena, K., et al., (2006). NMR analysis of a tau phosphorylation pattern. J. Am. Chem. S oc. 128, 3575­3583. 18. Jansen, S., Melková, K., Trošanová, Z., Hanáková, K., Zachrdla, M., Nováček, J., Župa, E., Zdráhal, Z., et al., (2017). Quantitative mapping of microtubule­associated protein 2c (MAP2c) phosphorylation and regulatory protein 14­3­3£­binding sites reveals key differences between MAP2c and its homolog Tau. J. Biol. Chem. 292, 6715­ 6727. 19. Sluchanko, N.N., Seit­Nebi, A.S., Gusev, N.B., (2009). Effect of phosphorylation on interaction of human tau protein with 14­3­3J. Biochem. Biophys. Res. Commun. 379, 990­994. Z. Trosanovä, P. Lousa, A. Kozelekova', et al. 20. Obsil, T., Ghirlando, R., Klein, D.C., Ganguly, S., Dyda, F., (2001). Crystal structure of the 14-3-3{:serotonin Nacetyltransferase complex: a role for scaffolding in enzyme regulation. Ce//105, 257-267. 21. Aitken, A., (2011). Post-translational modification of 14-3-3 isoforms and regulation of cellular function. Seminars Cell Develop. Biol. 22, 673-680. 22. Megidish, T., Cooper, J., Zhang, L , Fu, H., Hakomori, S.I., (1998). A novel sphingosine-dependent protein kinase (SDK1) specifically phosphorylates certain isoforms of 14- 3-3 protein. J. Biol. Chem. 273, 21834-21845. 23. Woodcock, J.M., Murphy, J., Stomski, F.C., Berndt, M.C., Lopez, A.F., (2003). The dimeric versus monomeric status of 14-3-3fis controlled by phosphorylation of Ser58 at the dimer interface. J. Biol. Chem. 278, 36323-36327. 24. Gu, Y.-M., Jin, Y.H., Choi, J.K., Baek, K.H., Yeo, C.Y., Lee, K.Y., (2006). Protein kinase A phosphorylates and regulates dimerization of 14-3-3J. FEBS Lett. 580, 305- 310. 25. Powell, D.W., Rane, M.J., Chen, Q., Singh, S., McLeish, K. R., (2002). Identification of 14-3-3{as a protein kinase B/ Akt substrate. J. Biol. Chem. 277, 21639-21642. 26. Gerst, F., Kaiser, G., Panse, M., Sartorius, T., Pujol, A., Hennige, A.M., Machicao, F., Lammers, R., et al., (2015). Protein kinase C<5regulates nuclear export of FOX01 through phosphorylation of the chaperone 14-3-3J. Diabetologia 58, 2819-2831. 27. Kim, Y.S., Choi, M.Y., Kim, Y.H., Jeon, B.T., Lee, D.H., Roh, G.S., Kang, S.S., Kim, H J . , et al., (2010). Protein kinase Cdelta is associated with 14-3-3 phosphorylation in seizure-induced neuronal death. Epilepsy Res. 92, 30-40. 28. Powell, D.W., Rane, M.J., Joughin, B.A., Kalmukova, R., Hong, J.-H., Tidor, B., Dean, W.L., Pierce, W.M., et al., (2003). Proteomic identification of 14-3-3fas a mitogenactivated protein kinase-activated protein kinase 2 substrate: Role in dimer formation and ligand binding. Mol. Cell. Biol. 23, 5376-5387. 29. Zhou, J., Shao, Z., Kerkela, R., Ichijo, H., Muslin, A.J., Pombo, C , Force, T., (2009). Serine 58 of 14-3-3{is a molecular switch regulating ASK1 and oxidant stressinduced cell death. Mol. Cell. Biol. 29, 4167-4176. 30. Kanno, T., Nishizaki, T., (2011). Sphingosine induces apoptosis in hippocampal neurons and astrocytes by activating caspase-3/-9 via a mitochondrial pathway linked to SDK/14-3-3 protein/Bax/cytochrome c. J. Cell. Physiol. 226, 2329-2337. 31. Civiero, L , Cogo, S., Kiekens, A., Morganti, C , Tessari, I., Lobbestael, E., Baekelandt, V., Taymans, J.M., et al., (2017). PAK6 phosphorylates 14-3-3yto regulate steady state phosphorylation of LRRK2. Front. Mol. Neurosci. 10, 417. 32. Gökirmak, T., Denison, F.C., Laughner, B.J., Paul, A.-L.L., Ferl, R.J., (2015). Phosphomimetic mutation of a conserved serine residue in Arabidopsis thaliana 14-3- 3cosuggests a regulatory role of phosphorylation in dimerization and target interactions. Plant Physiol. Biochem. 97, 296-303. 33. Sluchanko, N.N., Chernik, I.S., Seit-Nebi, A.S., Pivovarova, A.V., Levitsky, D.I., Gusev, N.B., (2008). Effect of mutations mimicking phosphorylation on the structure and properties of human 14-3-3f. Arch. Biochem. Biophys. 477, 305-312. 34. Denison, F.C., Gökirmak, T., Ferl, R.J., (2014). Phosphorylation-related modification at the dimer Journal of Molecular Biology 434 (2022) 167479 interface of 14-3-3codramatically alters monomer interaction dynamics. Arch. Biochem. Biophys. 541, 1-12. 35. Woodcock, J.M., Goodwin, KL., Sandow, J.J., Coolen, C , Perugini, M.A., Webb, A.I., Pitson, S.M., Lopez, A.F., et al., (2018). Role of salt bridges in the dimer interface of 14-3- 3fin dimer dynamics, N-terminal a-helical order, and molecular chaperone activity. J. Biol. Chem. 293, 89-99. 36. Jones, D.H., Ley, S., Aitken, A., (1995). Isoforms of 14-3-3 protein can form homo- and heterodimers in vivo and in vitro: implications for function as adapter proteins. FEBS Lett. 368, 55-58. 37. Gogl, G., Tugaeva, K.V., Eberling, P., Kostmann, C , Trave, G., Sluchanko, N.N., (2021). Hierarchized phosphotarget binding by the seven human 14-3-3 isoforms. Nature Commun. 12, 1677. 38. Wang, M., Herrmann, C.J., Simonovic, M., Szklarczyk, D., Mering, C , (2015). Version 4.0 of PaxDb: Protein abundance data, integrated across model organisms, tissues, and cell-lines. Proteomics 15, 3163-3168. 39. Jandovä, Z., Trosanovä, Z., Weisovä, V., Oostenbrink, C , Hritz, J., (1866, 2018,). Free energy calculations on the stability of the 14-3-3fprotein. Biochim. Biophys. Acta (BBA) Proteins Proteomics, 442-450. 40. Duhr, S., Braun, D., (2006). Why molecules move along a temperature gradient. Proc. Natl. Acad. Sei. USA 103, 1968-1972. 41. Duhr, S., Braun, D., (2006). Thermophoretic depletion follows Boltzmann distribution. Phys. Rev. Lett. 96, 168301-168304. 42. Eyring, H., (1935). The activated complex in chemical reactions. J. Chem. Phys. 3, 107-115. 43. Ghorbani, S., Fossbakk, A., Jorge-Finnigan, A., Flydal, M. I., Haavik, J., Kleppe, R., (2016). Regulation of tyrosine hydroxylase is preserved across different homo- and heterodimeric 14-3-3 proteins. Amino Acids 48, 1221- 1229. 44. Kleppe, R., Toska, K., Haavik, J., (2001). Interaction of phosphorylated tyrosine hydroxylase with 14-3-3 proteins: evidence for a phosphoserine 40-dependent association. J. Neurochem. 77, 1097-1107. 45. Skjevik, A.A., Mileni, M., Baumann, A., Halskau, O., Teigen, K., Stevens, R.C., Martinez, A., (2014). The Nterminal sequence of tyrosine hydroxylase is a conformationally versatile motif that binds 14-3-3 proteins and membranes. J. Mol. Biol. 426, 150-168. 46. Sluchanko, N.N., Gusev, N.B., (2012). Oligomeric structure of 14-3-3 protein: What do we know about monomers? FEBS Lett. 586, 4249-4256. 47. Bustos, D.M., (2012). The role of protein disorder in the 14- 3-3 interaction network. Mol. BioSyst. 8, 178-184. 48. Tinti, M., Johnson, C , Toth, R., Ferrier, D.E., Mackintosh, C , (2012). Evolution of signal multiplexing by 14-3-3binding 2R-ohnologue protein families in the vertebrates. Open Biol. 2, 120103. 49. Sluchanko, N.N., Sudnitsyna, M.V., Chernik, I.S., SeitNebi, A.S., Gusev, N.B., (2011). Phosphomimicking mutations of human 14-3-3faffect its interaction with tau protein and small heat shock protein HspB6. Arch. Biochem. Biophys. 506, 24-34. 50. Sluchanko, N.N., Sudnitsyna, M.V., Seit-Nebi, A.S., Antson, A.A., Gusev, N.B., (2011). Properties of the monomeric form of human 14-3-3fprotein and its interaction with tau and HspB6. Biochemistry 50, 9797- 9808. 12 Z Trošanová, P. Louša, A. Kozeleková, et al. 51. Sluchanko, N.N., Uversky, V.N., (1854). Hidden disorder propensity of the N­terminal segment of universal adapter protein 14­3­3 is manifested in its monomeric form: novel insights into protein dimerization and multifunctionality. Biochim. Biophys. Acta (BBA) Proteins Proteomics 2015, 492­504. 52. Boston, P.F., Jackson, P., Thompson, R.J., (1982). Human 14­3­3 protein: radioimmunoassay, tissue distribution, and cerebrospinal fluid levels in patients with neurological disorders. J. Neurochem. 38, 1475­1482. 53. Ellis, R., (2001). Macromolecular crowding: an important but neglected aspect of the intracellular environment. Curr. Opin. S truct. Biol. 11, 114­119. 54. Obsilova, V., Nedbalkova, E., Silhan, J., Boura, E., Herman, P., Vecer, J., Sulc, M., Teisinger, J., et al., (2008). The 14­3­3 protein affects the conformation of the regulatory domain of human tyrosine hydroxylase. Biochemistry 47, 1768­1777. 55. Alblová, M., Šmídové, A., Dočekal, V., Veselý, J., Herman, P., Obšilová, V., Obšil, T., (2017). Molecular basis of the 14­3­3 protein­dependent activation of yeast neutral trehalase Nth1. Proc. Natl. Acad. S ci. US A 114, E9811­ E9820. 56. Sluchanko, N.N., Beelen, S., Kulikova, A.A., Weeks, S.D., Antson, A.A., Gusev, N.B., Strelkov, S.V., (2017). Structural basis for the interaction of a human small heat shock protein with the 14­3­3 universal signaling regulator. Structure 25, 305­316. Journal of Molecular Biology 434 (2022) 167479 57. Tugaeva, K.V., Titterington, J., Sotnikov, D.V., Maksimov, E.G., Antson, A.A., Sluchanko, N.N., (2020). Molecular basis for the recognition of steroidogenic acute regulatory protein by the 14­3­3 protein family. FEBS J. 287, 3944­ 3966. 58. Sluchanko, N.N., (2018). Association of multiple phosphorylated proteins with the 14­3­3 regulatory hubs: problems and perspectives. J. Mol. Biol. 430, 20­26. 59. Hritz, J., Byeon, l.­J.L , Krzysiak, T., Martinez, A., Sklenář, V., Gronenborn, A., (2014). Dissection of binding between a phosphorylated tyrosine hydroxylase peptide and 14­3­ 3£: a complex story elucidated by NMR. Biophys. J. 107, 2185­2194. 60. L ouša, P., Nedozrálová, H., Župa, E., Nováček, J., Hritz, J., (2017). Phosphorylation of the regulatory domain of human tyrosine hydroxylase 1 monitored using non­uniformly sampled NMR. Biophys. Chem. 223, 25­29. 61. Hoops, S., Sahle, S., Gauges, R., L ee, C , Pahle, J., Simus, N., Singhai, M., Xu, L, et al., (2006). COPASI­a COmplex PAthway Simulator. Bioinformatics 22, 3067­ 3074. 62. Johnson, K.A., Simpson, Z.B., Blom, T., (2009). FitSpace Explorer: an algorithm to evaluate multidimensional parameter space in fitting kinetic data. Anal. Biochem. 387, 30^1. 63. Sali, A., Blundell, T.L ., (1993). Comparative protein modelling by satisfaction of spatial restraints. J. Mol. Biol. 234, 779­815. 13 Paper 2 Louša P., Nedozrálová H., Župa E., Nováček J., Hritz J.* * Corresponding author Phosphorylation of the regulatory domain of human tyrosine hydroxylase 1 monitored using non­uniformly sampled N M R Biophysical Chemistry 2017, 223, 25­29 doi: 10.1016/j.bpc.2017.01.003 PL designed the experiments, performed the NMR experiments, processed and assigned the spectra, analysed the kinetic data, wrote and reviewed the manuscript. 97 Biophysical Chemistry 223 (2017) 25­29 E L S E V I E F Contents lists available at ScienceDirect Biophysical Chemistry journal homepage: http://www.elsevier.com/locate/biophyschem Phosphorylation of the regulatory domain of human tyrosine hydroxylase 1 monitored using non­uniformly sampled NMR Petr Louša, Hana Nedozrálová, Erik Župa, Jiří Nováček, Jozef Hritz * CEIWC MU, Masaryk University, Kamenice 753/5,625 00 Brno, Czech Republic I CrossMark H I G H L I G H T S G R A P H I C A L A BS T R A C T • Disordered part of regulatory domain of human tyrosine hydroxylase 1 was assigned. • Transient alpha-helices are present next to phosphorylation sites S40 and S19. • The secondary structure does not change after phosphorylation. • The phosphorylation kinetic rates were measured efficiently using time resolved NMR. pS40 RD­hTHl A R T I C L E I N F O Article history: Received 30 December 2016 Accepted 24 January 2017 Available online 27 January 2017 Keywords: Human tyrosine hydroxylase IDP NMR Phosphorylation Kinetics Time­resolved NMR Non­uniform sampling SSP A B S T R A C T Human tyrosine hydroxylase 1 (hTHl) activity is regulated by phosphorylation of its regulatory domain (RDhTHl ) and by an interaction with the 14-3-3 protein. The RD-hTHl is composed of a structured region (66- 169) preceded by an intrinsically disordered protein region (IDP, hTHl_65) containing two phosphorylation sites (S19 and S40) which are highly relevant for its increase in activity. The NMR signals of the IDP region in the non-phosphorylated, singly phosphorylated (pS 40) and doubly phosphorylated states (pS19_pS40) were assigned by non-uniformly sampled spectra with increased dimensionality (5D). The structural changes induced by phosphorylation were analyzed by means of secondary structure propensities. The phosphorylation kinetics of the S40 and S19 by kinases PKA and PRAK respectively were monitored by non-uniformly sampled time-resolved NMR spectroscopy followed by their quantitative analysis. © 2017 Elsevier B.V. All rights reserved. 1. Introduction Tyrosine hydroxylase (TH) is an enzyme which converts L­tyrosine to L-DOPA. This reaction is the rate limiting step in the biosynthetic pathway producing important catecholamine neurotransmitters: dopamine, noradrenaline and adrenaline [1,2]. The human enzyme isoform 1 is a tetramer each consisting of three domains: an N­terminal regulatory Corresponding author. E­mail address: jozef.hritz@ceitec.muni.cz (J. Hritz). domain (RD­hTHl, 1­169 aa), a catalytic domain (170­450 aa), and a short C­terminal tetramerization domain (451­497 aa) [3]. The first 65 residues of the regulatory domain form an intrinsically disordered protein region (IDP, hTHl_65) important for regulation. The activity of hTHl is controlled by the phosphorylation of its IDP region (S19, S31, S40) and by the interaction with 14­3­3 protein [4]. Phosphorylation sites SI 9 and S40 are the most relevant phosphorylation sites regarding 14­3­3 binding [5,6]. Recently, the NMR structure of the ordered region (65­159) of the dimeric regulatory domain of rat tyrosine hydroxylase (the sequence http://dx.doi.org/l 0.1016/j.bpc.2017.01.003 0301­4622/© 2017 Elsevier B.V. All rights reserved. 26 P. Lousa et at / Biophysical Chemistry 223 (2017) 25-29 identity between rat and human RD-hTHl is 81.8%) was determined by Zhang et al. [7]. Authors claimed that the problems with unstable sample and low signal dispersion of the 1DP region prevented structural characterization of the full RD. We overcame these problems by modifying the preparation of full length RD-hTHl (1-169) sample in non- and phosphorylated states and by applying non-uniform sampling (NUS) NMR approaches allowing much faster data collection in comparison with the uniformly sampled experiments. In the past, the rate of the phosphorylation of hTHl was monitored semi-quantitatively by Toska et al. [8] using radioactively labeled ATP. Such methodology has several drawbacks, especially the necessity of working with radioactive material, laborious sample preparation and low temporal resolution. The alternative methodology for a monitoring of protein phosphorylation is NMR spectroscopy where individual NMR spectra are collected sequentially during the course of a phosphorylation reaction [9]. The time resolution of such an approach is on the order of the measurement time of one particular NMR spectra. Another possibility for monitoring is to measure one long 2D spectrum and afterwards analyze lineshapes of signals modulated by the reaction kinetics [10]. In this way, the time resolution can be reduced at cost of significant increase in difficulty of analysis. The recent advances in NMR allow us to overcome these limitations by applying non-uniform sampled time-resolved NMR spectroscopy to monitor the reaction course with time resolution down to seconds [11]. This approach was already successfully utilized for the monitoring of the phosphorylation of cytoplasmic domain of human B cell receptor protein CD79b [11]. In this study, we present the resonance assignment of the 1DP region (hTHl_65) of RD-hTHl in the non-, singly- and doubly-phosphorylated states using non-uniformly sampled NMR experiments with increased dimensionality. Next, we analyze the structural changes induced by the phosphorylation of S40 and SI 9 by PKA and PRAK kinase respectively and the kinetics of these processes. 2. Experimental 2.1. Protein expression and purification The regulatory domain (residues 1-169) of human tyrosine hydroxylase 1 (RD-hTHl) in a pET15b plasmid containing a TEVcleavable His-tag was expressed in E. coli BL21(DE3)R1L cells. For preparation of 1 5 N-labeled and 1 3 C,1 5 N-labeled samples, the cells were cultured in M9 medium with 1 5 NH4 C1, ampicillin (100 mg/ml), chloramphenicol (35 mg/ml) and with addition of 1 3 C6 -glucose for the double labeled sample. Cells grew at 37 °C until OD6oo - 0.8 then were induced with 0.5 mM 1PTG and further cultured at 18 °C for -18 h. Cells were harvested and homogenized in 50 mM Tris pH = 8,150 mM NaCl, 3 mM NaN3 . Cell lysate was centrifuged for 1 h at 21,040g. Supernatant was then applied on a Ni2 + affinity column (HisTrap HP, GE Healthcare) equilibrated in 50 mM Tris pH = 8, 500 mM NaCl, 3 mM NaN3 . Sample was eluted by gradient of elution buffer (equilibration buffer + 1 M imidazole) at its -70% concentration. Sample was then gel filtrated on a Superdex 75 column (HiLoad 16/600 Superdex 75 pg, GE Healthcare) equilibrated in 50 mM Tris pH = 8,100 mM NaCl, 3 mM NaN3 . The eluted sample was treated with TEV protease (proteimprotease ratio 20:1) at 4 °C overnight and then dialysed into phosphate buffer (20 mM sodium phosphate buffer pH = 6, 3 mM NaN3 ). The His-tag cleaved protein was loaded on a cation exchange column (Resource S, GE Healthcare) equilibrated in phosphate buffer (pH = 6.0), and sample was eluted by a gradient of elution buffer (phosphate buffer, pH = 6.0 + 1 M NaCl) at conductivity 23-45 m S - c m - 1 . Fractions containing our protein sample were again dialysed into phosphate buffer (pH = 6.0) and then concentrated for final gel filtration on a Superdex 75 column equilibrated in phosphate buffer. 22. NMR samples and phosphorylation The NMR backbone assignment of the 1DP region (1 -65 aa) was performed on [1 5 N, 1 3 C] labeled samples of non- and doubly phosphorylated RD-hTHl at 1.0 mM concentration in 20 mM sodium phosphate buffer pH = 6 with 8% D2 0 and 3 mM NaN3 . The phosphorylation kinetic studies were performed with [15 N] samples at 0.3 mM concentration in phosphorylation buffer containing 50 mM sodium phosphate buffer pH = 6,10 mM ATP, 10 mM MgCl2 with 8% D2 0 and 3 mM NaN3 . For phosphorylation of S40, PKA (the catalytic subunit of cAMP-dependent protein kinase, New England BioLabs Inc.) was used in 0.5 ug/ml (13.2 nM) concentration. Afterwards, S19 phosphorylation using PRAK (p38 regulated/activated protein kinase, obtained from University Dundee, Scotland) was performed. The concentration of PRAK was 0.11 mg/ml (2.02 (JM). Both phosphorylation reactions were monitored over the course of 40 h. 2.3. NMR experiments The assignment of non-phosphorylated RD-hTHl was performed using a Bruker 850 MHz US2 spectrometer equipped with cryogenic triple-resonance probe head (5 mm CPTC11H/19F-13C/15N/D). The kinetic measurements as well as assignment of phosphorylated RD-hTHl were performed using a Bruker 600 MHz spectrometer equipped with a cryogenic triple-resonance probe head (5 mm CPQC11H-31P/13C/ 15N/D). Both probe heads are equipped with z-axis gradient coils. All measurements were done at a temperature of 293.2 K. For the assignment of non-phosphorylated and doubly S19_S40phosphorylated RD-hTHl, 3D HNCO [12], 5D HN(CA)CONH and 5D HabCabCONH [13,14] were measured. All experiments were carried out with non-uniform sampling of the indirectly detected domains. The time schedule was generated using Poisson disk sampling on a grid, introducing distance constraints between points. The density of points was set according to a Gaussian distribution (a = 0.5). The phosphorylation was monitored using 2D HSQC experiments. In both cases, the maximal evolution times were set to 102 and 64 ms, respectively. The time resolution was achieved using non-uniform sampling in the indirect domain. The sampling schedule comprised 16,000 points in total for both phosphorylations. The size of the schedules exceeded the size of regular Nyquist grid 125 times. The total measuring time was 40 h. 2.4. Data processing The non-uniform 3D HNCO spectra were processed using a Multidimensional Fourier Transform [15], while 5D HN(CA)CONH and 5D HabCabCONH spectra were processed using a Sparse Multidimensional Fourier Transform (SMFT) algorithm [16]. The direct dimension was square cosine weighted and zero-filled to 4096 complex points, followed by a standard FFT. The 5D spectra were processed by the SMFT algorithm using the program reduced, using fixed frequencies of 1 3 C , 1 5 N and 1 H N identified in the 3D HNCO spectrum, providing sets of 2D slices. The assignment and visualization of the NMR spectra was performed in the software Sparky 3.115 [18]. The secondary structure propensities were calculated by the program SSP [19] using chemical shifts of 1 H a , 1 3 C a and 1 3 C P . The data from RefDB [17] were used for random coil referencing. These values were used for calculation of the phosphorylated state as well. The time-resolved HSQC spectra were processed using a coprocessed Multidimensional Decomposition (co-MDD) [11]. The initial window size was set to 64 points. In the case of PRAK phosphorylation, the window size was incremented by a factor of 1.05 to reduce the fitting errors in subsequent analyses. The processing yielded 250 individual frames for PKA and 52 frames for PRAK phosphorylation. P. Lousa et at / Biophysical Chemistry 223 (2017) 25-29 27 3. Results and discussion The 1 5 N - and 1 3 C,1 5 N-labeled samples of RD-hTHl (region 1-169) protein were expressed and purified as described in the Methods section. Final purity of prepared non-, singly and doubly phosphorylated variants was verified by MALD1-TOF-MS spectroscopy (Fig. SI in Suppl. mat.). We want to emphasize that in order to study the 1DP region of RD-hTHl, very high purity of samples is needed because of proteolytic degradation. 3.1. NMR assignment The assignment of the disordered part (1DP, first 65 residues) within the RD-hTHl was performed by employing non-uniform high-dimensional 5D spectra HabCabCONH and HN(CA)CONH. NUS 5D NMR spectroscopy allows measurement of well resolved spectra of the regions with fast tumbling (i.e. slow relaxation) in a very efficient way. The HN(CA)CONH spectrum provided sequential information leading to linkage of several long fragments. This assignment was then verified by amino acid classification based on chemical shifts of 1 H a , 1 H P , 1 3 C°\ 1 3 C P derived from cross-sections of HabCabCONH spectra. Using this approach, we were able to assign all amide resonances in region 1-65 and most resonances of H°\ Hp , C a , Cp , and carbonyl atoms (including prolines) except for those directly preceding prolines, i.e. M l , T3, T8, S31, and V60 and except for Hp , C p atoms of L21 and 142 (Tables SI, S2, S3 in Suppl. mat.). In comparison to previously published NMR assignment of RD of ratTH [7], hereby presented assignment forRD of human TH is more complete, especially in the region around second phosphorylation site (S40). The described assignment procedure was applied to the non- and double-phosphorylated variants of RD-hTHl. The single phosphorylated variant (pS40) was assigned based on the very close similarities in region 1-30 with respect to the non-phosphorylated variant, and region 30-65 was rather similar to the doubly phosphorylated variant. Fig. 1A presents the superposed assigned HSQC spectra of the 1DP region of RD-hTHl in the non-, singly- (pS40) and doubly- (pS19_pS40) phosphorylated states. Naturally, the largest changes in the chemical shifts (Fig. IB) are observed for the phosphorylated residue itself and its neighbors in the primary sequence. The effect on longer distances seems to be more 110 115 120 125 130 A 562 \ 540 1 J R f f i t t £ ^ non-phosphorylated O PS 4 ° Q pS19_pS40 9.0 8.5 7.6 C02 - 1 H (ppm) Q10 A l l K12 G13 F14 R15 R16 A17 V18 S19 E20 L21 D22 A23 K?4 Q25 A26 E27 A28 I29 M30 S31 r R33 F34 I35 G36 R37 R38 Q39 S40 L41 I42 E43 D44 A45 R46 K47 h48 R49 E50 A51 A52 V53 A54 A55 A56 A57 Ab8 A59 r A60 S62 E63 G65 B o.i 0.2 0.3 0.4 0.5 Chemical shift change (ppm) Fig. 1. Superposition of HSQC spectra of differently phosphorylated RD-hTHl (IDP region). The spectra of non-phosphorylated (black), single (pS40) phosphorylated (green) and doubly (pS19_pS40) phosphorylated spectrum (red) are shown in Panel A. The assignment labels are black for peak positions in non-phosphorylated RD-hTHl, green for signals that changed significantly during S40 phosphorylation (by PKA) and red for signals that changed during the subsequent SI 9 phosphorylation by PRAK. Changes of chemical shifts of RD-hTHl observed in HSQC spectra during phosphorylation of S40 (green) and of S19 (red) are shown in Panel B. The asterisks denote the position of phosphorylation sites. 28 P. Lousa et al / Biophysical Chemistry 223 (2017) 25-29 Table 1 Differences of chemical shifts of nuclei used in SSP calculation for residues close to phosphorylation sites. A17 V18 S19 E20 L21 R38 Q39 S40 L41 142 AoCa (ppm) 0.23 0.23 -0.08 -0.51 -0.17 -0.40 0.61 -0.04 -0.26 -0.04 A8CP (ppm) -0.07 -0.19 1.78 0.25 0.14 0.03 -0.19 1.73 0.09 - A8Ha (ppm) -0.03 -0.05 -0.01 0.00 -0.01 0.03 -0.09 -0.03 0.00 -0.02 pronounced around the pS40 site with change observed even at M30. This can be due to the transient secondary structure elements described in the following section. 32. Secondary structural changes induced by phosphorylation The secondary structure propensities (SSP) for non- and doubly phosphorylated RD-hTHl were determined by the approach described in Methods. The phosphoserines were excluded from this analysis due to large change in chemical shift of their beta carbon (Table 1). The SSP values close to zero mean that the particular amino acids are in a random-coil conformation; the larger positive and negative values indicate the alpha-helical or beta-sheet conformations, respectively [19]. Fig. 2A indicates slight alpha-helical propensities in region 17-23 and more pronounced in the region 35-55 for both non- and doubly phosphorylated variants of RD-hTHl. The intensities of signals in 3D HNCO spectra in Fig. 2B further support the existence of alpha-helical structure in the region 35-55. The intensity is related to relaxation properties of residues. The residues with slower tumbling, i.e. inside structured regions, relax faster which manifests as lowered intensity of the resulting signal. The highest intensity is found for the flexible N-terminal residues as expected, while the lowest intensity corresponds to the alpha-helical region suggested by SSP. Both secondary structure propensities (Fig. 2A) as well as the peak intensities (Fig. 2B) over residues within the 1DP region are very similar for non- and doubly phosphorylated RD-hTHl. There is slight increase 10 20 50 60 B phosphorylated RD-hTH l l . RD-hTH 10 20 30 40 50 60 Fig. 2. Secondary structure propensities (A) and intensities (B) of signals in 3D HNCO spectra of non-phosphorylated (black bars) and doubly-phosphorylated (pS19_pS40, grey bars) RD-hTHl (IDP region). The asterisks denote phosphorylation sites. in alpha-helical propensity for region 48-58 and slight decrease of alpha-helical propensity in region 18-22 around the second phosphorylation site (pS19). This observation indicates negligible conformational changes within the RD-hTHl due to the phosphorylation of S19 and S40. This is in agreement with the proposed molecular mechanism of hTHl activation suggesting that the phosphorylation at S40 and S19 induce conformational changes of the whole RD with respect to the catalytic domain rather than within the RD itself [20]. From a methodological point of view, we found it quite surprising that in contrast to amide groups, the chemical shifts of nuclei that are usually used for the SSP analysis ( 1 H a , 1 3 C a and 1 3 C P ) are relatively insensitive to the phosphorylation of serines (Table 1). The standard deviation of chemical shifts for SSP analysis is on the order of 1 ppm, which is much more than the differences in Table 1. The only significant difference is found only for the phosphorylated serine itself and its beta carbon. This insensitivity of chemical shifts of regions around phosphoserines in IDP regions allows us to perform SSP analysis using non-phosphorylated reference values also for other phosphorylated proteins without the need for intrinsic chemical shift referencing [21]. 3.3. Phosphorylation kinetics In Fig. 3 we present the time progress of phosphorylation for involved serine and surrounding well resolved residues. We plot together decrease in intensity of the original peak and increase in intensity of newly formed peak after phosphorylation. First, S40 was phosphorylated using PKA. The progress of phosphorylation was monitored by changes in the intensity of the HSQC spectrum, which provided the information to derive a 0t h order rate constant of k = 0.128 ± 0.003 mM • h"1 (Fig. 3A). After 10 h when S40 was fully phosphorylated, PRAK was added, and phosphorylation of SI 9 was observed. The intensity changes can be described by 1s t order kinetics. Fitting provides a rate constant of k = 0.556 ± 0.008 h"1 (Fig. 3B). After the normalization of rate constants to micromolar kinase concentrations, we obtained these values: 9.7 mM-h_ 1 /|JM for PKA and 0.27 h_ 1 /|JM for PRAK. 4. Conclusions Advanced techniques of non-uniform sampled NMR experiments were applied to gain information about the intrinsically disordered region of the regulatory domain of human tyrosine hydroxylase 1 (RD-hTHl) and its phosphorylated variants. The assignment of non-phosphorylated and doubly phosphorylated (pS19_pS40) samples was performed using non-uniformly sampled 5D NMR spectroscopy. These two sets of chemical shifts were also sufficient to perform assignment of singly phosphorylated (pS40) RD-hTHl. Although the phosphorylation has quite significant impact on amidic chemical shifts of phosphoserine and its neighboring residues, it has negligible effect on aliphatic chemical shifts with the exception of C p of the phosphoserine itself. Therefore, the aliphatic chemical shifts could be employed to monitor alpha-helical propensity in the region between residues 35 and 55 without the need of phosphoprotein reference. Phosphorylation has no significant impact on the secondary structure propensities of RD-hTHl. P. Louša et al / Biophysical Chemistry 223 (2017) 25-29 29 A B time [h] time [h] time [h] time [h] Fig. 3. Progress of phosphorylation by PKA kinase is shown in Panel A. The intensities of signals belonging to non-phosphorylated protein (squares) and S40-phosphorylated (triangles) are plotted against reaction time. Other residues close to S40 were not used due to severe overlap, leading to errors in analysis. Panel B: Progress of phosphorylation by PRAK Intensities belonging to singly (pS40) (squares) and to doubly (pS19_pS40) phosphorylated sample (triangles) are plotted against reaction time. Other residues were strongly influenced by signal overlap and therefore not used in subsequent analysis. The kinetics of phosphorylation of S40 by PKA and of SI 9 by PRAK were determined using non-uniformly sampled time-resolved NMR spectroscopy. The phosphorylation of S40 was determined as a 0t h order reaction with normalized rate constant value 9.7 mM-h_ 1 /|JM, while phosphorylation of SI 9 was observed as 1s t order reaction with value 0.27 h_ 1 /|JM. The described time-resolved NMR spectroscopy approach has general applicability over large time scales of various posttranslational processes and can be easily employed in other protein systems. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.bpc.2017.01.003. Acknowledgment We thank Dr. Tanvir Shaikh for critical reading of the manuscript. The project is financed from the SoMoPro 11 programme. The research leading to this invention has acquired a financial grant from the People Programme (Marie Curie action) of the Seventh Framework Programme of EU according to the REA Grant Agreement No. 291782. The research is further co-financed by the South-Moravian Region. The article/paper reflects only the author's views and the Union is not liable for any use that may be made of the information contained therein. In addition, this work was also supported by the Czech Science Foundation (15- 34684L). The C11SB research infrastructure project LM2015043 funded by MEYS CR is gratefully acknowledged for the partial financial support of the measurements at the Josef Dadok National NMR Centre and Proteolytics Core Facilities, CE1TEC - Masaryk University. References [1] T. Nagatsu, M. Levitt, S. Udenfriend, Tyrosine hydroxylase: the initial step in norephinephrine biosynthesis, J. Biol. Chem. 239 (1964) 2910-2917. [2] P.B. Molinoff, J. Axelrod, Biochemistry of catecholamines, Anna Rev. Biochem. 40 (1971)465-500. [3] S.C. Daubner, D.L. Lohse, P.F. Fitzpatrick, Expression and characterization of catalytic and regulatory domains of rat tyrosine hydroxylase, Protein Sci. 2 (1993) 1452-1460. [4] P.F. Fitzpatrick, Tetrahydropterin-dependent amino acid hydroxylases, Anna Rev. Biochem 68 (1999) 355-381. [5] J. Hritz, I.-J. Byeon, T. Krzysiak, A. Martinez, V. Sklenár, A.M. Gronenborn, Dissection of binding between a phosphorylated tyrosine hydroxylase peptide and 14-3-3zeta: a complex story elucidated by NMR, Biophys. J. 107 (2014) 2185-2194. [6] R. Kleppe, K. Toska, J. Haavik, Interaction of phosphorylated tyrosine hydroxylase with 14-3-3 proteins: evidence for a phosphoserine 40-dependent association, J. Neurochem. 77 (2001) 1097-1107. [7] S. Zhang, T. Huang, U. Ilangovan, AP. Hinck, P.F. Fitzpatrick, The solution structure of the regulatory domain of tyrosine hydroxylase, J. Mol. Biol. 426 (2014) 1483-1497. [8] K. Toska, R. Kleppe, C.G. Armstrong, N.A Morrice, P. Cohen, J. Haavik, Regulation of tyrosine hydroxylase by stress-activated protein kinases, J. Neurochem. 83 (2002) 775-783. [9] FX Theillet, H.M. Rose, S. Liokatis, A Binolfi, R. Thongwichian, M. Stuiver, P. Selenko, Site-specific NMR mapping and time-resolved monitoring of serine and threonine phosphorylation in reconstituted kinase reactions and mammalian cell extracts, Nat. Protoc. 8 (2013) 1416-1432. [10] I. Landrieu, L. Lacosse, A Leroy, J.-M. Wieruszeski, X. Trivelli, A. Sillen, N. Sibille, H. Schwalbe, K. Saxena, T. Langer, G. Lippens, NMR analysis of a tau phosphorylation pattern, J. Am Chem Soc. 128 (2006) 3575-3583. [11] M. Mayzel, J. Rosenlôw, L. Isaksson, V.Y. Orekhov, Time-resolved multidimensional NMR with non-uniform sampling, J. Biomol. NMR 58 (2014) 129-139. [12] L.E. Kay, M. Ikura, R. Tschudin, A. Bax, Three-dimensional triple-resonance NMR spectroscopy of isotopically enriched proteins, J. Magn. Reson. 89 (1990) 496-514. [13] K. Kazimierczuk, A Zawadzka-Kazimierczuk, W. Kozminski, Non-uniform frequency domain for optimal exploitation of non-uniform sampling, J. Magn. Reson. 205 (2010) 286-292. [14] V. Motackova, J. Novacek, A Zawadzka-Kazimierczuk, K. Kazimierczuk L. Zidek, H. Sanderova, L. Krásny, W. Kozminski, V. Sldenar, Strategy for complete NMR assignment of disordered proteins with highly repetitive sequences based on resolutionenhanced 5D experiments, J. Biomol. NMR 48 (2010) 169-177. [15] J. Stanek, W. Kozminski, Iterative algorithm of discrete Fourier transform for processing randomly sampled NMR data sets, J. Biomol. NMR 47 (2010) 65-77. [16] K. Kazimierczuk, A Zawadzka, W. Kozminski, Narrow peaks and high dimensionalities: exploiting the advantages of random sampling, J. Magn. Reson. 197 (2009) 219-228. [17] H. Zhang, S. Neal, D.S. Wishart, RefDB: a database of uniformly referenced protein chemical shifts, J. Biomol. NMR 25 (2003) 173-195. [18] T.D. Goddard and D. G. Kneller, SPARKY 3, University of California, San Francisco, USA. [19] JA Marsh, V.K. Singh, Z Jia, J.D. Forman-Kay, Sensitivity of secondary structure propensities to sequence differences between a- and 7-synudein: implications for fibrillation, Protein Sci. 15 (2006) 2795-2804. [20] S.C. Daubner, T. Lee, S. Wang Tyrosine hydroxylase and regulation of dopamine synthesis, Arch. Biochem. Biophys. 508 (2011) 1-12. [21 ] K. Modig, V.W.Júrgensen, K. Lindorff-Larsen, W. Fieber, H. Bohr, F.M. Poulsen, Detection of initiation sites in protein folding of the four helix bundle ACBP by chemical shift analysis, FEBS Lett. 581 (25) (2007) 4965-4971. Paper 3 Zapletal V., Mládek A., Melková K., Louša P., Nomilner E., Jaseňáková Z., Kubáň V., Makovická M., Laníková A., Zídek L., Hritz J.* * Corresponding author Choice of force field for proteins containing structured and intrinsically disordered regions Biophysical Journal 2020, 118, 7, 1621­1633 doi: 10.1016/j.bpj.2020.02.019 PL performed, processed and analysed the NMR experiments of RD­hTHl, wrote and reviewed the manuscript. 103 BiophysicalJournal A r t i c l e Choice of Force Field for Proteins Containing Structured and Intrinsically Disordered Regions Vojtěch Zapletal,1,2 Arnošt Mládek,2 Kateřina Melková,1,2 Petr Louša,2 Erik Nomilner,1 Zuzana Jaseňáková,1,2 Vojtěch Kubáň,1,2 Markéta Makovická,1 Alice Laníková,1 Lukáš Žídek,1,2 and Jozef Hritz2 * 'National Centre for Bi omolecular Research, Faculty of Sci ence and 2 Central European Institute of Technology, Masaryk Uni versi ty, Brno, Czech Republi c ABSTRACT Bio mo lecular fo rce fields o ptimized fo r glo bular pro teins fail to pro perly repro duce pro perties o f intrinsically disordered proteins. In particular, parameters of the water model need to be modified to improve applicability of the force fields to both ordered and disordered proteins. Here, we compared performance of force fields reco mmended for intrinsically disordered proteins in molecular dynamics simulatio ns o f three proteins differing in the content of ordered and disordered regio ns (two proteins consisting o f a well-structured domain and of a disordered regio n with and without a transient helical mo tif and one disordered protein containing a region of increased helical propensity). The obtained molecular dynamics trajectories were used to predict measurable parameters, including radii of gyration o f the pro teins and chemical shifts, residual dipo lar couplings, paramagnetic relaxatio n enhancement, and NMR relaxatio n data of their individual residues. The predicted quantities were compared with experimental data obtained within this study o r published previo usly. The results sho wed that the NMR relaxation parameters, rarely used for benchmarking, are particularly sensitive to the choice o f force-field parameters, especially those defining the water model. Interestingly, the TIP3P water model, leading to an artificial structural collapse, also resulted in unrealistic relaxatio n properties. The TIP4P-D water model, combined with three biomolecular force-field parameters for the protein part, significantly impro ved reliability of the simulations. Additio nal analysis revealed only one particular force field capable o f retaining the transient helical mo tif observed in NMR experiments. The benchmarking pro to co l used in our study, being more sensitive to imperfections than the commonly used tests, is well suited to evaluate the performance of newly developed force fields. SIGNIFICANCE We compared the performance of several force fields in molecular dynamics simulatio ns of three proteins differing in the content of ordered and disordered regions. Fro m the obtained trajectories, we predicted a set of measurable quantities and compared them with their experimental values. Amo ng the predicted parameters, NMR relaxation data were particularly sensitive to the choice of force field parameters, especially those defining the water model. The presented benchmarking pro to co l will help to select force fields that reliably simulate properties of physio lo gically important intrinsically disordered proteins. INTRODUCTION The fact that many proteins of biological relevance contain considerably large intrinsically disordered regions (IDRs), contradicting the classical structure­function paradigm, has been accepted by the structural biology community during the past two decades (I­3). Intrinsically disordered proteins (IDPs) represent a relatively diverse class of molecules differing in their biophysical properties. One polySubmitted April 17, 2019, and acceptedfor publication February 5, 2020. *Correspondence: jozef.hritz@ceitec.muni.cz Editor: Rohit Pappu. https://doi.Org/10.1016/j.bpj.2020.02.019 © 2020 Biophysical Society. peptide chain often contains fully structured domains together with IDRs. M oreover, IDRs are not random polymer chains but exhibit various degree of partial ordering. They typically contain a short segment with an increased propensity to form secondary transient structures, described in the literature as prestructured motifs (4), preformed structural elements (5), or molecular recognition features (6). Such hybrid systems present methodological challenges because an IDR tethered to a well­ordered domain is a molecule consisting of regions with highly diverse dynamics. As a result, it is difficult for experimental and computational methods to accurately capture both types of behavior in these systems. Biophysical Journal 118, 1621-1633, Apri l 7, 2020 1621 Zapletal et al. NMR represents a method of choice for studies of proteins with IDRs at atomic resolution, as disordered systems are difficult to investigate using single-crystal x-ray diffraction or single-particle reconstruction of cryo-electron microscopic images. Currently available NMR methods provide sufficient resolution to overcome the narrow distribution of chemical shifts of IDPs (7,8). However, it should be noted that the resolution improvement of IDP-targeted NMR experiments relies on the slow relaxation of IDRs. Therefore, the sensitivity of such experiments is often too low for the rapidly relaxing signals of amino acids in the well-ordered regions of hybrid proteins. Molecular dynamics (MD) simulations could, in principle, serve as an ideal tool to study behavior of hybrid proteins at an atomic level. Moreover, most of the NMR parameters can be reliably predicted from structural models, which allows for direct comparison of MD results with experimental data. In practice, several problems complicate M D simulations of IDPs. The energy landscapes of IDRs are expected to be weakly funneled such that the search for a specific functionally competent conformation could be extremely inefficient (9). In this study, we examined the applicability of MD simulations to hybrid proteins and assessed their reliability by predicting several measurable parameters from the obtained trajectories and comparing them with experimental data. Moreover, we provide an example of prediction of measurable parameters as a guide to select optimal setup for future experiments, such as positions with paramagnetic relaxation enhancement (PRE) labels. We tested the currently available force fields Amber99SB-1LDN (A99), CHARMM22* (C22*), and CHARMM36m (C36m) in combination with explicit solvent models T1P3P, T1PS3P, and T1P4P-D. Chemical shift, residual dipolar coupling (RDC), PRE, relaxation rate, and smallangle x-ray scattering (SAXS) experimental data were used for validation, ft should be pointed out that our goal was not to describe the properties of the studied proteins as faithfully as possible but to look for features that would distinguish the performance of various force fields in MD simulations of a microsecond time range. The hybrid proteins investigated in this study included 1) 5 subunit of RNA polymerase from Bacillus subtilis (<5RNAP), 2) regulatory domain of human tyrosine hydroxylase (RD-hTH), and 3) a fragment consisting of residues 159-254 of rat microtubule-associated protein 2c (MAP2c1 5 9 - 2 5 4 ). 5RNAP makes RNA polymerase sensitive to the concentration of initiating nucleoside triphosphates, which is important for rapid changes in gene expression (10). ft consists of two domains of a similar size. The N-terminal half of <5RNAP folds into a well-ordered, mostly a-helical domain, whereas the C-terminal half is disordered and highly negatively charged, with the exception of a lysine-rich motif 9 6 K A K K K K A K K 1 0 4 and C-terminal Lysl73 (11). Experimental data showed that the lysine stretch makes transient electrostatic contacts with various residues in the acidic C-terminal region (1). No sign of formation of transient a-helical structures was observed in the C-terminal domain, which prefers extended backbone conformations, presumably because of the electrostatic repulsion of the acidic side chains. RD-hTH catalyzes hydroxylation of L-tyrosine to L-3,4-dihydroxyphenylalanine (L-DOPA) and is a key and rate-limiting enzyme in biosynthesis of important catecholamine neurotransmitters (12,13). Its N-terminal region (~40% of its sequence) is disordered but contains a segment with ~80% propensity to form four turns of a-helix (14). M A P 2 c 1 5 9 - 2 5 4 corresponds to the central region of MAP2c, where proteins regulating the microtubule-stabilizing activity of MAP2c bind in a phosphorylation-dependent manner (15,16). M A P 2 c 1 5 9 " 2 5 4 is mostly disordered but exhibits an ~20% propensity to form four turns of an a-helix presumably important for intermolecular interactions (16,17). MATERIALS AND METHODS NMR spectroscopy NMR assignments of <5RNAP and of the disordered part of RD-hTH were published previously (14,18). M A P 2 c 1 5 9 - 2 5 4 w a s assigned as described for full-length MAP2c (19). PRE data of <5RNAP were published previously (11). RDCs were calculated as a difference between splitting observed in in-phase/anti-phase (IPAP) spectra (20) obtained for proteins in stretched 5% polyacrylamide gel and in isotropic medium. The RD-hTH data were acquired at 20°C on a 600 and 850 MHz Bruker Avance III spectrometer (Bruker, Billerica, MA), and the M A P 2 c 1 5 9 - 2 5 4 data were acquired at 27° C on a 600 MHz Bruker Avance III spectrometer. The backbone amide 1 5 N relaxation data were published previously for 5RNAP (21) and measured using standard pulse sequences (22) for 1.0 m M [1 5 N]-RD-hTH and 0.34 mM [ 1 5 N ] - M A P 2 c 1 5 9 - 2 5 4 at 27°C on 850 and 950 M H z Bruker Avance III spectrometers, respectively. The interscan delays were set to 1.5, 2, 6, and 25 s for R,, R2, and steady-state heteronuclear Overhauser effect (ssNOE) measurements without and with ssNOE transfer, respectively. The relaxation delays for the R\ experiment were 11.2, 16.8, 28.0, 44.8, 61.6, 95.2, 196, and 308 ms, and the delays for the R2 experiment were 0, 14.4, 28.8, 43.2, 57.6, 86.4, and 129.6 ms. Rt and R2 data were fitted using a two-parameter exponential. The errors were obtained using the bootstrap procedure (23). The ssNOE values were calculated as a ratio between signal intensities obtained from spectra with or without ' H saturation. The errors were derived from background noise levels in each individual spectrum. SAXS The SAXS data sets were collected using a BioSAXS-1000 (Rigaku, Tokyo, Japan) instrument with an x-ray beam wavelength of 1.54 A at 27°C. The distance between the sample and the detector (PILATUS 100K; Dectris, Baden-Daettwil, Switzerland) was 0.48 m, covering a scattering vector (q = 4irsin(6)/X) range from 0.009 to 0.65 A ~ . For solvent and sample, one two-dimensional image was collected with 1 h exposure time per image. Radial averaging of two-dimensional scattering images and the solvent subtractions were performed using SAXSLab3.0.0rl (Rigaku). All data sets were truncated to a maximal scattering vector of 0.3 A - 1 for further analysis. Radii of gyration were determined using PRIMUS from ATSAS v2.7.2 (24) with Guinier analysis (implemented in the PRIMUS Guinier Wizard) in the range 2-52 for <5RNAP and in the range 9-30 for M A P 2 c 1 5 9 " 2 5 4 ; the globular particle type was used for both proteins. The molecular form factor analysis (25) was performed online at http://sosnick.uchicago.edu/ SAXSonlDPs. 1622 Biophysical Journal 118, 1621-1633, April 7, 2020 Force Field for Disordered Regions Computational details The M D simulations were performed using Amber99SB-ILDN (26), CHARMM22* (27), and CHARMM36m (28) force-field parameters for the protein atoms and the TIP3P, TIPS3P (29,30), or TIP4P-D (31) water models. The proteins were solvated using a rhombic dodecahedral box of waters with a minimal distance between the box walls and solute of 2 nm. The charge of the system was neutralized by adding C P and N a + ions, and the concentration of salt was adjusted to 100 mM. A l l simulations were performed under periodic boundary conditions. Before the M D runs, in vacuo and solvent energy minimizations with the steepest descent algorithm were carried out. The lengths of bonds with hydrogen atoms were constrained using the LINCS algorithm (32). An integration time step of 2 fs was used. A cutoff of 1.0 nm was applied for the LennardJones interactions and short-range electrostatic interactions. Long-range electrostatic interactions were calculated by particle mesh Ewald summation with a grid spacing of 0.12 nm and a fourth-order interpolation (33). The four-step 8-ns-long equilibration protocol consisted of the following parts: 1) 2-ns relaxation of water molecules at 300 K during the N V T equilibration with restrained (1000 kj m o P 1 nnP ) solute coordinates, 2) 2-ns NVT (300 K) run with restrained (1000 kj m o P 1 nnP2 ) backbone atom coordinates, 3) 2-ns NpT (300 K, 1 atm) run with restrained (1000 kj moP nnP ) backbone atom coordinates, and 4) 2 ns of unrestrained NpT simulation (300 K, 1 atm). The length of the follow-up production NpT simulations (300 K, 1 atm) was 200 ns. The simulations using the TIP4P-D water model were further prolonged to 1 fis, and additional independent simulations (500 ns for (5RNAP and M A P 2 c 1 5 9 " 2 5 4 , 400 ns for RD-hTH) starting from different initial conditions were run for all listed force fields and protein systems. The temperature and pressure were maintained using the Berendsen coupling scheme (34) during the equilibration steps; the production NpT simulations were performed using the velocity rescaling thermostat with a stochastic term (35) and the Parrinello-Rahman barostat algorithms (36). Atomic coordinates were recorded every 1 ps. Predictions of NMR parameters Chemical shifts were calculated using the prediction algorithm SPARTA+ (37) for each structure, and averaged secondary chemical shifts (SCSs) were calculated by subtracting the random-coil values (38). RDC calculations using a local alignment window were performed. For calculations using a local alignment window, the RDC, calculated using the program PALES (39), for the central amino acid of the local 15-amino-acid segment was calculated for each conformer (40). The resulting RDC profile along the primary sequence was calculated by averaging each value over the whole trajectory and multiplying by the corresponding scaled absolute value of the generic baseline to account for long-range effects (41). The scale was chosen so that the lowest root mean-square deviations (RMSDs) from the experimental values were obtained in the disordered regions. PRE was calculated for the spin label used in the experiments, i.e., for thiol-reactive methanethiosulfonate (MTSL) attached to a side chain of cysteine introduced by site-directed mutagenesis. Sterically allowed MTSL sidechain conformations were sampled using previously published rotameric distributions (42) and built explicitly for each spin-label site of each individual structure backbone. 600 side-chain conformers were calculated, and the sterically allowed conformers were retained. Relaxation effects were averaged over these conformers as described by Salmon et al. (41). The strategy described previously (43) was used to calculate relaxation rates. The autocorrelation function C,(T) was calculated from subtrajectories of each simulation, probing two timescales (TM A X < 5 or 50 ns). Each trajectory was thus divided into the blocks of 10 and 100 ns, respectively, and averaged. For each averaged block, C,(T) was described as a sum of N = 512 exponentials e~T /T < whose amplitudes A, were obtained using a Tikhonov regularization procedure (44). For each averaged block, the spectral densities are then defined as L —' 1 -4- (t) T • / — 1 1 t w ' c,l and used to predict spin relaxation rates. RESULTS AND DISCUSSION Impact of selected water model Our first goal was to examine the performance of various models of water in simulations of three hybrid proteins used as test molecules in this study. It is well described in the literature (31) that the water models typically used in MD simulations (e.g., TIP3P) significantly underestimate London dispersion interactions. To prevent this problem, TIPS3P and TIP4P-D water models were also tested, and the reliability of the behavior of IDPs or the IDRs was analyzed. Modified TIP3P containing nonzero van der Waals parameters was introduced originally in 1998 (30) and shown to provide more realistic results for IDPs when combined with C36m (28). Simulations using this model typically result in extended states when other models of water tend to produce ensembles that are structurally too compact relative to experiments (31). Our preliminary 200-ns MD runs, using the TIP3P water model in combination with A99 (26), confirmed the artificial behavior of TIP3P. The inferior performance of TIP3P was manifested most clearly by the calculated radius of gyration (Rg) of M A P 2 c 1 5 9 - 2 5 4 . All simulations started with conformations having Rg close to the experimental value of 2.5 nm obtained from SAXS. During the initial 100 ns of simulations with TIP3P, Rg dropped to ~1.5 nm regardless of the protein force field used (Fig. 1 b). This result confirms that attractive interactions leading to formation of compact structures are unnaturally enhanced when TIP3P is used, as was also reported previously in (31,45). We also ran 200-ns simulations with the TIPS3P model (30) in combination with the C22* (26) and C36m (28) force fields. Remarkably, TIPS3P did not prevent the artificial compaction of M A P 2 c 1 5 9 - 2 5 4 (Fig. 1 b; a comparison of trajectories calculated using C36m with TIP3P and TIPS3P is presented in Fig. SI). A similar trend was observed also for RD-hTH (Fig. 1 c), although a direct comparison with the experiment was not possible because of the RD-hTH dimerization in real samples. Finally, we performed MD simulations with the TIP4P-D water model combined with A99, C22*, and C36m. In agreement with the literature (31) and in a sharp contrast with the simulations run with TIP3P and TIPS3P, Rg of M A P 2 c 1 5 9 - 2 5 4 did not drop below the value calculated from the SAXS data (Fig. 1 e). In the case of RD-hTH, TIP4P-D also prevented unexpected compaction but only in combination with C36m (see the discussion of the impact of the protein force field parameters in the following section). Biophysical Journal 118, 1621-1633, April 7, 2020 1623 Zapletal et al. 1.0 L= I I l= I I I I I I I I L _ 0 100 200 0 100 200 300 400 500 600 700 800 900 Time / ns Time / ns FIGURE 1 Simulated Rg of (5RNAP (a and d), M A P 2 c 1 5 9 - 2 5 4 (b and e), and RD-hTH (c and/) obtained using the TIP3P (TIPS3P in the case of C22* and C36m) water model (a-c) and TIP4P-D water model (d-f). Experimental data are shown in gold, and values calculated using A99, C22*, and C36m are shown in green, red, and blue, respectively. To see this figure in color, go online. Interestingly, the water model influenced simulated Rg of 5RNAP less than Rg of M A P 2 c 1 5 9 " 2 5 4 and RD-hTH. We observed the artificial collapsed structure only rarely in simulations of 5RNAP with TIP3P or TIPS3P (see A99 data in Fig. 1 a). The most likely explanation is that the C-terminal IDR of 5RNAP is very strongly negatively charged, preventing its collapsed structure regardless of applied water model (Fig. S2). To examine the influence of water models in the simulations of <5RNAP more closely, we also calculated NMR relaxation parameters from the M D trajectories. As discussed below in more details, the prediction of relaxation parameters is a challenging task. Therefore, we hoped that a comparison of calculated and experimental relaxation parameters might reveal more subtle effects. Indeed, we observed much lower predicted values of ssNOE for residues 134-173 of 5RNAP (Fig. 2 a). Such values indicate an artificially high disorder of the highly acidic region further than 30 residues from the positively charged lysine stretch (residues 96-104). Importantly, this effect was observed not only for the original TIP3P model but also for the modified version TIPS3P (red and blue traces in Fig. 2 a). Limitations of TIPS3P were already noticed by Huang et al., who reported that this model with C36m provided correct Rg for the arginine-serine but underestimated the Rg of the Thermotoga maritima cold-shock protein and of the N-terminal domain of HIV-1 integrase (28). A further modification of TIPS3P improved the Rg prediction for the 100 120 140 160 100 120 140 160 Residue number Residue number FIGURE 2 ssNOE (a and d) and relaxation rates Tx (b and e) and R\ (c and/) of the C-terminal IDR of <5RNAP obtained with the 5-ns sliding window using the TIP3P (TIPS3P in the case of C22* and C36m) water model (a-c) and TIP4P-D water model (d-f). Experimental data are shown in gold, and values calculated using A99, C22*, and C36m are shown in green, red, and blue, respectively. To see this figure in color, go online. cold-shock proteins but worsened the agreement for the other two proteins (28). In conclusion, significant differences between water models were already observed during the first 200 ns of the simulations. Considering the results for TIP3P, TIPS3P, and TIP4P-D, we continued our study only with the TIP4PD model, extending the MD simulations to the microsecond range. Impact of force fields on global shape of proteins We first examined the global shape of studied proteins in terms of Rg and SAXS curves. The /?g-values calculated from 1-fis (Fig. 1) and 0.5-^s (Figs. S3 and S4) A99, C22*, and C36m trajectories of <5RNAP (simulated using TIP4P-D) oscillated between 3 and 6 nm. The experimental f?g-value, obtained by the Guinier analysis of the SAXS curve, was (3.45 ± 0.3) nm. The Rg alone did not reveal any systematic difference among the force fields, except for an extending event observed around 550 ns in the C36m simulation. The /?„-values calculated from M D trajectories of MAP2c 2 5 4 oscillated between 2 and 4.5 nm, close to the experimental value of (2.5 ± 0.3) nm (Fig. 1 e). Figs. S5 and S6 show the comparisons between the predicted and measured SAXS profiles for <5RNAP and M A P 2 c 1 5 9 " 2 5 4 and for individual force fields. The lowest Q-factors, indicating the best fit (Table 1), were obtained 1624 Biophysical Journal 118, 1621-1633, April 7, 2020 Force Field for Disordered Regions TABLE 1 Comparison of Calculated RMSD from Experimental Values and Normalized scores Metric <5RNAP RD-hTH M A P 2 c 1 5 9 - 2 5 4 A99 C22* C36m A99 C22* C36m A99 C22* C36m RMSD: zl<5C7ppm 0.92 0.88 0.82 0.92 0.84 0.64 0.74 0.44 0.55 0.86 0.78 0.79 0.63 0.52 0.43 0.45 0.45 0.41 Zl<5C(0)/ppm 0.74 0.67 0.62 0.96 0.90 0.56 0.84 0.55 0.65 zl<5N/ppm 1.71 1.92 1.67 3.55 3.40 2.88 1.54 1.90 1.07 Q for £>(NHN ) 0.30 0.33 0.33 1.21 1.12 0.89 0.76 0.69 0.74 Local PRE'' by MTSL at L I IOC 0.25 0.41 0.26 n.d.b n.d. n.d. n.d. n.d. n.d. Local PRE by MTSL at L132C 0.07 0.17 0.16 n.d. n.d. n.d. n.d. n.d. n.d. Local PRE by MTSL at L151C 0.12 0.13 0.14 n.d. n.d. n.d. n.d. n.d. n.d. Local PRE by MTSL at L168C 0.06 0.06 0.07 n.d. n.d. n.d. n.d. n.d. n.d. PRE (all data) 0.08 0.10 0.09 n.d. n.d. n.d. n.d. n.d. n.d. ssNOE 0.15 0.17 0.20 0.11 0.38 0.25 0.15 0.20 0.16 R2ot j y s - 1 1.50 1.56 2.46 2.62 2.32 0.93 1.13 1.32 1.09 0.12 0.15 0.14 0.19 0.24 0.13 0.21 0.19 0.23 Q for SAXS (0.1 nrrT1 < q < 1 nrrT1 ) 0.057 0.040 0.084 n.d. n.d. n.d. 0.070 0.048 0.075 Q for SAXS (all data) 0.060 0.045 0.085 n.d. n.d. n.d. 0.084 0.066 0.083 Score: Jcs 1.11 1.08 1.00 1.46 1.32 1.00 1.42 1.22 1.11 S RDC 1.00 1.09 1.08 1.33 1.20 1.00 1.06 1.00 1.36 s P R E (all data) 1.00 1.25 1.11 n.d. n.d. n.d. n.d. n.d. n.d. ^relax 1.00 1.13 1.37 1.74 2.57 1.44 1.04 1.16 1.07 % M R 1.00 1.13 1.44 1.54 1.89 1.22 1.06 1.12 1.07 SSAXS (all data) 1.32 1.00 1.87 n.d. n.d. n.d. 1.28 1.00 1.26 Sail 1.08 1.08 1.21 1.55 1.78 1.17 1.24 1.15 1.11 Rg/nm 4.13 4.03 4.66 n.d. n.d. n.d. 2.60 2.94 2.78 Rg penalty 0.13 0.10 0.28 n.d. n.d. n.d. 0 0.06 0 ^combined 1.21 1.18 1.49 1.55 1.78 1.17 1.24 1.21 1.11 "Calculated for PRE of residues 84-104 because of the indicated spin label. Input experimental data not determined. for the C22* force field. Inspection of individual SAXS profiles revealed that C36m and A99 overestimated SAXS intensities for M A P 2 c 1 5 9 " 2 5 4 and underestimated SAXS intensities for <5RNAP in the medium- and low-q region (0.1 nm"1 < q < 1 run- 1 ; see Table 1), whereas C22* predicted the medium- and low-q SAXS intensities well. We also fitted the experimental and predicted SAXS profiles to molecular form factors (MFFs) developed by Riback et al. (25). Fitting the experimental M A P 2 c 1 5 9 " 2 5 4 data to an MFF provided Rg = (2.87 ± 0.3) nm and a value of the Flory exponent typical for a fully unfolded polypeptide (v = 0.59 ± 0.02). Very similar values were obtained by fitting the profile simulated by C22*, Rg = (2.952 ± 0.003) nm and v = 0.603 ± 0.001. As expected, the MFF did not fit well the SAXS profile of 5RNAP, containing a large well-ordered domain. The effect of a force field on the global conformations of RD-hTH was examined as well (Fig. 1 f). Compact structures with Rg ~ 2 nm were formed after 500 ns in simulations with A99 and C22* force fields but not with C36m. Therefore, it seems that CH36m in combination with TIP4P-D most efficiently prevents the collapse of structures of RD-hTH and MAP2c, representing proteins not exhibiting extraordinary electrostatic repulsion. Impact of force fields on long-range contacts After evaluating the effect of different force-field parameters on the global shapes of the studied proteins, we compared the ability of the force fields to properly reproduce contacts between residues further apart in the sequence. We started by inspecting <5RNAP as a test system for which long-range contacts are observed experimentally, yet no transient helicity is observed in the C-terminal IDR. Residue pairwise distance maps represent an efficient way to evaluate contacts formed during individual M D runs. Rectangular boxes in the maps presented in Fig. 3 highlight distances 1) between the lysine-rich stretch 9 6 K A K K K K A K K 1 0 4 and highly acidic residues in the C-terminal region and 2) between the ordered N-terminal and disordered C-terminal domains. Differences between distances represented by colors document that the relative average distances varied depending on the force field used. For example, the lysine tract interacted with the residues in the vicinity of Glul20 more strongly in the A99 simulation than in the runs with the C22* and C36m force fields (blue rectangles in Fig. 3). The average distances between the N-terminal domain and vicinity of Glul20 (red rectangles in Fig. 3) also differed. We want to emphasize that the reliability of distance maps is influenced not only by the capabilities of individual force fields applied for Biophysical Journal 118, 1621-1633, April 7, 2020 1625 Zapletal et al. 10 20 30 40 50 60 7C SO 90 100110120130140150160170 10 20 30 40 50 60 70 80 90 100113120130140150160170 10 20 30 40 50 60 70 80 90 10C110120130140150160170 b d f 10 20 30 40 50 60 7C 80 90 100110120130140150160170 10 20 30 40 50 60 70 80 90 100110120130140150160170 10 20 30 40 50 60 70 80 90 10C110120130140150160170 Residue number Residue number Residue number FIGURE 3 Maps describing the distances between residue pairs of <5RNAP simulated using A99 (a and b), C22* (c and d), and C36m (e and/) force fields in combination with TIP4P-D. The scales shown at the right indicate color coding of mean inter-residue C"-C" distances (lower row) and of populations of events when the distance between any pair of residue atoms was shorter than 0.4 nm (upper row). To see this figure in color, go online. the 5RNAP in combination with the TIP4P-D water model but also by limited sampling within the trajectories of the cumulative length of 1.5 ,us. To directly compare the contacts observed during the MD calculations with the experimental data, we simulated PRE of individual residues of <5RNAP for a series of spinlabel positions examined in a previous experimental study (11). The results showed that A99 realistically reproduced experimentally observed contacts of the lysine stretch 9 6 K A K K K K A K K 1 0 4 with labels placed at LI 10C and L132C in the highly acidic C-terminal sequence (Fig. 4 , a and b). The experimentally detected contact with L151C was also predicted. In agreement with the distance map, C22* overestimated contacts of the lysines with the closest label at LI 10C and underestimated contacts with more distant labels, including L132C. C36m did not overestimate contacts with the label at LI 10C but underestimated contacts with the label at L132C similarly to C22* and A99 (see RMSD for PRE in the lysine stretch for individual spin labels, listed in Table 1). In conclusion, correct prediction of electrostatic contacts between residues apart in the sequence is a challenging task. In the case of <5RNAP, A99 with TIP4P-D performed best, being able to predict reliably distances between amino acids separated by up to 50 residues in the sequence. Impact of protein force-field parameters on local conformations In the next step, we analyzed the accuracy of description of local backbone conformations of <5RNAP in the simulations. For this purpose, we compared experimental values of several NMR parameters reflecting the local backbone conformation with the values calculated from snapshots of the M D simulations. First, we checked values of RDC. In principle, RDCs depend both on local conformation and on the overall shape of the molecule, determining the distribution of orientations of the molecule in a partially aligned environment (20). However, it is very demanding to achieve a good sampling of conformations and orientations to faithfully reproduce experimental data without any prior information. Therefore, we used the knowledge of long-range electrostatic contacts in the <5RNAP molecule, obtained experimentally as PRE, and applied the local averaging window ( 4 0 ) to predict RDC values. The predicted RDC values are plotted in Fig. 5 a. The prediction varied especially in the vicinity of the lysine stretch, in 1626 Biophysical Journal 118, 1621-1633, April 7, 2020 Force Field for Disordered Regions 70 80 90 100 110 ' Residue number 130 140 1 I 160170 FIGURE 4 Simulated PRE of <5RNAP with the spin label at L I 10C (a), L132C (b), L151C (c), and E168C (d). Experimental data are shown in gold, and values calculated from M D simulations using the TIP4P-D water model combined with the A99, C22*, and C36m force fields are shown in green, red, and blue, respectively. To see this figure in color, go online. which A99 and C36m achieved better agreement with the experiment than C22*. The second examined NMR parameter reflecting local conformation was the chemical shift. In comparison with RDC, the chemical shifts are measured with higher precision, and their prediction does not require a prior knowledge of molecular orientation. SCSs (deviations of predicted chemical shifts from their random-coil values (38)) are compared with the experimental data in Fig. 5, b-e. Values provided by all force fields agree well with the experimental chemical shifts, except for some mismatch of data predicted by C22* for the lysine stretch. In summary, all force fields predicted local conformation of the extended disordered region of <5RNA reasonably well. The slightly worse prediction in the lysine-rich motif by C22* most likely reflects less accurate description of electrostatic contacts by the force field. Simulation of transient helical regions In the next step, we tested how different force fields describe transient a-helical elements in (partially) disordered proteins. A propensity to adopt a-helical secondary structure represents another level of complexity of IDP conformations not present in the mostly extended C-terminal domain of 5RNAR To explore its effect on the simulations, we inspected M D trajectories of RD-hTH and M A P 2 c 1 5 9 - 2 5 4 These proteins were chosen so that they differ in a-helical propensity. Experimental values of chemical shifts (17,19) indicate that RD-hTH and M A P 2 c 1 5 9 - 2 5 4 form a-helices in the regions 40-53 and 200-216 with ~80 and 20% propensity, respectively (cf. Fig. 5). We performed a set of simulations with different force fields, starting from structures containing ideal a-helices in the experimentally identified a-helical regions, and observed the stability of the helices (Figs. S10S12). During M D simulations using A99 and C22* force fields, the initially present a-helix unfolded in less than 80 ns (Fig. S12). Huang et al. (28) reported that C36m optimized with TIPS3P 1) correctly simulated transient a-helices and 2) provided correct Rg for the arginine-serine peptide but not for other two tested IDPs. Therefore, we were curious whether C36m would keep its ability to maintain the transient a-helices with TIP4P-D. For RD-hTH, C36m with TIP4P-D maintained the a-helical conformation in the whole l-/xs trajectory and for 220 ns in an independent 400-ns run (Fig. SI3). In the case of M A P 2 c 1 5 9 - 2 5 4 , the transient a-helix unfolded after ~35, 85, and 830 ns in three independent runs starting from conformations including the helix at the beginning (Fig. S9). The lower stability of the M A P 2 c 1 5 9 - 2 5 4 helix in simulations is in agreement with its low (20%) population observed experimentally. Formation of the well-defined helix was not observed in the remaining 465, 415, and 170 ns of the trajectories or during two 0.5-,u,s runs started from conformations without the helix. However, temporary formation of the helix (for ~50 ns) was observed in simulations of M A P 2 c 1 5 9 " 2 5 4 using C22* (Fig. S12 a) and A99 (Fig. S12 d). It should be emphasized that the calculated trajectories do not fully sample the equilibrium canonical ensembles. Therefore, we do not expect to observe quantitative agreement between the experimental and simulated populations of transient a helices. The ability of the force fields to reliably describe transient a-helices was directly reflected by the predicted SCS values (Fig. 5). Outside of the helical region, a good agreement of the predicted and experimental chemical shifts was obtained for all force fields tested. Deviations from the experimental values were observed for the A99 and C22* simulations in the region where the originally present a-helix unfolded. For the C36m simulations, SCSs typical for a-helices were obtained in the regions where the a-helix was modeled (Fig. 5, g-i and l-n). Quantitatively, the values of predicted versus SCSs corresponded to 87% populations of the helix in the simulations vs. ~80% population in the real sample of RD-hTH. In the case of M A P 2 c 1 5 9 - 2 5 4 , data from all trajectories were combined in a ratio corresponding to the experimentally estimated 20% population of the helix. In conclusion, the ability of C36m and TIP4P-D to keep the transient a-helices during the simulation and to prevent the artificial collapse of IDR structures suggests that this combination is most robust for our studied proteins. This is noteworthy considering that C36m was not optimized in combination with TIP4P-D, but with TIPS3P and Biophysical Journal 118, 1621-1633, April 7, 2020 1627 Zapletal et al. 5RNAP 1 MAP2C 1 5 9 - 2 5 4 : 159 m 120 140 Residue number 20 40 Residue number 200 220 Residue number FIGURE 5 Values of RDC (a,/, and k) and SCS (b-e, g-j, and l-o) in C-terminal IDR (residues 85-173) of <5RNAP (a-e), N-terminal IDR (residues 1-65) of RD-hTH (f—f), and M A P 2 c 1 5 9 - 2 5 4 (k-o) simulated with the TIP4P-D water model. Experimental data are shown in gold, and values calculated using A99, C22*, and C36m are shown in green, red, and blue, respectively. The random-coil limits are shown in gray. The M D simulations started from structures with a-helices modeled for residues 40-53 of RD-hTH and 200-216 of M A P 2 c 1 5 9 - 2 5 4 . For the sake of clarity, the statistical errors are not displayed here but separately in Figs. S7-S9. To see this figure in color, go online. its variants (28). We believe that the C36m and TIP4P-D combination is a promising candidate for a benchmarking on a wider range of proteins, as was done, e.g., by Robustelli et al. (45). It should be noted that the tests discussed so far utilize a limited set of the most readily available experimental values. We already discussed that NMR relaxation data distinguished the performance of TIP3P and TIP4P-D water models better than Rg in the A99 and C36m simulations of 5RNAP (details are presented in the section Impact of Selected Water Model). In the next section, we examine whether extending the benchmarking to NMR relaxation helps to discriminate the ability of force fields to describe dynamic properties of hybrid proteins. Simulation of NMR relaxation data To test the examined force fields within the MD framework more thoroughly, we compared their abilities to reproduce NMR relaxation rates. The NMR signal relaxes because of the local magnetic fields that fluctuate as a result of stochastic molecular motions. The stochastic reorientation of vectors describing the interactions contributing to relaxation is described by the correlation function. In real samples, large numbers of molecules with an almost isotropic distribution of orientations are measured. Consequently, the correlation functions have a simple analytical form (series of exponential functions). In simulations, the sufficient sampling of the orientations is difficult to achieve (in principle, data of large sets of independent trajectories should be averaged), which deteriorates the calculated correlation function regardless of the force field employed. In our analyses of limited numbers of trajectories, the ensemble averaging was approximated by calculating averages of correlation functions of different regions of the trajectories. The simulations should be also sufficiently long because the correlation function reflects only the effects of motions on the timescale covered by the simulation. The obtained predictions must be interpreted carefully to not confuse the artifacts of the force fields with the effects of the sampling scheme used. The 5RNAP molecule is particularly well suited to test the ability of force fields to predict relaxation rates because its domains greatly differ in their stochastic motions. Dynamics of the well-ordered N-terminal domain is dominated by overall tumbling and probes the ability of the force fields to reproduce hydrodynamic properties. In contrast, internal motions contribute most significantly to the relaxation rates 1628 Biophysical Journal 118, 1621-1633, April 7, 2020 FIGURE 6 ssNOE (a, d, and g) and relaxation rates Tx (b), R2 (e and h), and R, (c,f, and i) of (5RNAP (a-c), N-terminal IDR (residues 1-65) of RD-hTH (d-f), and M A P 2 c 1 5 9 " 2 5 4 (g-i). Experimental data are shown in gold, and values calculated using A99, C22*, and C36m are shown in green, red, and blue, respectively. Solid and dashed lines indicate data obtained from correlation functions calculated using 5- and 50-ns windows, respectively. For the sake of clarity, the statistical errors are not displayed here but separately in Figs. S13-S15. To see this figure in color, go online. of residues in the disordered C-terminal domain. Analysis of the simulated trajectories showed that correlation functions calculated for 5-ns sliding time windows (solid lines in Fig. 6, a-c) describe relaxation in the C-terminal IDR sufficiently well but fail to match the experimental data in the well-ordered region, where slower motions, most notably the overall tumbling, dominate the dynamics. Smooth profiles of the calculated NMR relaxation parameters, resembling the experimental profiles in the IDR, document that the use of a short sliding window allowed us to average a sufficient number of correlation functions and to capture most important modes of motion. To obtain relaxation rates close to the experimental values also in the well-ordered N-terminal region, the time window was extended to 50 ns, and averages of 12 correlation functions were calculated (dashed lines in Fig. 6, a-c). Prediction of the cross-correlated Tx rate, which is most sensitive to slow motions, was most informative (Tx was used to monitor the slow motions instead of the more frequently measured R2 rates because I'x could be obtained with a higher accuracy than R2 for <5RNAP, which has an extremely poor dispersion of chemical shifts in its C-terminal IDR (21)). Predicted and experimental values were comparable for tested force fields, albeit scattered for individual residues because of the contribution of slow motions that were not sufficiently averaged in the small set of the 30 independent correlation functions. The match of the experimental and (average) simulated /^-values is in the line with the fact that T1P4P-D reproduces the water diffusion coefficient more reliably than the T1P3P model (31). The general agreement was also good in the disordered region for all three force field, but significant differences were observed when different regions of the <5RNAP sequence were compared. The trend of ssNOE values, most sensitive to fast motions, was best reproduced by A99 (green line in Fig. 6 a). Also, the predictions by C H A R M M force fields reflected their abilities to predict long-range contacts. Somewhat higher ssNOE values (and elevated Tx) in the vicinity of residues 100 and 125 indicated that C36m slightly overestimated partial ordering in the most rigid regions of the C-terminal domain, where longrange contacts were predicted (see the section Impact of the Force Fields on Long-Range Contacts). This is in agreement with the difficulty of predicting contacts between residues far in the sequence, as discussed above. C22* overestimated ssNOE around residues 115 and 90, in agreement with its tendency to prefer contacts of the lysine stretch with residues closer in the sequence. In conclusion, we calculated NMR relaxation rates from MD trajectories using two time windows (50 and 5 ns) to cover different timescales of motions in the ordered and disordered regions, respectively. The results independently confirmed the observed moderate differences in predicting long-range contacts and showed that all force fields describe well the hydrodynamic properties for the T1P4P-D water model. For <5RNAP with a well-ordered domain and highly charged disordered domains not forming transient helical Biophysical Journal 118, 1621-1633, April 7, 2020 1629 Zapletal et al. structures, A99 performed best. The C36m force field described long-range electrostatic contacts slightly worse, but its accuracy was acceptable. C22* somewhat overestimated electrostatic contacts between residues close in the sequence, which resulted in noticeable, but not dramatic, deviations of simulated N M R parameters from the experiment. NMR relaxation data were also predicted for RD-hTH and M A P 2 c 1 5 9 - 2 5 4 In the case of RD-hTH, the comparison was possible only in the N-terminal IDR because the protein dimerizes in real samples. As a consequence, relaxation of the well-ordered portion of RD-hTH is incomparable (faster) with that simulated for the monomer. Moreover, broadening of the peaks of the well-ordered region did not allow us to obtain sufficiently sensitive N M R spectra under the conditions used. Comparison of the simulated and experimental relaxation rates in the N-terminal IDR of RD-hTH (Fig. 6 , d-f) and M A P 2 c 1 5 9 - 2 5 4 (Fig. 6, g-i) led to the same general conclusions as for 5RNAR However, the simulation of RD-hTH allowed us to address a particular feature not manifested by 5RNAP, namely the effect of formation of the transient a-helix on the calculated relaxation rates. In the experimental data, the propensity to form an a-helix is reflected by elevated R2 values. Comparison of the relaxation rates calculated from the C36m trajectories (in which the helix was present during most of the simulation time) with those obtained from the A99 and C22* runs (in which the helix quickly unfolded) revealed that the presence of the a-helix in the simulated structures is needed to reproduce the values of R2 (Fig. 6). The 5-ns sliding window was sufficient to match the experimental R2 profile in the C36m simulations, indicating that the dynamics of the IDR of RD-hTH is dominated by relatively short correlation times. Less frequent conformational changes occurring in the a-helical region were not sampled sufficiently and resulted in high standard deviations of R2. Outside of the transient a-helix, the data obtained with the 50-ns window matched the experimental profile well. In the case of M A P 2 c 1 5 9 - 2 5 4 , which has a much lower a-helical propensity, the increase of R2 is hardly visible in the experimental data. Similarly to RD-hTH, relaxation data predicted using the 5-ns window matched the experimental values reasonably well. Quantitative comparison of the force-field reliability To express the discussed differences of the force field performance quantitatively, we applied the metrics developed by Robustelli et al. (45) and calculated normalized forcefield scores based on the RMSDs of the experimentally obtained parameters from the corresponding values predicted from the simulations (Table 1). The calculated combined force-field scores (sc o m bined) show that none of the tested force fields provided superior prediction of all parameters for all proteins. In the case of <5RNAP with a highly charged IDR forming no transient a-helices, relative accuracy of predicting experimental data varied for different parameters. Predicted chemical shifts, reflecting the local backbone conformation, was similar for all force fields (best for C36m). A99 best reproduced NMR data sensitive to long-range intramolecular interactions (RDC, PRE, and relaxation data, described by j N M R ) . The chemical shift RMSD values are within the reported standard deviation of the SPARTA+ predictor (2.45, 1.09, 0.94, and 1.14 ppm for backbone N, C(O), Ca , and C'3 nuclei ( 3 7 ) ) . The quality factor Q of RDC is comparable with typical RMSDs of data predicted from x-ray structures with 2-A resolution (46). PRE and relaxation data deviated from the experimental ranges by less than 10%, with the exception of Tx, which is particularly difficult to predict, as discussed above. C22* provided the best prediction of parameters describing the overall shape of 5RNAP (Rg and SAXS profiles), with the Q-factor of SAXS data equal to 5%. If the individual scores are averaged with the same weights, the overall scores sa l l are better for C22* and A99 than for C36m. In the case of RD-hTH, with the experimentally determined 80% propensity to form an a-helix in its IDR, C36m predicted all experimental data much better than A99 or C22*. The superior performance of C36m, reflected by low sa l l , can be clearly attributed to its ability to maintain the experimentally observed transient a-helix. The quantitative parameters showed that C36m predicted the experimental data well (compared with the reference values discussed above), with the exception of RDC. In the case of M A P 2 c 1 5 9 - 2 5 4 which lacks a well-ordered domain and exhibits only ~20% propensity to form an a-helix, all force fields predicted the experimental data with similar accuracy, reflected by small differences between their scores. It documents that the ability of C36m to maintain transient a-helices loses its significance if the populations of the helical structures are low. The performance of all force fields was good, based on comparison with the reference values discussed in the details of the <5RNAP simulations. In conclusion, the quantitative comparison confirmed the superior performance of C36m in the case when a transient a-helix was present. A99 predicted most accurately parameters influenced by long-range electrostatic interactions (%MR = 1for<5RNAP). Prediction of suitable spin-label positions Calculation of already measured experimental parameters is important for benchmarking of the simulations as illustrated above. However, prediction of so far unknown measurable values is also very useful because it can facilitate experimental design. Selection of residues for placement of 1630 Biophysical Journal 118, 1621-1633, April 7, 2020 Force Field for Disordered Regions paramagnetic labels can serve as an example. Preparation of paramagnetically labeled samples is a time-consuming procedure, including site-directed mutagenesis, expression, purification, and paramagnetic spin labeling of the protein before the NMR PRE measurements. A choice of the label position not providing information about long-range contacts thus represents a considerable waste of time and sources. Reliable M D simulations can be used for in silico prediction how much structural information a label in a certain position can provide and how much such labeling perturbs the native structural ensemble. An example of such prediction for all solvent accessible positions within the RD-hTH is presented in Fig. 7. Based on the simulations using C36m and TIP4P-D, we calculated PRE profiles for all possible positions. To localize solvent accessible positions reporting on longrange contacts with the disordered region, we defined the following score. Integrals of the areas between PRE profiles and a threshold of 0.8 were summed in the disordered region (residues 1-70), except for ± 10 residues in the vicinity of the spin label. In addition, the score was set to zero for residues with solvent accessibility lower than 0.6 according to Fraczkiewicz and Braun (47) because the label should be freely accessible on the surface of protein and should not interfere with any protein conformation. The analysis presented in Fig. 7 a indicates four to five areas well suited for attachment of spin labels for future PRE experiments. Predicted PRE profiles for labels placed in the centers of the suggested areas are presented in Fig. 7, b-e. a 1 III 1 11 i b ^ V \ , ' / v ' b ^ V \ , ' / v ' c Jhf — d \ f \ fe I i i MX0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 Residue number FIGURE 7 Predicted preferential positions of spin labels in RD-hTH (a) and simulated PRE profiles for four selected spin-label positions (b-e). To see this figure in color, go online. CONCLUSIONS The goal of our study was to assess the applicability of currently available M D approaches to hybrid proteins consisting of ordered and disordered regions. We tested the performance of the A99, C22*, and C36m force fields in combination with the TIP3P and TIP4P-D water models. The performance was examined for mostly disordered M A P 2 c 1 5 9 - 2 5 4 containing a low population of prestructured a-helix, for <5RNAP consisting of comparably large well-ordered and disordered domains without transient a-helices, and for RD-hTH consisting of ordered and disordered domains with a highly populated a-helical prestructured motif. Considering the functional importance of transient helical elements, we paid particular attention to the ability to preserve the a-helix during the simulation. The generated structural ensembles were used for predicting a variety of NMR and SAXS parameters and subsequently compared with the experimental data. The TIP4P-D water model performed substantially better than TIP3P or TIPS3P. The differences between force fields were less distinct. The optimal (most universal) combination was CHARMM36m with the TIP4P-D water model, which most efficiently prevented artificial collapse of disordered regions and retained transient a-helical structure elements within the disordered regions. The study also showed that the performance of different force fields and models can vary depending on the actual physical properties of investigated IDRs. Considering the fast pace of force-field development (45,48), the major value of this study is not identification of the best force field available at the time of testing but presentation of a generally applicable benchmarking approach, including N M R relaxation rates. An important feature of the approach is the combination of checked parameters that report on abilities of the force fields to reproduce a broad range of physical properties of the studied molecules. SUPPORTING MATERIAL Supporting Material can be found online at https://doi.Org/10.1016/j.bpj. 2020.02.019. AUTHOR CONTRIBUTIONS J.H. and L.Z. designed the research. V.Z., A.M., PL., A.L., E.N., and V.K. carried out calculations, performed the experiment, and analyzed the data. Z.J., M.M., and K . M . prepared the samples and performed the experiment. V.Z., A.M., PL., J.H., and L.Z. wrote the article, and all authors reviewed the manuscript. ACKNOWLEDGMENTS This research was funded by the Ministry of Education, Youth, and Sport of the Czech Republic (MEYS CR), grant numbers LTC17078 Biophysical Journal 118, 1621-1633, April 7, 2020 1631 Zapletal et al. (Inter­Excellence I nter­Cost), LTAUSA18168 (I nter­Excellence I nter­Action), and LQ1601 (National Sustainability Programme I I Project CEI TEC 2020). Computational resources were provided by CESNET (LM2015042) and the CERIT Scientific Cloud (LM2015085) under the program "Projects of Large Research, Development, and Innovations I nfrastructures" funded by MEYS C R and by IT4Innovations National Supercomputing Center (LM2015070) under the program "Large I nfrastructures for Research, Experimental Development and Innovations" funded by MEYS CR. The Czech I nfrastructure for Integrative Structural Biology (CI I SB) research infrastructure project LM2018127 funded by MEYS C R is gratefully acknowledged for the partial financial support of the measurements at the Josef Dadok National N M R Centre and at X­Ray Diffraction and BioSAXS Core Facilities, CEI TEC­Masaryk University. REFERENCES 1. Uversky, V. N. 2002. Natively unfolded proteins: a point where biology waits for physics. Protein Sci. 11:739­756. 2. Dunker, A. K., C. J. Brown, ..., Z. Obradovič. 2002. Intrinsic disorder and protein function. Biochemistry. 41:6573­6582. 3. Tompa, P. 2011. Unstructural biology coming of age. Curr. Opin. Struct. Biol. 21:419^25. 4. Chi, S.­W., D.­H. Kim, ..., K. H. Han. 2007. Pre­structured motifs in the natively unstructured preS 1 surface antigen of hepatitis B virus. Protein Sci. 16:2108­2117. 5. Fuxreiter, M., I . Simon, P. Tompa. 2004. Preformed structural elements feature in partner recognition by intrinsically unstructured proteins. J. Mol. Biol. 338:1015­1026. 6. Vacic, V., C J. Oldfield, ..., A. K. Dunker. 2007. Characterization of molecular recognition features, MoRFs, and their binding partners. J. Proteome Res. 6:2351­2366. 7. Nováček, J., L. Židek, and V. Sklenář. 2014. Toward optimal­resolution NMR of intrinsically disordered proteins. J. Magn. Reson. 241:41­52. 8. Nowakowski, M., S. Saxena, ..., W. Kožmiňski. 2015. Applications of high dimensionality experiments to biomolecular NMR. Prog. Nucl. Magn. Reson. Spectrosc. 90­91:49­73. 9. Papoian, G. A. 2008. Proteins with weakly funneled energy landscapes challenge the classical structure­function paradigm. Proc. Natl. Acad. Sci. USA. 105:14237­14238. 10. Rabatinová, A., H. Šanderová, L. Krásný. 2013. The 5 subunit of RNA polymerase is required for rapid changes in gene expression and competitive fitness of the cell. J. Bacteriol. 195:2603­2611. 11. Papoušková, V., P. Kadeřávek, L. Žídek. 2013. Structural study of the partially disordered full­length 5 subunit of RNA polymerase from Bacillus subtilis. ChemBioChem. 14:1772­1779. 12. Nagatsu, T., M . Levitt, and S. Udenfriend. 1964. Tyrosine hydroxylase. The initial step in norepinephrine biosynthesis. J. Biol. Chem. 239:2910­2917. 13. Molinoff, P. B., and J. Axelrod. 1971. Biochemistry of catecholamines. Annu. Rev. Biochem. 40:465­500. 14. Louša, P., H. Nedozrálová, J. Hritz. 2017. Phosphorylation of the regulatory domain of human tyrosine hydroxylase 1 monitored using non­uniformly sampled NMR. Biophys. Chem. 223:25­29. 15. Jansen, S., K. Melková, L. Žídek. 2017. Quantitative mapping of microtubule­associated protein 2c (MAP2c) phosphorylation and regulatory protein 14­3­3^­binding sites reveals key differences between MAP2c and its homolog Tau. J. Biol. Chem. 292:6715­6727. 16. Melková, K., V. Zapletal, L . Žídek. 2018. Functionally specific binding regions of microtubule­associated protein 2c exhibit distinct conformations and dynamics. J. Biol. Chem. 293:13297­13309. 17. Melková, K , V. Zapletal,..., L. Žídek. 2019. Structure and functions of microtubule associated proteins Tau and MAP2c: similarities and differences. Biomolecules. 9:E105. 18. Motáčková, V , J. Nováček, ..., V. Sklenář. 2010. Strategy for complete NMR assignment of disordered proteins with highly repetitive sequences based on resolution­enhanced 5D experiments. J. Biomol. NMR. 48:169­177. 19. Nováček, J., L. Janda,..., V. Sklenář. 2013. Efficient protocol for backbone and side­chain assignments of large, intrinsically disordered proteins: transient secondary structure analysis of 49.2 kDa microtubule associated protein 2c. J. Biomol. NMR. 56:291­301. 20. Ottiger, M . , F. Delaglio, and A. Bax. 1998. Measurement of J and dipolar couplings from simplified two­dimensional N M R spectra. J. Magn. Reson. 131:373­378. 21. Srb, P., J. Nováček, L. Žídek. 2017. Triple resonance 1 5 N N M R relaxation experiments for studies of intrinsically disordered proteins. J. Biomol. NMR. 69:133­146. 22. Korzhnev, D. M., M . Billeter, ..., V. Y. Orekhov. 2001. N M R studies of Brownian tumbling and internal motions in proteins. Prog. Nucl. Magn. Reson. Spectrosc. 38:197—266. 23. Efron, B. 1979. Bootstrap methods: another look at the jackknife. Ann. Stat. 7:1­26. 24. Petoukhov, M . V., D. Franke, D. I. Svergun. 2012. New developments in the ATSAS program package for small­angle scattering data analysis. J. Appl. Cryst. 45:342­350. 25. Riback, J. A., M . A. Bowman, ..., T. R. Sosnick. 2017. Innovative scattering analysis shows that hydrophobic disordered proteins are expanded in water. Science. 358:238­241. 26. Lindorff­Larsen, K , S. Piana, D. E. Shaw. 2010. I mproved sidechain torsion potentials for the Amber ff99SB protein force field. Proteins. 78:1950­1958. 27. Piana, S., K. Lindorff­Larsen, and D. E. Shaw. 2011. How robust are protein folding simulations with respect to force field parameterization? Biophys. J. 100:L47­L49. 28. Huang, J., S. Rauscher,..., A. D. MacKerell, Jr. 2017. CHARMM36m: an improved force field for folded and intrinsically disordered proteins. Nat. Methods. 14:71­73. 29. Jorgensen, W. L. 1981. Quantum and statistical mechanical studies of liquids. 10. Transferable intermolecular potential functions for water, alcohols, and ethers. Application to liquid water. J. Am. Chem. Soc. 103:335­340. 30. MacKerell, A. D., D. Bashford, ..., M . Karplus. 1998. All­atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B. 102:3586­3616. 31. Piana, S., A. G. Donchev, D. E. Shaw. 2015. Water dispersion interactions strongly influence simulated structural properties of disordered protein states. J. Phys. Chem. B. 119:5113­5123. 32. Hess, B., H. Bekker, J. G. E. M . Fraaije. 1997. LI NCS: a linear constraint solver for molecular simulations. J. Comput. Chem. 18:1463­1472. 33. Essmann, U., L. Perera, L. G. Pedersen. 1995. A smooth particle mesh Ewald method. J. Chem. Phys. 103:8577­8593. 34. Berendsen.H.J.CJ.P.M.Postma, ...J.R.Haak. 1984. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81:3684—3690. 35. Bussi, G., D. Donadio, and M . Parrinello. 2007. Canonical sampling through velocity rescaling. J. Chem. Phys. 126:014101. 36. Parrinello, M., and A. Rahman. 1981. Polymorphic transitions in single crystals: anew molecular dynamics method. J. Appl. Phys. 52:7182­7190. 37. Shen, Y , and A. Bax. 2010. SPARTA+: a modest improvement in empirical NMR chemical shift prediction by means of an artificial neural network. J. Biomol. NMR. 48:13­22. 38. Nielsen, J. T., and F. A. A. Mulder. 2018. POTENCI: prediction of temperature, neighbor and pH­corrected chemical shifts for intrinsically disordered proteins. J. Biomol. NMR. 70:141­165. 39. Zweckstetter, M . 2008. NMR: prediction of molecular alignment from structure using the PALES software. Nat. Protoc. 3:679­690. 40. Nodet, G., L. Salmon, M . Blackledge. 2009. Quantitative description of backbone conformational sampling of unfolded proteins at 1632 Biophysical Journal 118, 1621­1633, April 7, 2020 Force Field for Disordered Regions amino acid resolution from N M R residual dipolar couplings. J. Am. Chem. Soc. 131:17908-17918. 41. Salmon, L., G. Nodet, ..., M . Blackledge. 2010. NMR characterization of long-range order in intrinsically disordered proteins. J. Am. Chem. Soc. 132:8407-8418. 42. Sezer, D., J. H. Freed, and B. Roux. 2008. Simulating electron spin resonance spectra of nitroxide spin labels from molecular dynamics and stochastic trajectories. J. Chem. Phys. 128:165106. 43. Salvi, N., A. Abyzov, and M . Blackledge. 2016. Multi-timescale dynamics in intrinsically disordered proteins from N M R relaxation and molecular simulation. J. Phys. Chem. Lett. 7:2483-2489. 44. Urbaiiczyk, M., D. Bernin,..., K. Kazimierczuk. 2013. Iterative thresholding algorithm for multiexponential decay applied to PGSE NMR data. Anal. Chem. 85:1828-1833. 45. Robustelli, P., S. Piana, and D. E. Shaw. 2018. Developing a molecular dynamics force field for both folded and disordered protein states. Proc. Natl. Acad. Sci. USA. 115:E4758-E4766. 46. Bax, A. 2003. Weak alignment offers new N M R opportunities to study protein structure and dynamics. Protein Sci. 12:1-16. 47. Fraczkiewicz, R., and W. Braun. 1998. Exact and efficient analytical calculation of the accessible surface areas and their gradients for macromolecules. J. Comput. Chem. 19:319-333. 48. Song, D., R. Luo, and H.-F. Chen. 2017. The IDP-specific force field ffl4IDPSFF improves the conformer sampling of intrinsically disordered proteins. J. Chem. Inf. Model. 57:1166-1178. Biophysical Journal 118, 1621-1633, April 7, 2020 1633