pH dependence of charge multipole moments in proteins

TL;DR

This study investigates how pH variations affect the charge multipole moments of proteins, revealing changes in both magnitude and orientation that impact protein interactions and assembly processes.

Contribution

It introduces a multipole expansion approach to analyze pH-dependent charge distributions on proteins, highlighting orientation changes of multipoles with pH.

Findings

pH influences the magnitude of protein charge multipoles

pH affects the spatial orientation of multipole axes

Implications for protein interaction modeling

Abstract

Electrostatic interactions play a fundamental role in the structure and function of proteins. Due to ionizable amino acid residues present on the solvent-exposed surfaces of proteins, the protein charge is not constant but varies with the changes in the environment -- most notably, the pH of the surrounding solution. We study the effects of pH on the charge of four globular proteins by expanding their surface charge distributions in terms of multipoles. The detailed representation of the charges on the proteins is in this way replaced by the magnitudes and orientations of the multipole moments of varying order. Focusing on the three lowest-order multipoles -- the total charge, dipole, and quadrupole moment -- we show that the value of pH influences not only their magnitudes, but more notably and importantly also the spatial orientation of their principal axes. Our findings imply important consequences for the study of protein-protein interactions and the assembly of both proteinaceous shells and patchy colloids with dissociable charge groups.

Figures

Tables

Table 1. Table 1: Summary of results for the p ​ H 𝑝 𝐻 pH dependence of the magnitudes of multipole components when cysteine acidity is not considered. Listed are the isoelectric point, p ​ I 𝑝 𝐼 pI , where the monopole moment (total charge) vanishes, and the magnitudes of the dipole and quadrupole moment in this point (subscript I 𝐼 I ). We also list the magnitudes of the monopole, dipole, and quadrupole moments at neutral p ​ H = 7 𝑝 𝐻 7 pH=7 (subscript n 𝑛 n ). Lastly, we list the quadrupole ratio ⟨ | Q x ​ x / Q z ​ z | ⟩ p ​ H subscript delimited-⟨⟩ subscript 𝑄 𝑥 𝑥 subscript 𝑄 𝑧 𝑧 𝑝 𝐻 \langle|Q_{xx}/Q_{zz}|\rangle_{pH} , averaged across the entire range of p ​ H 𝑝 𝐻 pH values. The magnitudes of the multipole moments correspond to total charge in the case of the monopole, Cartesian norm in the case of the dipole, and Q = [ ( Q x ​ x − Q y ​ y ) 2 + Q z ​ z 2 ] 1 / 2 𝑄 superscript delimited-[] superscript subscript 𝑄 𝑥 𝑥 subscript 𝑄 𝑦 𝑦 2 superscript subscript 𝑄 𝑧 𝑧 2 1 2 Q=[(Q_{xx}-Q_{yy})^{2}+Q_{zz}^{2}]^{1/2} in the case of the quadrupole.

PDB	CYS	$pI$	${\mu }_I$ $\left[e_0\times \text{nm}\right]$	$Q_I$ $\left[e_0\times {\text{nm}}^2\right]$	$q_n$ $\left[e_0\right]$	${\mu }_n$ $\left[e_0\times \text{nm}\right]$	$Q_n$ $\left[e_0\times {\text{nm}}^2\right]$	${⟨\left|Q_{xx}/Q_{zz}\right|⟩}_{pH}$
2lyz	no	$11.08$	$1.85$	$3.86$	$8.10$	$2.67$	$2.92$	$1.09$
1e7h	no	$4.78$	$2.27$	$34.68$	$-23.45$	$7.86$	$38.01$	$0.93$
2blg	no	$4.47$	$5.77$	$2.01$	$-7.78$	$7.43$	$2.64$	$0.91$
2ms2	no	$9.78$	$4.88$	$2.34$	$2.00$	$4.91$	$2.78$	$1.02$

Table 2. Table S1: Number of charges present on the solvent-exposed surfaces of proteins used in our study as a function of the RSA cut-off c 𝑐 c , with or without considering the cysteine (CYS) acidity. The number of charges present on the surface decreases with increasing cut-off. Outlined in bold is the cut-off of c = 0.25 𝑐 0.25 c=0.25 , used in the majority of the main text.

	0.1	0.2	0.25	0.3	0.4	0.5
2lyz (no CYS)	30	30	30	28	26	23
2lyz (with CYS)	34	33	33	31	27	24
1e7h (no CYS)	164	162	157	152	145	119
1e7h (with CYS)	192	188	182	173	162	127
2ms2 (no CYS)	23	23	22	22	19	17
2ms2 (with CYS)	25	25	24	24	21	19
2blg (no CYS)	50	49	46	45	42	36
2blg (with CYS)	51	50	47	46	43	37

Table 3. Table S2: Intrinsic p ​ K a 𝑝 subscript 𝐾 𝑎 pK_{a} values of amino acid functional groups in bulk dilute aqueous solutions. Values taken from Ref.  (Haynes, 2016 ) .

	ASP	GLU	TYR	ARG	HIS	LYS	CYS
$p{K}_a$	3.71	4.15	10.10	12.10	6.04	10.67	8.14

Table 4. Table S3: Circumscribed radii of the proteins used in the study: lysozyme (2lys), human serum albumin (1e7h), β 𝛽 \beta -lactoglobulin (2blg), and subunit A of the phage MS2 capsid protein (2ms2). The radii were obtained as the largest distance of any C α subscript 𝐶 𝛼 C_{\alpha} atom from the center-of-mass of a protein.

	2lyz	1e7h	2blg	2ms2
$R$ [nm]	0.987	1.97	1.06	1.12

Table 5. Table S4: Summary of results for the p ​ H 𝑝 𝐻 pH dependence of the magnitudes of multipole components, comparing our simple model, where every amino acid group is assigned the same p ​ K a 𝑝 subscript 𝐾 𝑎 pK_{a} value [cf. Table  S2 ], with two models introducing local effects to the p ​ K a 𝑝 subscript 𝐾 𝑎 pK_{a} values of individual amino acid residues. From the latter two models, the first one uses empirical predictions of PROPKA software, which includes site-site interactions  Søndergaard et al. ( 2011 ) . The second model takes into account the local p ​ K a 𝑝 subscript 𝐾 𝑎 pK_{a} shifts imparted by the electrostatic potential of the protein in the presence of 1:1 salt with bulk concentration c 0 = 1 subscript 𝑐 0 1 c_{0}=1 M, derived using linearized DH equation with charge regulation boundary condition. Listed are the isoelectric point, p ​ I 𝑝 𝐼 pI , where the monopole moment (total charge) vanishes, and the magnitudes of the dipole and quadrupole moment in this point (subscript I 𝐼 I ). We also list the magnitudes of the monopole, dipole, and quadrupole moments at neutral p ​ H = 7 𝑝 𝐻 7 pH=7 (subscript n 𝑛 n ). Lastly, we list the quadrupole ratio ⟨ | Q x ​ x / Q z ​ z | ⟩ p ​ H subscript delimited-⟨⟩ subscript 𝑄 𝑥 𝑥 subscript 𝑄 𝑧 𝑧 𝑝 𝐻 \langle|Q_{xx}/Q_{zz}|\rangle_{pH} , averaged across the entire range of p ​ H 𝑝 𝐻 pH values. The magnitudes of the multipole moments correspond to total charge in the case of the monopole, Cartesian norm in the case of the dipole, and Q = [ ( Q x ​ x − Q y ​ y ) 2 + Q z ​ z 2 ] 1 / 2 𝑄 superscript delimited-[] superscript subscript 𝑄 𝑥 𝑥 subscript 𝑄 𝑦 𝑦 2 superscript subscript 𝑄 𝑧 𝑧 2 1 2 Q=[(Q_{xx}-Q_{yy})^{2}+Q_{zz}^{2}]^{1/2} in the case of the quadrupole.

PDB	model	$pI$	${\mu }_I$ $\left[e_0\times \text{nm}\right]$	$Q_I$ $\left[e_0\times {\text{nm}}^2\right]$	$q_n$ $\left[e_0\right]$	${\mu }_n$ $\left[e_0\times \text{nm}\right]$	$Q_n$ $\left[e_0\times {\text{nm}}^2\right]$	${⟨\left|Q_{xx}/Q_{zz}\right|⟩}_{pH}$
2lyz	no CYS	$11.08$	$1.85$	$3.86$	$8.10$	$2.67$	$2.92$	$1.09$
2lyz	CYS	$10.53$	$1.80$	$3.18$	$7.90$	$2.61$	$2.82$	$1.24$
2lyz	PROPKA	$11.38$	$2.60$	$3.85$	$8.70$	$4.43$	$3.75$	$1.16$
2lyz	$c_0=1$ M	$11.08$	$2.60$	$5.16$	$8.38$	$2.75$	$3.08$	$1.08$
1e7h	no CYS	$4.78$	$2.27$	$34.68$	$-23.45$	$7.86$	$38.01$	$0.93$
1e7h	CYS	$4.78$	$2.27$	$34.69$	$-25.14$	$8.08$	$39.17$	$0.95$
1e7h	PROPKA	$5.12$	$4.40$	$31.23$	$-19.27$	$8.61$	$33.14$	$0.88$
1e7h	$c_0=1$ M	$4.78$	$3.91$	$53.41$	$-24.41$	$8.61$	$39.47$	$0.93$
2blg	no CYS	$4.47$	$5.77$	$2.01$	$-7.78$	$7.43$	$2.64$	$0.91$
2blg	CYS	$4.47$	$5.77$	$2.01$	$-7.85$	$7.49$	$2.58$	$0.97$
2blg	PROPKA	$4.78$	$3.95$	$4.28$	$-8.62$	$5.41$	$5.77$	$1.19$
2blg	$c_0=1$ M	$4.38$	$10.46$	$4.26$	$-8.00$	$7.55$	$2.70$	$0.92$
2ms2	no CYS	$9.78$	$4.88$	$2.34$	$2.00$	$4.91$	$2.78$	$1.02$
2ms2	CYS	$8.88$	$4.65$	$2.64$	$1.87$	$4.87$	$2.68$	$1.07$
2ms2	PROPKA	$9.03$	$4.22$	$1.49$	$1.00$	$4.35$	$2.19$	$1.24$
2ms2	$c_0=1$ M	$9.68$	$5.55$	$2.47$	$2.01$	$4.92$	$2.78$	$1.03$

Equations

$$\mathrm{RSA}({\mathrm{AA}}_k)⩾0.25\Rightarrow {\mathrm{AA}}_k\in \mathrm{exposed}.$$

$$q_k^{\pm }=\frac{\pm 1}{1+e^{\pm \ln 10(pH-p{K}_a^{(k)})}}$$

$$q_k^{\pm }=\frac{\pm 1}{1+e^{\pm \ln 10(pH-p{K}_a^{(k)})\mp \beta e_0\psi ({\mathbf{r}}_k)}},$$

$$\sigma (\Omega )$$

$${\sigma }_{lm}=\sum _{k\in \mathrm{AA}}{q}_k{Y}_{lm}^{\ast }({\Omega }_k).$$

$$q=\sum _{k\in \mathrm{AA}}{q}_k,$$

$$\mathbf{\mu }$$

$$\mathcal{Q}$$

$$q_k^{\pm }=\frac{\pm e_0}{1+e^{\pm \ln 10(pH-p{K}_a^{(k)})\mp \beta e_0\psi ({\mathbf{r}}_k)}},$$

$$p{K}_a^{(k)}\to p{K}_a^{(k)}+\frac{\beta e_0\psi ({\mathbf{r}}_k)}{\ln 10}.$$

$${\nabla }^2{\psi }_1(r,\Omega )=0;r\le R\mathrm{and}{\nabla }^2{\psi }_2(r,\Omega )={\kappa }^2{\psi }_2(r,\Omega );r\ge R,$$

$${\psi }_1(R,\Omega )={\psi }_2(R,\Omega )\forall \Omega ,$$

$$\varepsilon_{p}\varepsilon_{0}\frac{\partial\psi_{1}(r,\Omega)}{\partial r}\Big{|}_{r=R}-\varepsilon_{w}\varepsilon_{0}\frac{\partial\psi_{2}(r,\Omega)}{\partial r}\Big{|}_{r=R}=\sigma(\Omega)=\sigma^{(0)}(\Omega)+\sigma^{(1)}(\Omega)\beta e_{0}\psi(R,\Omega)\quad\forall\Omega.$$

$$q_k^{\pm }\approx \frac{\pm 1}{1+e^{\pm \ln 10\left(pH-p{K}_a^{(k)}\right)}}\mp \frac{e^{\pm \ln 10\left(pH-p{K}_a^{(k)}\right)}}{{\left[1+e^{\pm \ln 10\left(pH-p{K}_a^{(k)}\right)}\right]}^2}\times \beta e_0\psi (R,{\Omega }_k)=q_k^{\pm ,(0)}+q_k^{\pm ,(1)}\beta e_0\psi (R,{\Omega }_k).$$

$$\sigma (\Omega )={\sigma }^{(0)}(\Omega )+{\sigma }^{(1)}(\Omega )\beta e_0\psi (R,\Omega )=\frac{e_0}{R^2}\left[\sum _{k\in \mathrm{AA}}{q}_k^{(0)}\delta (\Omega -{\Omega }_k)+\sum _{k\in \mathrm{AA}}{q}_k^{(1)}\delta (\Omega -{\Omega }_k)\beta e_0\psi (R,\Omega )\right],$$

$${\sigma }_{lm}^{(0)}=\sum _{k\in \mathrm{AA}}{q}_k^{(0)}{Y}_{lm}^{\ast }({\Omega }_k)\mathrm{and}{\sigma }_{lm}^{(1)}=\sum _{k\in \mathrm{AA}}{q}_k^{(1)}{Y}_{lm}^{\ast }({\Omega }_k).$$

$${\psi }_1=\sum _{l,m}{A}_{lm}{\left(\frac{r}{R}\right)}^l{Y}_{lm}(\Omega )\mathrm{and}{\psi }_2=\sum _{l,m}{B}_{lm}{k}_l(\kappa r)Y_{lm}(\Omega ),$$

$$f_u(\kappa R)B_{uv}={\sigma }_{uv}^{(0)}+\beta e_0\sum _{l,m}\sum _{l^{\prime },m^{\prime }}{\sigma }_{lm}^{(1)}{B}_{l^{\prime }{m}^{\prime }}{k}_{l^{\prime }}(\kappa R)T(uv|lm|l^{\prime }{m}^{\prime }).$$

$$f_l(\kappa R)=\frac{R^2\kappa {\epsilon }_w{\epsilon }_0}{e_0}\left[\frac{{\epsilon }_p}{{\epsilon }_w}\frac{l{k}_l(\kappa R)}{\kappa R}-k_l^{\prime }(\kappa R)\right],$$

$$T(uv|lm|l^{\prime }{m}^{\prime })=\sqrt{\frac{(2u+1)(2l+1)(2l^{\prime }+1)}{4\pi }}\times \left(\begin{array}{ccc}u& l& l^{\prime }\\ 0& 0& 0\end{array}\right)\left(\begin{array}{ccc}u& l& l^{\prime }\\ v& m& m^{\prime }\end{array}\right).$$

$${\overline{pK}}_a=\frac{1}{N}\sum _{k\in \mathrm{AA}}\frac{|\beta e_0\psi ({\mathbf{r}}_k)|}{\ln 10}.$$

Full text

INTRODUCTION

Electrostatic interactions are an important part of the long- and short-range interactions in the biological environment. Their understanding is usually based on the framework of the generalized DLVO (Derjaguin-Landau-Verwey-Overbeek) theory of interactions between colloids in ionic solutions, where the canonical electrostatic and van der Waals components are supplemented by the solvent structure effects of either hydration or hydrophobic type (Leckband and Israelachvili, 2001). This general decomposition of the interactions remains valid also for proteins in aqueous solution, where it is in addition augmented by short-range recognition and docking interactions as well as specific ion effects (Leckband and Sivasankar, 1999; Piazza, 2004).

While the van der Waals interactions are a functional of the dielectric response function (Simonson, 2003; Woods et al., 2016), the electrostatic interactions within and between the proteins are based on their charge distributions. The electrostatic interactions dominate many aspects of protein behavior, and can be modelled on different levels of detail (Perutz, 1978; Warshel et al., 2006; Gitlin et al., 2006). One of the defining differences between the specific electrostatics of proteins and the generic electrostatics of colloids that needs to be taken into account is the existence of ionizable amino acid residues in proteins. The interactions and energetics of these residues enable local charge separation, implying a protein-specific distribution of charges (Nozaki and Tanford, 1967; Warshel et al., 2006) and influencing the protein-protein electrostatic interactions (Lund and Jönsson, 2005; Barosso da Silva et al., 2006; Barroso da Silva and Jönsson, 2009). The mechanism of charge separation in turn creates a distinction between the undissociated chargeable groups buried inside the proteins on the one hand and the solvent-exposed and dissociable surface charges on the other (Finkelstein and Ptitsyn, 2016). This distinction can be blurred, as internal ionizable groups can be to some extent dielectrically shielded even in the strongly hydrophobic protein core (Gong et al., 2008; Isom et al., 2010).

Identification of the dissociable, solvent-exposed amino acid residues is the first step in obtaining the description of charge distributions in proteins. Afterwards, one needs to take into account the proper description of the dissociation mechanism for the deprotonated carboxylate of aspartic and glutamic acid, deprotonated hydroxyl of the tyrosine phenyl group, the protonated amine group of arginine and lysine, and the protonated secondary amine of histidine (Gitlin et al., 2006; Jensen, 2008). Since the dissociation process and the local electrostatic field are coupled via the charge regulation mechanism (Markovich et al., 2016), the dissociation rate depends on the local $pH$ that can be obtained only self-consistently (Longo et al., 2012; Nap et al., 2014). In principle, only detailed ab initio simulations can provide a detailed quantification of the partial charges buried inside or exposed on a protein surface. These simulations are, however, usually hampered by the sheer size and number of atoms one needs to invoke in order to achieve a necessary amount of realism for the calculations (Lee et al., 2013; Adhikari et al., 2014; Eifler et al., 2016).

A proper quantification of the electrostatic interactions in proteins requires the encoding of not only the magnitudes of the charges but also of their distribution in space (Gramada and Bourne, 2006; Hoppe, 2013). The latter can be represented to any desirable accuracy by a multipole expansion of the charge density, where each term in the multipole series describes a deviation with a specific symmetry from the zeroth-order, spherically symmetric distribution (Paulini et al., 2005). There are many variants of the multipole expansion, formulated either on the level of amino acids or on the level of complete proteins (Yuan et al., 2014a, b; Fletcher and Popelier, 2016), with the most straightforward being the irreducible spherical representation of the multipoles (Gramada and Bourne, 2006). This representation is obtained by mapping the charge distribution on the original solvent-accessible protein surface onto a sphere circumscribed to the protein (Hoppe, 2013; Postarnakevich and Singh, 2009; Arzenšek et al., 2015).

The multipole expansion provides a bridge between a coarse-grained description of the charge density and its detailed microscopic description, the level of detail depending on the multipole order of the expansion used. As such, even multipoles of lower order can provide a signature of charge distributions in molecules (Gramada and Bourne, 2006; Platt and Silverman, 1996; Nakamura and Wada, 1985; Kim et al., 2006), and a small number of higher-order multipoles can account for almost all of the electrostatic field in the aqueous solvent (Hoppe, 2013). In the presence of ionic screening this otherwise standard result is modified, as the effects of charge anisotropy and higher-order multipole moments then extend to the far-field region. In contrast to the standard multipole expansion, the screened electrostatic potential retains the full directional dependence of all multipole moments, an important difference which is often overlooked (Rowan et al., 2000; Kjellander and Ramirez, 2008; Lošdorfer Božič and Podgornik, 2013; Kjellander, 2016).

With a few exceptions (Hoppe, 2013), the $pH$ dependence of charge distributions in proteins has been standardly studied mostly on the level of the spherically symmetric total charge, that is, the zeroth multipole moment. This approach, however, completely neglects the orientational dependence of the interactions. Due to the peculiarity of the multipole expansion for screened electrostatic interactions, the orientational effects are present even at the lowest multipolar order and can thus – contrary to the existing analytical models – never really be ignored (Rowan et al., 2000; Kjellander and Ramirez, 2008; Lošdorfer Božič and Podgornik, 2013; Kjellander, 2016). In addition, for higher-order multipoles the two components of the multipole expansion, the magnitude of the multipole moments and their spatial distribution, can hardly be separated. We thus investigate not only the $pH$ dependence of the magnitudes but also of the directions of the principal axes of the charge dipole and quadrupole moments in proteins, a phenomenon missed also in cases where the higher multipole terms were included (Li et al., 2015; McManus et al., 2016; Grant, 2001; Matthew et al., 1985). It is in fact the orientational dependencies that are of particular importance for the local electrostatic interactions in numerous contexts, even when the coupling between $pH$, protonation states, and protein conformation is not taken into account. The results of our model should thus be relevant for the studies of interaction and assembly in ordered protein assemblies (such as a proteinaceous virus shell (Šiber et al., 2012)), between charged Janus colloids (Hieronimus et al., 2016), in the general context of patchy globular proteins, colloids, and polyelectrolytes (Yigit et al., 2015a, b, 2017; Bianchi et al., 2011, 2014; Stipsitz et al., 2015), or possibly as a driving mechanism for local packing symmetry transitions in the proteinaceous capsid engineering in the presence of supercharged proteins (Sasaki and Hilvert, 2016).

MATERIALS AND METHODS

Protein dataset

For our study we chose four protein structures from the RCSB Protein Data Bank (PDB) (Berman et al., 2000): hen egg-white lysozyme (2lyz), human serum albumin (1e7h), bovine $\beta$-lactoglobulin (2blg), and phage MS2 capsid protein (2ms2, subunit A). These proteins are biologically well-studied, with known structure and role. As such, they are often used to examine the role of electrostatic contributions in protein systems and to study protein assembly and aggregation, where higher-order charge multipoles play a role (Hoppe, 2013; Nap et al., 2014; Lund and Jönsson, 2005; Yigit et al., 2017; Cardinaux et al., 2007). $\beta$-lactoglobulin, for instance, is known to form dimers or tetramers depending on the $pH$. The capsid protein of MS2 also first assembles into trimers and from there into full capsids made of 180 proteins, imparting to the capsids a $pH$-dependent stability. As our interest lies in exploring the $pH$ dependence of multipole moments and their orientations that could relate to higher-order assembly, we study here only the monomer forms of each protein. In addition, the proteins in our dataset have globular geometry, easily approximated by a sphere, and most of them are fairly small, with lysozyme and MS2 capsid protein being composed of 129 amino acids (AA), and $\beta$-lactoglobulin of 162 AA. Human serum albumin (HSA) is the largest of the four, consisting of 585 AA. While we will focus on the four proteins in the dataset to study the general properties of our model, we note that our approach can also be straightforwardly implemented for any globular protein.

To obtain the surface charges on each protein at a given $pH$ value, we first determine which AA residues are solvent-accessible, and compute the charge on them by using the canonical static dissociation constant $pK_{a}$ value pertaining to each AA type. By projecting the positions of the charged residues onto a sphere, we obtain the surface charge distribution of each protein, which we then use to compute the electrostatic multipole moments. Details of each of these steps are laid out in the following Subsections, and a sketch of the model is shown in Fig. 1.

Relative solvent accessibility of amino acid residues

Relative solvent accessibility (RSA) of an amino acid indicates its degree of burial in the three-dimensional protein structure, and is important in determining which AA residues can contribute to the surface charge of the protein. We obtain the RSA of each AA in a protein with the help of the UCSF Chimera software (Pettersen et al., 2004; Sanner et al., 1996) by normalizing the solvent-exposed surface area of each residue in the protein structure with the surface area of the same type of residue in a reference state (Bendell et al., 2014). The classification of AA residues as “buried” or “exposed” is then done on the basis of an RSA cut-off $c$, which typically ranges between 5% and 30% (Pollastri et al., 2002; Yuan et al., 2002; Chen and Zhou, 2005; Wu et al., 2017). We opt for a cut-off of $c=0.25$, thus defining as exposed all amino acid residues $\mathrm{AA}_{k}$ with an RSA value greater than 25%:

[TABLE]

The choice of this (arbitrary) cut-off influences the number of charged residues that will be taken into account in our calculation of the surface charge of the protein. How the number of charges varies with the selection of the RSA cut-off is shown in Table S1 in the Supporting Material. Similar variations in the number of surface charges could also occur, for instance, due to a different choice of the normalization values (reference states) in the calculation of the RSA (Tien et al., 2013).

For consistency, we use the RSA cut-off of $c=0.25$ throughout the paper, unless specified otherwise. However, we also test our predictions for several other values of the cut-off and show that it bears no influence on the conclusions obtained in our work.

pH dependence of charged amino acid residues

Once we know which AAs are solvent-exposed and thus dissociable, we can determine their charge at a given $pH$ value. The charged residues we consider are the aspartic acid (ASP), glutamic acid (GLU), tyrosine (TYR), arginine (ARG), lysine (LYS), and histidine (HIS), while the case of cysteine (CYS) protonation is considered separately in the Discussion and Supporting Material. To obtain the charge on each of the residues as a function of $pH$ we use the acid-base dissociation constants $pK_{a}$ from Ref. (Haynes, 2016) (and listed in Table S2 in the Supporting Material). The degree of dissociation of each amino acid AAk as a function of its $pK_{a}^{(k)}$ and $pH$ is then given by the Henderson-Hasselbach equation:

[TABLE]

for bases ($q^{+}_{k}>0$) and acids ($q^{-}_{k}<0$), respectively. We note that eq. (2) can be expanded to include the local electrostatic potential at the position of each charge, $\psi(\mathbf{r}_{k})$, inducing a $pK_{a}$ shift (Ninham and Parsegian, 1971; Krishnan, 2017):

[TABLE]

where $\beta=1/k_{B}T$, $T$ is the room temperature, $k_{B}$ the Boltzmann constant, and $e_{0}$ is the elementary charge. The $pK_{a}$ shift due to local electrostatic potential is hard to neglect, especially at low salt concentrations. However, reformulating the Tanford-Kirkwood model (Tanford and Kirkwood, 1957; Barroso da Silva et al., 2001) by using a linearized Debye-Hückel (DH) theory on a low dielectric constant sphere with charge regulation boundary condition, we show that for large (physiological) salt concentrations the effect of the $pK_{a}$ shift due to the electrostatic potential is relatively small and eq. (2) can be used. Details of the DH model are given in the Supporting Material [eqs. (S1)-(S12)], and its consequences commented on in the Discussion.

Electrostatic effects are not the only possible factor inducing a local $pK_{a}$ shift of a given amino acid in a protein, as hydrogen bonding and desolvation effects often play an important role. Numerous methods exist for the prediction of $pK_{a}$ values on different levels of detail, ranging from various molecular dynamics models and ab initio quantum mechanical approaches to empirical models, which trade the detailed description of a system for a very fast computational time (Alexov et al., 2011). In order to further compare our results, obtained by assigning the same $pK_{a}$ value to each amino acid residue group, with a more realistic model including site-site interactions, we use the PROPKA3.1 software (Søndergaard et al., 2011) to predict local $pK_{a}$ values of each amino acid residue. PROPKA is a widely-used empirical software which uses three-dimensional structure of proteins to estimate desolvation effects and intraprotein interactions based on the positions and chemical nature of the groups proximate to the $pK_{a}$ sites. The $pK_{a}$ values predicted by PROPKA do not include salt effects as an input parameter, and thus provide another layer of contrast to the $pK_{a}$ values shifted due to electrostatic effects, described above. The results obtained in this way are presented in the Discussion and the Supporting Material, where we show that for the purposes of our study, the simple model that we use fares well even when compared to a more realistic model.

Surface charge distribution and electrostatic multipoles

With the approach outlined in the previous Subsections we can obtain, at any value of $pH$, the positions of charged residues for a given protein, $\bm{r}_{k}=(x_{k},y_{k},z_{k})$, and the (fractional) charge they carry, $q_{k}$. To obtain a surface charge distribution we then project them onto a spherical surface, $\bm{r}_{k}=(R,\Omega_{k})=(R,\vartheta_{k},\varphi_{k})$, so that their positions are characterized only by their solid angle, $\Omega_{k}$. Here, the radius of the projecting sphere, $R$, can be any characteristic dimension of the protein, its circumscribed radius most often being used for this purpose (Hoppe, 2013). The circumscribed radii of the studied proteins are given in Table S3 in the Supporting Material, and fall in the range of $R\sim 1$-$2$ nm.

The surface charge distribution of the discrete charges can then be written simply as

[TABLE]

if rewritten in the form of an expansion in terms of the irreducible spherical representation of multipoles; $\delta(x)$ is the Dirac delta function and $Y_{lm}(\Omega)$ are the spherical harmonics. From eq. (4) we can then obtain the multipole expansion coefficients $\sigma_{lm}$ as

[TABLE]

Contrary to the case of an unscreened Coulomb potential, the Debye screening limits the infinite sum over $l$ in the multipole expansion [eq. (4)] to only the first several terms, leading to a coarse grained description as the details of the charge distribution are smeared below the Debye screening length cut-off. Therefore, even in the presence of higher order symmetries (e.g., octahedral or icosahedral), the lower-order multipoles up to and including the quadrupole provide a good signature of the charge distributions in various molecules, accounting for most of the electrostatic field (Gramada and Bourne, 2006; Platt and Silverman, 1996; Nakamura and Wada, 1985; Kim et al., 2006; Hoppe, 2013).

In view of this, and limiting ourselves to the limit of strong screening, we will focus on the three lowest-order multipoles, the total charge (monopole with rank $\ell=0$), the dipole moment ($\ell=1$), and the quadrupole moment ($\ell=2$). The total charge is independent of the choice of coordinates and can be obtained simply as

[TABLE]

where $k$ runs over all of the charged amino acids in the protein. The total charge as well as the higher-order multipoles will be given in units of elementary charge, $e_{0}$, unless noted otherwise.

For simplicity, in dealing with the dipole and quadrupole moment we will also rescale all our calculations with the characteristic radius of the protein $R$, obtaining an orientational distribution on a unit sphere. The positions of the charges become unit vectors $\bm{n}_{k}$, expressed in spherical coordinates as $\bm{n}_{k}=(1,\vartheta_{k},\varphi_{k})$. As we will also be interested in the orientation of the dipole vector and the eigenvectors of the quadrupole tensor, it will be easier for us to deal with them in Cartesian coordinates. The expressions for the dipole vector $\bm{\mu}$ and the quadrupole tensor $\mathcal{Q}$ are then (Gray and Gubbins, 1984)

[TABLE]

or in an obvious component notation ($\mu_{i}$ and $Q_{ij}$) that we do not write down explicitly. The proper units for the multipoles can be obtained by multiplying the expressions in eqs. (7) and (8) with the corresponding power of the characteristic protein radius $R^{\ell}$ (i.e., $\ell=1$ for the dipole and $\ell=2$ for the quadrupole). Cartesian components of both dipole and quadrupole moments can also be very easily transformed into a spherical form in which they might be more suitable for analytical or numerical calculations (cf. Ref. (Gray and Gubbins, 1984)).

We choose the center-of-mass of each protein for the origin of the coordinate system in which we compute the surface charge distribution and the electrostatic multipoles. Since the monopole moment (the total charge) of the surface charge distribution is always non-zero except at the isoelectric point, the higher multipoles are dependent on the choice of the origin. However, the transformations of the dipole and quadrupole moment to other coordinate systems are very simple and given in Ref. (Gray and Gubbins, 1984).

Dipole and quadrupole principal axes

Surface charge distributions on the proteins, and therefore their multipole expansions, will change with $pH$. This will influence not only the magnitudes of the individual multipoles but also their orientation in space. Consequently, we will be interested in the orientations of the dipole and quadrupole distributions with respect to the original (reference) coordinate system. For each protein, this coordinate system is derived from its PDB entry; the exact orientation of the original coordinate system with respect to the protein structure will not be of interest to us, since we will be interested in relative changes of the dipole and quadrupole orientations with $pH$. While the monopole moment is rotationally invariant, we can always find a rotation of the original coordinate system by keeping the protein structure fixed so that either (i) the dipole vector has a non-zero component only in the $z$ direction, or (ii) the quadrupole tensor is diagonal, with the largest eigenvalue aligned along the $z$ axis. We will refer to these two $z$ axes in the rotated coordinate systems as the dipole and quadrupole principal $z$ axes, respectively.

The dipole moment, being a vector, has three independent components: in the original coordinate system they are the $\mu_{x^{\prime}}$, $\mu_{y^{\prime}}$, and $\mu_{z^{\prime}}$. Upon rotation into the dipole coordinate system, the dipole moment has a non-zero component only along its principal $z$ axis, $\bm{\mu}=(0,0,\mu_{z})$, where $\mu_{z}=(\mu_{x^{\prime}}^{2}+\mu_{y^{\prime}}^{2}+\mu_{z^{\prime}}^{2})^{1/2}$. The two remaining parameters are now the two angles needed to align the original $z^{\prime}$ axis into the new, principal $z$ axis. (Since the only non-zero component is in the $z$ direction, the position of the new $x$ and $y$ axes is irrelevant.)

The quadrupole moment is a symmetric tensor, thus having five independent components. In the quadrupole coordinate system specified by its eigenvectors the tensor becomes diagonal, and we can always order its eigenvalues by value so that the largest one is oriented along the principal $z$ axis, $Q_{xx}\leqslant Q_{yy}\leqslant Q_{zz}$. Since the quadrupole tensor is traceless, we can also express one of its eigenvalues with the other two, e.g., $Q_{yy}=-(Q_{xx}+Q_{zz})$. This leaves us with three more independent components which we can identify with, for instance, the three Euler angles needed to rotate the original coordinate system into the one defined by the quadrupole eigenvectors. Due to the symmetry of the quadrupole, we will also restrict the location of its principal $z$ axis only to the upper hemisphere of the circumscribed sphere.

Quadrupole ratio

In the case of the quadrupole moment we will also be interested in the projections of the surface charge distributions along the $y$ and $x$ axes in the quadrupole coordinate system. Because of our choice of the ordering of the eigenvalues, $Q_{zz}$ will always take on the largest positive value and $Q_{xx}$ the largest negative value. Thus, a ratio of these two eigenvalues, $|Q_{xx}/Q_{zz}|$, will provide us with information on what proportion of the quadrupole distribution is aligned with the principal $z$ and $x$ axis, respectively. (The $y$ axis information can be again omitted due to $\mathrm{Tr}\mathcal{Q}=0$.)

When the quadrupole ratio is close to $0.5$, the quadrupole distribution is axial and oriented predominantly along the $z$ axis. Likewise, when the ratio is close to $2$, the distribution is axial but with an opposite sign and oriented along the $x$ axis. On the other hand, when this ratio is close to $1$, the distribution is represented symmetrically in the $z$ and $x$ directions while vanishing along the $y$ axis. In this scenario, the distribution is not axial but corresponds better to a planar one. The ratio of the two quadrupole eigenvalues thus gives us an insight into the spatial distribution of the quadrupole, which can have an impact on, for instance, the interaction and assembly of molecules with a pronounced quadrupole moment. Some examples of the relation between the quadrupole ratio and the spatial distribution of the quadrupole moment are shown in Figs. S2 and S3 in the Supporting Material.

RESULTS

pH dependence of multipole magnitudes

Figure 2 shows the $pH$ dependence of the first three multipole components in the cases of lysozyme and HSA; similar plots for the MS2 capsid protein and $\beta$-lactoglobulin are shown in Fig. S4 in the Supporting Material. The total charge on the proteins standardly decreases from positive to negative as the $pH$ increases, crossing the point of zero charge at the isoelectric point, $pI$. The isoelectric points of the proteins are interesting since the monopole moment vanishes and higher-order multipoles take on a greater importance. When cysteine acidity is not considered, the predicted isoelectric point of lysozyme is $pI=11.08$, an alkaline $pI$ that is in very good agreement with the values found in the literature ($pI\gtrsim 11$ (Wetter and Deutsch, 1951; Strang, 1984)), even though we used a fairly simple method to obtain it. The $pI$ values we obtain for other proteins are listed in Table 1 and also correspond well with the experimental values: $pI\sim 4.7$–$5.6$ for HSA (Evenson and Deutsch, 1978; Langer et al., 2003) and $pI\sim 5.1$ for $\beta$-lactoglobulin (Majhi et al., 2006; Mercadante et al., 2012). The $pI$ value obtained for the phage MS2 capsid protein at a first glance disagrees with the very low reported value of the phage MS2, $pI\sim 2.2$–$3.9$ (Michen and Graule, 2010; Dika et al., 2015). The notable discrepancy stems most probably from the presence of the genome in the interior of the virion and the permeability of the capsid to external flow (Dika et al., 2011) – in studies where the capsid protein alone was considered, the obtained isoelectric point was at a $pI>8$ (Langlet et al., 2008; Penrod et al., 1996; Schaldach et al., 2006). Our results thus match the ones found in the literature, as we consider only a single capsid protein and not the formed virion (with or without the genome), even though they do not match the experimental results for the formed phage.

In contrast to the total charge, the dipole moment does not change its sign, which is expected, as we always orient the dipole coordinate system so that it points in the $z$ direction. The $pH$ dependence of the dipole magnitude is nonetheless non-monotonic. In a similar fashion, the quadrupole components vary non-monotonically across the range of $pH$ values, with the $Q_{zz}$ and $Q_{xx}$ values always being positive and negative, respectively, due to our definition of their direction. The value of $Q_{yy}$, as a consequence of $\mathrm{Tr}\mathcal{Q}=0$, crosses zero when the $Q_{xx}$ component becomes larger than the $Q_{zz}$ component, or vice versa.

While we use the normalized multipole moments (in units of $[e_{0}/R^{\ell}]$ in the figures to allow for an easier comparison, we also provide in Table 1 the multipole magnitudes in the more intuitive units of $[e_{0}\times\mathrm{nm}^{\ell}]$. We list the dipole and quadrupole magnitudes at the isoelectric point of each protein, as well as the multipole magnitudes at neutral $pH=7$. We can observe a significant decrease in the dipole moment of HSA at its isoelectric point when compared to the dipole moment at neutral $pH$. A closer inspection of Fig. 2 indeed shows that the dipole moment of HSA exhibits a minimum close to its isoelectric point. HSA also attains by far the largest charge and quadrupole moment of the four proteins at neutral $pH$; however, its dipole moment is comparable to that of the others. Another observation we can draw from Table 1 is that dipole moments seem to play a bigger role in the cases of $\beta$-lactoglobulin and phage MS2 capsid protein when compared to their quadrupole moments, while something of the opposite is true for HSA and in a lesser way for lysozyme.

Quadrupole ratio: axial or planar distribution

We have seen that both $Q_{zz}$ and $Q_{xx}$ components of the quadrupole tensor exhibit non-monotonic variation with $pH$. To be able to interpret their relationship more easily we also plot their ratio, $|Q_{xx}/Q_{zz}|$, as a function of $pH$ (Fig. 3). When either of the two components dominates, we have $|Q_{xx}/Q_{zz}|\sim 0.5$ or $|Q_{xx}/Q_{zz}|\sim 2$, and the distribution is axial along the $z$ or $x$ axis, respectively. On the other hand, when the two components are comparable, we have $|Q_{xx}/Q_{zz}|\sim 1$, and the distribution is “planar” in the $x$-$z$ plane. In all the cases studied the quadrupole ratio varies quite a lot, although for most proteins (with the exception of $\beta$-lactoglobulin) it never reaches a fully axial distribution. Conversely, in certain ranges of $pH$ the quadrupole distributions of all the proteins are approximately planar. This is also mirrored by the average values of their quadrupole ratios, which are indeed close to 1 (Table 1).

These results indicate that the nature of quadrupole distributions in the proteins used in our work can be influenced by changing the $pH$: At certain values of $pH$, the quadrupole distributions are axial and thus oriented along a single axis (here, due to our definitions, along the $z$ or $x$ axis). At other values of $pH$, however, the quadrupole distribution becomes planar and symmetrically distributed in the $x$-$z$ plane. The nature of the quadrupole distribution can thus change quite drastically with $pH$ even without any concomitant changes in the conformation of the protein (O’Brien et al., 2011), which should influence, for instance, the interaction properties of these proteins.

pH dependence of dipole and quadrupole principal axes

Now that we have seen how the magnitudes of different multipole moments change with $pH$, we turn our attention to the $pH$ dependence of the orientations of the dipole and quadrupole principal axes when compared to the original (reference) coordinate system. In order to indicate the level of detail described by different multipoles, we show in Figs. 4a and 4b the multipole expansions of the surface charge distribution of $\beta$-lactoglobulin on a fine-grained level with maximum rank of $\ell=6$ and on a coarse-grained level with rank $\ell=2$, respectively. The distributions are mapped from a sphere to a plane using the Mollweide projection, which has the polar and azimuthal angles as coordinate axes (Snyder, 1987). Separately, we also show in Figs. 4c and 4d both the dipole and quadrupole distributions in the original coordinate system. All the plots are shown at $pH=7$, and a similar figure for HSA is shown in Fig. S5 in the Supporting Material.

From the multipole representation of the surface charge distribution where terms up to the order $\ell=6$ are included (Fig. 4a) it is not immediately obvious where the dipole and quadrupole axes are located. However, when we isolate both terms (Figs. 4c and 4d), this becomes more apparent. When the two distributions are combined together with the total charge (Fig. 4b), they describe the coarse-grained variation of the surface charge consistent with the Debye screening cut-off. Similar observations can be drawn also for the other proteins studied. For example, by comparing Figs. 4 and S5 we can observe a noticeable difference between $\beta$-lactoglobulin and HSA: In the former, the dipole distribution is dominant among the lower-order multipoles, while the exact opposite is true in the case of the latter.

To demonstrate next how the orientations of the principal $z$ axes of the dipole and quadrupole moments change with $pH$, we show in Fig. 5 snapshots of the multipole representation of the surface charge distribution of lysozyme (with $\ell_{\max}=6$) at three different values of $pH$. Similar snapshots for HSA are shown in Fig. S6 in the Supporting Material. For a complete comparison, Fig. 6 isolates the orientations of the dipole and quadrupole principal $z$ axes of all four proteins studied, and shows their variation in space over the entire range of $pH$ values while keeping the coordinate system of the protein fixed. (Note that the apparent jumps in the location of the quadrupole axis in the cases of HSA and $\beta$-lactoglobulin are a consequence of our confinement of the axis to the upper hemisphere for reasons of symmetry.)

As $pH$ is increased, the overall charge moves towards more negative values. At the same time, the positions of the dipole and quadrupole principal $z$ axes trace quite a path in space. In the case of the lysozyme (Fig. 6a), most of the variation in the axes’ orientation occurs after $pH>7$. On the contrary, in the case of HSA (Fig. 6b) the majority of the variation, especially large in the case of the dipole, happens up until that point, i.e., when $pH<7$. It is worth noting that the isoelectric points of the two proteins are at the opposite sides of the spectrum (Table 1). The axes of the MS2 capsid protein exhibit perhaps the least variation (Fig. 6c), while the quadrupole axis of $\beta$-lactoglobulin shifts in space to a large extent (Fig. 6d).

Taken together, we see that not only does $pH$ influence the magnitudes of the various multipole moments, but it influences even more so their orientation in space. The changes in the orientation do not happen uniformly with the changing $pH$, but are spaced unevenly across the range. Most of the variation seems to usually (but not always) occur in the $pH$ range near the isoelectric point of each protein.

DISCUSSION

The results of our model clearly show that $pH$ variation causes significant changes not only in the magnitudes of charge multipole moments in various proteins, but also in the orientations of the principal axes of their multipoles. In obtaining these results, we have, however, resorted to some necessary simplifications. First of all, while we have used the acid-base dissociation constants $pK_{a}$ of amino acid functional groups in bulk dilute aqueous solutions, the actual $pK_{a}$ values in proteins are site-dependent, i.e., influenced by the local ionic and structural environment (Hass and Mulder, 2015; Schaefer et al., 1997). The deviations range from small to significant, depending on the type of site considered (Gitlin et al., 2006; Langlet et al., 2008). Secondly, changes in $pH$ often also lead to conformational changes in proteins (O’Brien et al., 2011; Schaefer et al., 1997; Harrison et al., 2013; Goh et al., 2014), whereas we have used fixed structures obtained from PDB to extract the protein charge distributions. In the following, we discuss these simplifications and limitations of our model, and show that our conclusions should remain valid even as more details are included in the description of charge distributions in proteins.

In our study, we have intentionally de-coupled the effects of conformational changes from the $pH$-dependent changes in charges of a protein, as this allowed us to clearly observe and study the effects of $pH$ on the charge multipoles. In this way we were also able to treat the electrostatics of the proteins using a simple analytical model. Coupling the electrostatics with the conformational changes would, on the other hand, necessarily demand a different treatment – involving, e.g., studying protein dynamics in order to predict the changes in conformation. At the same time, joining the two models would make it harder to discern which changes in the magnitudes and orientations of the charge multipoles are a consequence of the conformational changes of a protein, and which are a consequence of its $pH$-dependent fractional charges. In the remaining discussion we will thus not delve any further into the influence of conformational changes, but will nonetheless examine the results of our model in light of the effects of local $pK_{a}$ shifts arising as a consequence of site-site interactions or electrostatic potential. We stress, however, that in real systems where both $pH$-dependent charge and conformational changes operate in parallel, we can expect some quantitative changes in the observed $pH$-dependent behavior of the protein charge multipoles. The discrepancy will naturally be larger in those proteins where conformational changes due to $pH$ are more significant.

To investigate the basic effect of the local $pK_{a}$ shifts induced by site-site interactions, we used PROPKA software for the empirical prediction of $pK_{a}$ values of the ionizable amino acids, as described in the Materials and Methods. The results are summarized in Figs. S7, S8, and S9, and in Table S4 in the Supporting Material. Even though the local $pK_{a}$ values of the amino acid residues now differ from the bulk $pK_{a}$ values used in our model (Table S1) upward to $2$ units of $pH$, the effect this has on the magnitudes of multipole moments is nonetheless small (Fig. S7). The biggest change can be observed in the case of the phage MS2 capsid protein in the range around neutral $pH$; however, the difference in the predicted isoelectric point is still less than $1$ unit of $pH$ (Table S4). Our model with $pK_{a}$ values derived from PROPKA also predicts less of a difference between the magnitudes of the dipole and quadrupole moments in the case of $\beta$-lactoglobulin, yet these differences are at most within a factor of $2$ from the predictions of the basic model. Similar differences can be observed with PROPKA-predicted $pK_{a}$ values in the orientations of the multipole moments, which can range from small in the case of lysozyme (Fig. S8) to bigger in the case of $\beta$-lactoglobulin (Fig. S9), yet never significantly change the qualitative behavior of the $pH$ dependence of the multipole orientations.

Another factor that we have not explicitly considered in our model and can nevertheless potentially have a significant effect on the local $pK_{a}$ values is the electrostatic potential. In principle, the local electrostatic potential at the position of each ionizable amino acids modifies its $pK_{a}$ value [eq. (3)]. Only when the potential is small is the local $pK_{a}$ shift due to it negligible. This holds, for instance, in the limit of high salt or large screening, relevant in general for physiological conditions. In order to estimate the $pK_{a}$ shift due to electrostatic potential in this limit we have extended our model by solving the Debye-Hückel equation for the electrostatic potential of the protein in the presence of charge regulation boundary condition. The details of this methodology are given in the Materials and Methods and in the Supporting Material. We observe that the average $pK_{a}$ shift imparted by the electrostatic potential is on the scale of $1$ $pH$ unit in the limit of high salt concentration ($c_{0}\sim 1$ M), as shown in Fig. S1; this is comparable with the $pK_{a}$ shifts predicted by PROPKA. In contrast to the latter, the shifts caused by local electrostatic potential are also $pH$ dependent due to the nature of the model. In the limit of high salt, the electrostatically-shifted $pK_{a}$ values have only a small effect: The predicted isoelectric points are within less than $1$ unit of $pH$ from the predictions of our basic model, and the $pH$-dependence of magnitudes and orientations of the multipole moments follow similar patterns as in the case where electrostatics is not included (Table S4 and Figs. S7, S8, and S9). Notably, the predictions based on the inclusion of the electrostatic potential are sometimes closer to those obtained using PROPKA-predicted $pK_{a}$ values, while in other cases closer to those obtained using our simple model.

It has to be noted that there are also limitations to the electrostatic model we have used for the prediction of local $pK_{a}$ shifts. The DH equation performs well in the limit of high salt, which was of interest to us. However, it is known to overpredict the values of the electrostatic potential in the limit of vanishing salt. More importantly, the DH equation is based on a continuum description of the electrostatics. While this is in line with the general model used in our present work, any ion-specific binding effects are thus neglected. Coupling the local description of charge interaction with a global multipole expansion of charge distribution is a difficult problem due to the difference of scales involved. A consistent inclusion of local and ion-specific effects would at this point require a drastically different model, based on an atomic description of the protein and the solvent, while our approach still retains its fundamental validity in the high salt limit.

Although we have not considered detailed structural environments of proteins in our work, we note that structural changes bring about a different configuration of charges on the surface of a protein. A similar re-configuration of charges can arise by the variation of the RSA cut-off $c$. To show that the choice of the RSA cut-off does not have a large impact on the qualitative results of our study, we plot in Fig. S10 in the Supporting Material the $pH$ dependence of the multipole moments of the four proteins studied for three different choices of the cut-off, ranging from $c=0.1$ to $c=0.5$. In each of the cases considered, going from $c=0.1$ to $c=0.5$ amounts to quite a significant loss of $24$–$34$% of the total number of surface charges (see Table S1). And yet, some apparent differences when the cut-off is increased notwithstanding, the qualitative behavior of the multipole moments remains much the same. This validates our approach, which fixed the cut-off to $c=0.25$, and the conclusions drawn from it: While small variations of the number of charges on the surface of a protein will necessarily change the underlying multipole expansion of the surface charge distribution, the qualitative behavior we observed when the $pH$ is varied will remain unaltered.

Lastly, we separately examine the possibility of cysteine protonation, which was not explicitly considered in our main model. Due to its $pK_{a}$ value, this should in principle push the charge distribution towards more negative values at $pH>7$. In order to see how the consideration of cysteine acidity influences our results, we show in Fig. S11 in the Supporting Material a comparison of the $pH$ dependence of the multipole magnitudes in the case of HSA and MS2 capsid protein. We can again observe that the presence of cysteine protonation has no qualitative effect on the behavior of the system. It does affect the location of the isoelectric point to some extent (cf. Table S4), although the difference is again less than 1 unit of $pH$. Other than that, the magnitudes of the multipole components and their $pH$ dependence remain relatively unchanged. Similarly, by comparing the orientations of the dipole and quadrupole axes between the cases where cysteine acidity is or is not considered (Figs. S12 and 6, respectively) we see that the presence of cysteine charges influences the orientational variation of the principal axes at $pH>7$, as expected, although the changes are usually minor.

CONCLUSIONS

In this work, we have studied the effects of $pH$ on charge multipole moments in proteins. We first obtained $pH$-dependent surface charge distributions of four globular proteins and expanded the distributions in terms of electrostatic multipoles. We have limited ourselves to lower-order multipoles (the total charge and the dipole and quadrupole moments), and studied the effect of $pH$ on both their magnitudes and the orientations of their principal axes in space. The value of $pH$ was shown to have a significant effect on both, particularly on the orientations of the multipole principal axes, which in some cases exhibited large variation. This variation was found to be a non-uniform function of $pH$, spaced rather unevenly across the solid angle, with most of the changes occurring in the $pH$ range near the isoelectric point of each protein.

We have also pointed out some limitations and necessary simplifications of our approach. While the precise determination of the fractional charge of amino acid residues in a protein due to site-site interactions and local electrostatic potential is difficult, we have nonetheless been able to explore the general influence of local $pK_{a}$ shifts. We have done this by extending our simple model with either empirically predicted $pK_{a}$ values or with $pK_{a}$ values shifted by an electrostatic potential derived from a detailed solution of the corresponding DH equation. In this way, we have shown when and how the qualitative conclusions derived in our paper hold even when the surface charges vary due to local $pK_{a}$ shifts, and have demonstrated that the observed behavior of multipoles described by our model should persist even as a more detailed picture becomes attainable.

The highly protein-specific $pH$-dependent changes in charge distributions are often intertwined with $pH$-dependent conformational changes of proteins and can occur concomitantly. In order to elucidate the direct influence of the $pH$ dependence of charge multipoles using a simple analytical model, we have neglected any conformational changes; this is, of course, an oversimplification. However, the main results of our model persisted even as the number of charges or their local $pK_{a}$ values were varied. This indicates that the $pH$ variation of the charge multipole moments of proteins can possibly also play a role as a driving force for conformational changes, and that the $pH$ variation of the direction of the multipole principal axes can affect the electrostatic part of the deformation energy. In this way, it can induce an orientational conformational change of the protein (O’Brien et al., 2011; Schaefer et al., 1997), in direct analogy to the structural phase transitions in general condensed matter context.

In addition, the orientational variation of the multipole principal axes remains of fundamental importance for protein-protein interactions also when screened electrostatic interactions in the presence of ionic solution are taken into account. There, the multipole moments play a significantly different role compared to the case of pure, unscreened Couloumb interactions, in the sense that the screened electrostatic potential retains the full directional dependence of all multipole moments even in the far-field limit. Consequently, the usual argument that at large separations only the monopole moment matters is untenable, and while this property of screened multipole electrostatics is not new (Rowan et al., 2000; Kjellander and Ramirez, 2008; Lošdorfer Božič and Podgornik, 2013; Kjellander, 2016), it has yet to penetrate the community involved in modelling protein-protein interactions.

Finally, we can speculate that higher-order multipole moments can provide a specific electrostatic signature of each protein, and could potentially be used to classify different protein types (Nakamura and Wada, 1985). This electrostatic signature could be further refined by defining and tracking the “orientations” of higher-order multipoles, possibly allowing for a classification scheme based on a small number of significant parameters coding for the protein symmetry and interactions.

The major lesson of our investigation is thus that $pH$ variation has a significant influence on the magnitudes and especially orientations of charge multipole moments in proteins. This finding was derived using a simple analytical model in a form compatible with models of charged patchy proteins, and at the same time amenable to the inclusion of additional effects, such as empirically predicted $pK_{a}$ values or $pK_{a}$ shifts due to local electrostatic potential. The results presented in this work should be of particular significance to protein assembly engineering, e.g., for viral capsids and enzyme nanocontainers, where $pH$ could drive symmetry transitions of the entire assembly induced by local, orientationally-dependent electrostatic interactions.

SUPPORTING MATERIAL

Four tables, twelve figures, and thirteen equations of supporting material are available.

AUTHOR CONTRIBUTIONS

A.L.B. and R.P. designed the research; A.L.B. performed the research; A.L.B. and R.P. analyzed the data and wrote the article.

ACKNOWLEDGMENTS

We thank S. Čopar, G. Posnjak, and I. Rapošová for helpful discussions and comments on the manuscript. We are also grateful to the anonymous reviewers for their comments and suggestions to improve the quality of the manuscript. A.L.B. and R.P. acknowledge the financial support from the Slovenian Research Agency (research core funding No. (P1-0055)).

I NUMBER OF SURFACE CHARGES AS A FUNCTION OF RSA CUT-OFF

II ACID-BASE EQUILIBRIUM CONSTANTS

III CIRCUMSCRIBED RADII OF STUDIED PROTEINS

IV ROLE OF ELECTROSTATIC POTENTIAL IN CHARGE REGULATION

Local electrostatic potential $\psi(\mathbf{r})$ can play a role in the regulation of charge on the amino acid residues, inducing a local $pK_{a}$ shift. The effect of the potential on the charge of an amino acid can be written as Ninham and Parsegian (1971)

[TABLE]

where $\beta=1/k_{B}T$, $T$ is the room temperature, $k_{B}$ the Boltzmann constant, and $e_{0}$ the elementary charge. The renormalized $pK_{a}$ values thus become

[TABLE]

In order to obtain the electrostatic potential, we use a reformulation of the canonical Tanford-Kirkwood model (Tanford and Kirkwood, 1957; Barroso da Silva et al., 2001) by treating the protein as a sphere with radius $R$ and dielectric constant $\varepsilon_{p}=4$, carrying a surface charge distribution $\sigma(\Omega)$ including the monopole, dipole, and quadrupole contributions ($\ell_{\mathrm{max}}=2$), and immersed in a $1:1$ salt solution with bulk concentration $c_{0}$ and with the dielectric constant of water $\varepsilon_{w}=80$. At large, physiologically relevant salt concentrations ($c_{0}\gtrsim 100$ mM), we can solve the Poisson-Boltzmann equation for the electrostatic potential in the linearized, Debye-Hückel (DH) approximation:

[TABLE]

inside and outside of the protein, respectively (for details, see for instance Ref. Lošdorfer Božič and Podgornik (2013)). Here, $\kappa=\sqrt{8\pi\ell_{B}c_{0}}$ is the inverse Debye length, and $\ell_{B}=\beta e_{0}^{2}/4\pi\varepsilon_{0}\varepsilon_{w}$ is the Bjerrum length. The two potentials are continuous on the boundary of the protein,

[TABLE]

while the jump in their derivatives is proportional to the charge on the protein, which is in turn coupled with the local potential

[TABLE]

Here, we have used the high-salt/low electrostatic potential approximation to expand the regulated charge in eq. (S1) to the linear order in $\beta e_{0}\psi\ll 1$:

[TABLE]

This in turn yields a surface charge expansion of the form

[TABLE]

where both terms can be expanded in terms of multipoles,

[TABLE]

We note that the first term of the expanded surface charge density, $\sigma^{(0)}(\Omega)$, corresponds to the original surface charge density used in the main text which is not corrected for the effects of electrostatic charge regulation.

The DH equation for the electrostatic potential inside and outside the protein [eq. (S3)] can be solved in the form

[TABLE]

where $k_{l}(x)$ are the modified spherical Bessel functions of the second kind. By equating the two potentials at $r=R$ and taking into account the boundary condition for the derivative of the potential at $r=R$, we obtain the following expression for the coefficients $B_{uv}$:

[TABLE]

Here, we have introduced a shorthand

[TABLE]

and used the standard form of the Wigner 3-j symbols (for details, see again Ref. Lošdorfer Božič and Podgornik (2013))

[TABLE]

Due to charge regulation on the surface of the protein, the different coefficients $B_{uv}$ are clearly coupled, and eq. (S10) cannot be solved further analytically. However, an iterative scheme with the initial value $B_{uv}^{0}=\sigma_{uv}^{(0)}/f_{u}(\kappa R)$ allows us to obtain the coefficients $B_{uv}$ with a desired accuracy in typically 5-10 iterations.

The solution for the coefficients $B_{lm}$ allows us to obtain the local electrostatic potential on the surface of the protein in the continuum approximation. Evaluating the potential at the positions of the amino acid residues enables us to write their re-normalized charges as given by eq. (S1) and to obtain the $pK_{a}$ shifts of individual amino acids [eq. (S2)]. We also note that the $pK_{a}$ shifts obtained in this way are $pH$-dependent due to the way the electrostatic potential is coupled with the fractional charge of the amino acids. In order to show the average effect the electrostatic potential has on the local $pK_{a}$ values, we plot in Fig. S1 the mean $pK_{a}$ shift, averaged over all $N$ ionizable amino acids in a protein,

[TABLE]

V GEOMETRICAL INTERPRETATION OF THE QUADRUPOLE EIGENVALUES RATIO

VI PH DEPENDENCE OF MULTIPOLE MAGNITUDES

VII PH DEPENDENCE OF DIPOLE AND QUADRUPOLE PRINCIPAL AXES

VIII COMPARISON WITH PROPKA AND ELECTROSTATIC CHARGE REGULATION

IX EFFECT OF RSA CUT-OFF VARIATION

X EFFECT OF CYSTEINE PROTONATION

Cysteine has a thiol functional end group which is a very weak acid; in addition, it is reactive and can form disulfide bonds. For these reasons it is often not considered as an acid at all, and we thus did not take cysteine protonation into account in the main text. Here, we consider the effects in the case when the acidity of cysteine is considered. Cysteine residues comprise between $2$–$14$% of our studied proteins’ AA composition at the RSA cut-off of $c=0.25$ (cf. also Table S1). Due to its $pK_{a}$ value (Table S2), cysteine’s region of influence is limited to $pH>7$ and its presence shifts charge distributions towards more negative values in that range (Nap et al., 2014).