pH environments

A protein never meets "pH 7". Pepsin works in a stomach at pH 1.5–2.5, a peptide transporter of the gut wall faces an acid film of about pH 6 on its outer side and the cytosol's pH 7.2 on its inner side, an antibody circulates at 7.4, and a lysosomal enzyme lives at 4.7. Each of these places sets the charge of every acidic and basic group of the protein, and the charges in turn set how the protein binds, folds, sticks to membranes and falls apart.

GL1X answers this for the structure of any MolNFT record. For every ionizable group of the protein it estimates the pKa from the structure: the pH at which the group holds its proton half of the time. From the pKa values it gives the group's charge in each environment you choose, the protein's net charge and isoelectric point, how its folding stability changes with pH, the groups that switch charge from one environment to another, and the pH at which the structure itself was determined. For a membrane protein it can place the membrane and give the outer and inner sides their own pH, so that a transporter's outer face can meet the gut lumen while its inner face meets the cytosol.

Where: the pH environments tool of the research workbench (/bio/workbench/{contract}/{token}?tool=ph, also in the left menu under Structure), the API (/api/bio/records/{contract}/{tokenId}/ph, /api/bio/ph/environments) and MCP (ph_environment).

These are estimates on fixed coordinates, made by an empirical model. They are good at telling which groups are unusual and why, and at comparing environments; they are not measurements. The section How the model was set and checked says how close the model comes to an established program, and Limits says what it leaves out.

What the analysis answers #

The chemistry in brief #

pKa, charge and the Henderson–Hasselbalch equation #

An acid such as the carboxyl of aspartate gives up a proton and becomes negatively charged (COOH → COO⁻ + H⁺); a base such as the amine of lysine takes up a proton and becomes positively charged (NH₂ + H⁺ → NH₃⁺). Which form dominates depends on the proton concentration, written as pH = −log₁₀[H⁺], and on the group's pKa. For a group on its own, the fraction that holds its proton follows the Henderson–Hasselbalch equation:

fraction protonated  f = 1 / (1 + 10^(pH − pKa))
charge of an acid    q = f − 1        (−1 when it has lost its proton)
charge of a base     q = f            (+1 when it holds its proton)

One pH unit away from the pKa the group is 91% in one form; two units away, 99%. An aspartate with a pKa of 3.67 holds its proton 99.3% of the time at pH 1.5 (charge −0.007) and has lost it 99.98% of the time at pH 7.4 (charge −1.00). A histidine with a pKa of 6.54 is 12% charged at 7.4, 73% at 6.1 and 99% at 4.7: it is the group that switches between blood, the intestinal wall and the lysosome.

A pKa is a free energy in disguise. The free energy of taking the proton away is ΔG = kT·ln10·(pKa − pH): at 25 °C, kT·ln10 = 1.364 kcal/mol, so every interaction that changes a group's free energy by 1.364 kcal/mol moves its pKa by one unit. This is how GL1X turns energies into pKa shifts and back.

Why a protein moves its groups' pKa #

In a small unfolded peptide each group sees only water, and its pKa is close to a fixed value for its kind: the model pKa. Folding changes the surroundings, and with them the free energy of the charged form:

• Burial (desolvation). Water stabilises a charge strongly. A charge pushed into the protein loses part of that water and pays for it, so the neutral form gains: buried acids keep their proton longer (their pKa rises), buried bases lose theirs sooner (their pKa falls).
• Hydrogen bonds. A backbone amide donating a hydrogen bond to a carboxylate stabilises the negative form, so the acid's pKa falls; a carbonyl accepting one from an ammonium stabilises the positive form, so the base's pKa rises.
• Other charges. A positive charge nearby favours a neighbouring acid's negative form (its pKa falls) and disfavours a neighbouring base's positive form (its pKa falls too); a negative charge does the opposite.

Measured pKa values in folded proteins spread widely around the model values: Grimsley et al. (2009) summarised 541 of them in 78 proteins and found Asp 3.5 ± 1.2, Glu 4.2 ± 0.9, His 6.6 ± 1.0, Cys 6.8 ± 2.7, Tyr 10.3 ± 1.2 and Lys 10.5 ± 1.1 (mean ± standard deviation). The spread is what a structure-based estimate tries to explain.

Coupled groups #

Charges interact with each other, so groups that sit close together do not titrate independently: whether one group holds its proton depends on whether its neighbour does. Two carboxylates 4 Å apart in a buried site repel each other when both are charged, so the first loses its proton at a low pH and the second only at a much higher one. Their titration curves spread out, and the Hill coefficient of each falls below 1. GL1X therefore titrates the groups together rather than one at a time (see Titration of coupled groups).

The environments #

GL1X offers the compartments below by name, each with the pH its source reports. Any pH from 0 to 14 can be given as a number instead, and the pH the structure was determined at can be added as the environment experiment (see The structure's own pH).

The digestive-tract values come from radiotelemetry capsules swallowed by healthy fasting adults (Evans et al. 1988): gastric pH ranged from 1.0 to 2.5 in every subject, the proximal small intestine read 6.6 during the first hour, pH rose to 7.5 in the terminal ileum, fell to 6.4 in the caecum and rose again to 7.0 in the left colon. The brush-border value is different in kind: it is the pH of the thin unstirred layer of fluid at the surface of the intestinal wall, measured with microelectrodes in rat jejunum (about 6.1 there, 6.5 in the midgut and 7.3 in the distal ileum; Lucas 1983). This acid microclimate is the proton gradient that proton-coupled transporters of the gut wall use to take up peptides, nutrients and drugs (Thwaites and Anderson 2007), and it is what the outer face of such a transporter meets.

Every environment is a single pH. Real compartments vary (a meal raises gastric pH, disease changes many of them), and the reported spread is there to show by how much; analyse the ends of the range as separate environments when it matters.

How GL1X estimates a pKa #

For each ionizable group the estimate is built in two stages. The intrinsic pKa collects everything that does not depend on the other groups' states:

pKa_intrinsic = pKa_model + ΔpK_desolvation + ΔpK_hydrogen_bonds + ΔpK_fixed_charges

The groups are then titrated together: their mutual interactions depend on which of them are charged at a given pH, so the apparent pKa of each group comes out of the joint titration, and the difference between the apparent and the intrinsic pKa is reported as its coupling shift.

The groups that titrate #

Residues are recognised by name, including the protonation-state names that simulation files use (ASH, GLH, HID, HIE, HSD, HSE, HSP, CYX, CYM, LYN, ARN); a residue named HIP that holds a phosphorus atom is the PDB's phosphohistidine, not a protonated histidine, and is left out. Every protein chain gets one N-terminal amine on its first residue and one C-terminal carboxyl on its last, unless the chain starts with an acetyl, formyl, pyroglutamate or formylmethionine group (ACE, FOR, PCA, FME) or ends with an amide cap (NH2, NME). When a chain's first or last residues are missing from the file (disordered tails), the terminus sits on the first or last residue the file has.

Some groups do not titrate:

• Cysteines in disulfides (two SG atoms within 2.5 Å) have no proton to lose. They are counted under not_titrated.
• Groups bonded to a ligand (a charged atom within 2.1 Å of a ligand's non-metal atom: a cysteine linked to a heme's vinyl group, a covalent inhibitor) are left out.
• Groups coordinating a metal (a charged atom within 2.8 Å of a metal atom, in an ion or in a cofactor such as heme) keep a fixed state: a coordinating carboxylate, thiolate or phenolate is deprotonated (−1), a coordinating histidine or amine is neutral. Their pKa is not reported.
• Groups whose atoms are missing (a side chain the crystallographers could not see) keep their model pKa and titrate alone: they are disordered and exposed to water, so their charge still counts in the totals.

Modified residues (phosphoserine, methyllysines, oxidised cysteines) are not titrated, nucleotide bases are not titrated, and ligands other than metal ions carry no charge (see Limits).

Model pKa values #

The model pKa of each kind is its pKa in a small unfolded peptide, where it sees only water. GL1X takes the values Thurlkill et al. (2006) measured in alanine pentapeptides (0.1 M KCl, 25 °C) and, for arginine, the value Fitch et al. (2015) measured for the guanidinium group:

These are not the same values the sequence tools use: the theoretical pI of Molecular features and the virtual gel use the pK set of Bjellqvist et al., which was fitted to where proteins stop in isoelectric focusing gels. The two sets answer different questions and give somewhat different pI values for the same sequence.

The temperature of the model is 25 °C, the temperature at which the model values were measured. Groups with large ionization enthalpies (histidine and the amines) have lower pKa values at body temperature; GL1X does not correct for it.

Burial and desolvation #

Burial is measured by counting atoms. Around each group's charge centre GL1X counts the heavy atoms of every other residue and ligand (waters left out) within 15 Å, and maps the count N onto a burial fraction:

burial B = 0                     when N ≤ 200
         = (N − 200) / 600       between
         = 1                     when N ≥ 800

Ionizable groups mostly sit near the surface: over the 5,054 groups of the 76 benchmark structures compared with PROPKA below, the median count was 285 atoms, the 10th percentile 159, the 90th 499 and the largest 718, so a count of 800 is deeper than any group in those proteins. The count looks far beyond the group's immediate neighbours on purpose: a charge is screened by the water within several ångströms of it, and how much water a deep cleft or a domain interface removes depends on the protein's shape on that scale.

The desolvation shift grows linearly with burial up to a largest value that depends on the kind:

ΔpK_desolvation = ± D_kind × B        (+ for acids, − for bases)

The values were set by calibration (see How the model was set and checked): they are the sizes that best reproduce the shifts of the reference program, a few pK units for a fully buried group. Arginine, whose charge is spread over three nitrogens, pays less to be buried, which its small value expresses; the C-terminus is almost always exposed, which leaves its value poorly determined.

Hydrogen bonds #

GL1X looks for hydrogen bonds between each group and the atoms around it that do not titrate:

The backbone hydrogen is placed as DSSP places it (Kabsch and Sander 1983): 1 Å from the nitrogen, opposite the previous residue's carbonyl. How fully a partner bonds a group atom is a number from 0 to 1:

backbone N–H:  1 when H···A ≤ 2.1 Å, 0 when ≥ 2.7 Å, linear between
               (2.5 and 3.1 Å when the acceptor is a sulfur),
               and only when the angle N–H···A is at least 110°
side chains:   1 when the heavy atoms are ≤ 2.8 Å apart, 0 when ≥ 3.5 Å
               (3.3 and 3.9 Å with a sulfur)

Each full bond shifts the pKa by a fixed amount:

A partner counts once, by its best contact with any of the group's atoms; the group's own residue is not a partner, and the total is held between −3 and +3 pK units. Hydrogen bonds between two titrating groups (a salt bridge between a carboxylate and an ammonium) are not counted here: their effect depends on both groups' states and belongs to the charge–charge interactions. Tyrosine is the exception: its OH is a partner, taken as the neutral phenol it is below pH 10, so a tyrosine hydrogen-bonded to a carboxylate shifts it here and the two groups also interact as charges.

Charges that do not titrate #

Three kinds of fixed charge act on the groups:

• Metal ions: every single-atom metal ion of the structure (an element on the ligand index's list of metal elements), with the charge its chemical component's formula states (Zn²⁺ from "Zn +2"), or, when the file gives none, the charge of the element's usual simple cation in water: +1 for the alkali metals, silver, gold and thallium; +2 for the alkaline-earth metals, manganese, iron, cobalt, nickel, copper, zinc, cadmium, mercury, lead, tin, palladium and platinum; +3 for aluminium, gallium, indium, scandium, yttrium, the lanthanides, chromium, bismuth, rhodium, actinium and the actinides from americium on; +4 for zirconium, hafnium and thorium. The metalloids (germanium, arsenic, antimony), polonium, and the metals whose usual forms in water are oxo ions or several (titanium, vanadium, niobium, tantalum, molybdenum, tungsten, technetium, rhenium, ruthenium, osmium, iridium, protactinium, uranium, neptunium, plutonium) get no default: they still hold the groups they coordinate, but add no charge. The full list is in Names and identifiers. A metal inside a cluster or a cofactor (the iron of a heme, the irons of an iron–sulfur cluster) gets no point charge, because its net charge is not known.
• Nucleic acid phosphates: one negative charge on each phosphorus of a DNA or RNA chain in the structure, so that the arginines and lysines of a DNA-binding protein feel the DNA.
• Coordinating groups: the fixed charge of each group held by a metal (above).

A fixed charge q at distance r shifts a group's pKa by

ΔpK = −q × W(r) / (kT·ln10)

with W the interaction of two unit charges defined in the next section. The total over all fixed charges within 20 Å is the group's fixed_charges shift; each one that moves it by 0.05 or more is listed among its partners.

Charge–charge interactions #

Two unit charges at a distance r interact with the energy

W(r) = 332.06 / (ε × r) × exp(−r × (1 − b) / λ_D)      kcal/mol

where 332.06 kcal·Å/mol converts the charges' Coulomb energy to chemical units, r is the distance between the groups' charge centres (at least 3 Å, pairs farther than 20 Å apart are dropped), b = √(B₁·B₂) is the pair's mean burial, ε is the dielectric between them and λ_D the Debye length of the solution.

The dielectric expresses how much the medium between two charges weakens their interaction. In water it is 78.4 at 25 °C; between two charges buried in a protein it is far smaller when they are close, since little water fits between them, and rises towards water's within several ångströms as screening material returns. GL1X writes this as

ε(r, b) = 78.4 − (78.4 − ε_core) × S(r / 3 Å)
ε_core  = 78.4 − (78.4 − 12) × b
S(R)    = (R² + 2R + 2) × e^(−R) / 2          (1 at r = 0, 0.92 at 3 Å, falling to 0)

so an exposed pair (b = 0) always sees water's 78.4. For a fully buried pair (b = 1), ε_core = 12 is the limit as r goes to zero; at the 3 Å the distances are floored at it sees 17.3, at 4, 6, 8 and 12 Å it sees 22.0, 33.5, 45.1 and 62.6, and water's 78.4 beyond. The value 12 was set by calibration.

The ions of the solution screen charges further. The Debye length is λ_D = 3.04/√I Å for an ionic strength I in mol/L at 25 °C (Debye and Hückel 1923): 7.85 Å at the physiological 0.15 M. The screening applies to the part of the pair exposed to the solution (the factor 1 − b), since ions do not enter the protein.

What this gives, in pK units of shift, for two opposite charges at 0.15 M:

A salt bridge on the surface is worth about half a pK unit; the same pair buried is worth several, which is why buried salt bridges are rare and important. The desolvation of the two charges works in the other direction, and the balance decides.

Titration of coupled groups #

With N groups, each holding its proton (x = 1) or not (x = 0), a protonation state of the whole protein has the free energy, in units of kT,

E(x) = Σᵢ xᵢ · ln10 · (pHᵢ − pKa_intrinsic,ᵢ) + Σᵢ<ⱼ wᵢⱼ · qᵢ(xᵢ) · qⱼ(xⱼ)

where qᵢ is group i's charge in that state (x − 1 for an acid, x for a base), wᵢⱼ = W(rᵢⱼ)/kT the interaction of the two groups, and pHᵢ the pH group i sees (the same for every group, except in membrane mode). The probability of a state is proportional to exp(−E), and each group's protonated fraction is its average over all states. Exact averages need all 2^N states, which is out of reach beyond about 25 groups, while a protein has dozens to thousands.

GL1X uses the idea of the hybrid method of Gilson (1993): groups that interact strongly are bound into clusters, and the weak interactions between clusters are treated in mean field.

1. Every pair of groups that interact by 1 kT or more (0.59 kcal/mol, 0.43 pK units) is a candidate bond. Bonds are taken from the strongest down, and a bond joins two clusters only while the joined cluster has at most 8 groups.
2. Within a cluster every joint state is enumerated (at most 2⁸ = 256), each state's energy following from the last one by flipping one group (states are visited in Gray-code order).
3. Each cluster sees the groups of the other clusters through their mean charges: the field φᵢ = Σⱼ wᵢⱼ⟨qⱼ⟩ over groups j outside the cluster.
4. The clusters are averaged in turn and the mean charges updated (each new value moved 60% of the way from the old one, which keeps the iteration stable) until no protonated fraction changes by more than 10⁻⁷.

The titration runs from pH 0 to 14 in steps of 0.1 (0.25 for structures with more than 3,000 groups), each pH starting from the solution at the previous one. On the 92 readable structure files of the benchmark folder (6,878 groups) the pKa values agree with a nearly exact reference (clusters of up to 12 groups bound by couplings of 0.5 kT) within 0.004 pK units root-mean-square and 0.11 at worst, at about a twenty-fifth of the reference's cost; a unit test checks the cluster solution against the exact average over all states on small random systems.

The apparent pKa and the Hill coefficient #

A group's apparent pKa is the pH at which its titration curve crosses half protonation, interpolated between the grid points. A curve that never crosses within 0–14 is reported as below or above the range (arginines, whose model value is 13.8, are usually above it). Strong coupling can make a curve cross half protonation more than once; the group is then marked irregular, and the crossing nearest the curve's midpoint by area (pH₀ + ∫ f dpH) is kept.

The Hill coefficient is the slope of the curve at its pKa:

n = − d log₁₀( f / (1 − f) ) / dpH

It is 1 for a group that titrates alone, below 1 when neighbours of the same charge spread its titration out, and above 1 when they make it cooperative.

Each group's answer separates the contributions: model_pka, the desolvation, hydrogen_bonds and fixed_charges shifts, the intrinsic_pka they add up to, the coupling shift (apparent − intrinsic), and the pka. Its partners list every hydrogen bond, the fixed charges that move it by 0.05 or more and up to eight coupled groups that do, each with its distance and shift. A coupled group's shift is what its mean charge does to this group at this group's pKa: −wᵢⱼ⟨qⱼ⟩/ln10.

A worked example: hen lysozyme #

Hen egg-white lysozyme (PDB 1AKI, 129 residues, crystallised at pH 4.48) has 32 ionizable groups; its eight cysteines form four disulfides and do not titrate. Its two catalytic carboxylates show the model at work:

Glu35 sits in the active-site cleft, half buried and with a single hydrogen bond, and keeps its proton up to pH 5.3: in the enzyme it is the general acid that must stay protonated. Asp66, at the surface but held by four hydrogen bonds, gives up its proton already near pH 2.6, a unit below its model value. Glu7's pKa falls by a full unit through coupling, more than half of it from Lys1, whose ammonium sits 3.7 Å away (the partner list gives −0.58). Tyr53 rises by 0.62 through coupling because Asp66's negative charge is 3.6 Å from its phenol. The protein is basic:

Its estimated isoelectric point is 11.0 (10.8 for the same chain unfolded). Moving from blood to the stomach, eleven groups switch charge (every carboxylate and His15), the net charge rises from +7.7 to +18.6 and the fold loses about 5.7 kcal/mol of electrostatic stability (see Stability over pH).

What the answer holds #

Net charge and isoelectric point #

The net charge in an environment is the sum of every group's mean charge plus the fixed charges of metal ions and nucleic acid phosphates (net_charge); protein_charge is the groups alone. The isoelectric point (summary.pi) is the pH at which the folded structure's net charge crosses zero; pi_unfolded is the same for the chains unfolded, every group at its model pKa with no interactions. The charge curves over pH 0–14 are returned for both states (curves.charge_folded, curves.charge_unfolded).

Net charge governs how a protein moves in an electric field, how it binds to charged surfaces (membranes, DNA, chromatography resins) and how soluble it is: proteins are least soluble near their isoelectric point.

Stability over pH #

Protons bind differently to the folded and the unfolded protein. If the folded structure holds fewer protons than its unfolded chain at some pH, lowering the pH favours unfolding; the thermodynamic link (Wyman linkage) is

dΔG_unfold / dpH = kT · ln10 · (Q_unfolded − Q_folded)

where Q is the net charge (each bound proton adds one positive charge). GL1X integrates this over the pH grid, with the unfolded chains titrating at their model pKa values and the folded structure as titrated above, and reports the electrostatic part of the stability relative to the most stable pH (curves.stability, at most 0; summary.stability_optimum). A value of −5.7 kcal/mol between blood and the stomach means the fold is 5.7 kcal/mol less stable at pH 1.5 for electrostatic reasons; the whole folding stability of many small proteins is of that order, so a loss of that size can unfold them. This is the driving force of acid denaturation, and in the stomach it is part of digestion: unfolded proteins are what pepsin cuts.

The approach follows Yang and Honig (1993). It gives the shape of the profile, not the absolute stability (which also depends on hydrophobic packing and other forces GL1X does not model), and it assumes the unfolded chains titrate as free peptides.

Groups that switch charge #

The environments are compared in order: the first is the reference, and each other one is compared with it (comparisons). A group switches when its mean charge changes by 0.5 or more between the two; each switch gives the group, its charge in both environments and its pKa. A comparison also gives the change of net charge and, in the uniform mode, of stability.

Switching groups are where a protein's response to its environment lives: the histidines that change charge between blood and an endosome, the carboxylates that pick up protons in the stomach, the groups of a transporter's outer face that answer the acid microclimate.

The structure's own pH #

A structure is a snapshot taken in the conditions of its experiment. GL1X reads them from the file:

The pH of the experiment (the crystal's, else the NMR sample's, else the microscopy buffer's) is offered as the environment experiment, and by default the answer compares blood with it. Lysozyme 1AKI was crystallised at pH 4.48: at blood's 7.4, His15, Glu35 and Asp52 carry a different charge than in the crystal. The coordinates still show the molecule at the crystal's pH; GL1X evaluates other environments on the same coordinates, and says so in the answer's notes. Modern crystal structures are often also measured frozen, near 100 K.

Membrane proteins: a pH on each side #

A membrane protein is exposed to two compartments at once, and they need not share a pH: the outer face of a gut-wall transporter meets the acid microclimate while its inner face meets the cytosol; a lysosomal membrane protein faces pH 4.7 inside the lysosome and 7.2 in the cytosol. With membrane=auto or on, GL1X places the membrane, assigns every group to a side, gives the outer side the environments' pH and the inner side a pH of its own (inside_ph, 7.2 by default), and adds the lipid's burial to the groups that face it.

Placing the membrane #

The hydrocarbon core of a lipid bilayer is a slab about 30 Å thick. The surface of a membrane protein that faces it is apolar; the surface outside it is polar and charged. GL1X finds the slab that best separates the two:

1. Each protein residue's solvent-exposed area comes from the Shrake–Rupley calculation of Structure analysis, and its class from its kind: +1 for Ala, Val, Leu, Ile, Met, Phe, Cys (and selenomethionine), +0.5 for Trp and Tyr (which line the membrane's interfaces), 0 for Gly and Pro, −0.25 for Ser and Thr, −0.5 for Asn, Gln and His, −1 for Asp, Glu, Lys and Arg.
2. For 2,000 directions spread evenly over a half sphere, for every position of the slab's centre along the direction (1 Å steps) and for every half-thickness from 10 to 18 Å, it sums over the residues whose alpha carbon lies in the slab the net apolar exposed area S = Σ class × area and the exposed area A.
3. It keeps the slab with the highest S / √(A + 2000 Ų). The score favours neither the thickest slab (S alone keeps growing as long as each added shell is apolar on balance) nor the thinnest (S/A favours a few apolar residues).
4. The slab is taken for a membrane when S is at least 1,200 Ų, the net apolar share S/A at least 0.25, and the protein reaches more than 3 Å beyond both faces. With membrane=auto a structure that fails these tests is analysed with one pH everywhere, and the answer says so; with membrane=on the best slab is used anyway.

The idea is that of TMDET (Tusnády et al. 2004), in a simpler score. The placement is reported in membrane: the normal, a point on the mid-plane, the half-thickness, the score and the apolar share.

The cytoplasmic side #

Which face is the cytoplasm comes from the positive-inside rule: the loops of membrane proteins on the cytoplasmic side are rich in lysines and arginines (von Heijne 1992). GL1X counts the Lys and Arg residues within 15 Å beyond each face; the face with more is the inside, and the membrane's normal points from it to the outside. When the two counts differ by less than 2 the orientation is not clear, the answer says so, and flip=true swaps the sides. Proteins of outer membranes (bacterial porins) and of organelles may not follow the rule; check the orientation against what is known of the protein.

Lipid burial #

A group inside the slab that faces the lipid loses water as surely as a group buried in protein. For each group within 3 Å of the slab's faces or inside them, GL1X casts 16 rays in the membrane's plane, each 30 Å long, and counts those that meet no atom of the structure other than the group's own residue (protein, nucleic acid or ligand; points are tested every 1.5 Å from 3 to 30 Å, and an atom within 2 Å of a point blocks the ray). A group on the lipid-facing surface sees out along about half the directions, one lining a water-filled pore along none; the lipid exposure L is twice the open fraction, capped at 1. The lipid's burial is full within 3 Å inside the faces and falls to none 3 Å outside them, and it combines with the protein's burial as two independent sources:

B_total = 1 − (1 − B_protein) × (1 − B_membrane(z) × L)

Fixed charges in the slab take the membrane's burial without the ray test.

Compartments and the outside scan #

Each group's height above the mid-plane decides its side: above the slab it is outside, below it inside, within it membrane (outer half) or membrane (inner half), and a group in the membrane takes the pH of the side it is nearer. In each environment the answer gives the charge of the outside, the inside and the membrane groups, and curves.outside_scan gives the charge of each part as the outer pH runs from 0 to 14 with the inner pH held fixed.

The group pKa values of a membrane protein are those of the structure with the lipid in place, titrated at one pH; the environments then apply the two pH values.

A worked example: an aspartate transporter #

GltPh is an archaeal homologue of the human excitatory amino acid transporters, the proteins that move glutamate across cell membranes coupled to sodium (Boudker et al. 2007). In the file 2NWL as the OPM database oriented it (Lomize et al. 2006), whose header gives a half-thickness of 13.4 Å, GL1X places the membrane with its normal within 1.3° of OPM's, its mid-plane 0.6 Å from OPM's and a half-thickness of 13 Å, and finds the cytoplasmic side by 26 lysines and arginines within 15 Å of its face against 5 on the other. (OPM marks its membrane planes with thousands of dummy atoms, DUM; GL1X leaves them out.)

The file holds 209 ionizable groups: 204 of the trimer's three chains and 5 of the three bound aspartate substrates, which the file writes as atom records of a chain D and GL1X therefore reads as a short chain. 53 are outside, 85 inside and 71 in the membrane. With the inside at the cytosol's 7.2, the net charge and the charge of the groups of each part are:

The inner face hardly changes, as it should not: its groups feel the outer pH only through their couplings. The outer face goes from slightly negative in blood to clearly positive in an acid microclimate, and the membrane's outer half follows: Asp48 of chain A, 5.9 Å above the mid-plane with a burial of 0.55, has an estimated pKa of 5.4 and goes from −0.97 at pH 7.4 to −0.77 at 6.1; His132, at the outer interface (pKa 5.6), goes from +0.03 to +0.28, and to +1 at pH 1.5. Groups like these are the candidates for a transporter's pH sensing or proton coupling; the analysis points to them, and experiments decide.

How the model was set and checked #

Calibration against PROPKA 3 #

The burial bounds, the desolvation sizes, the hydrogen-bond weights and the protein dielectric were set by comparing GL1X with PROPKA 3.5.1 (Li et al. 2005; Olsson et al. 2011), an established empirical pKa program, run locally as a reference and not shipped with GL1X. Measured pKa values were out of reach when the model was set, so it was calibrated against a program; how it fares against measured values was checked afterwards (below).

76 protein structures from the test sets of open-source structural biology packages were analysed by both programs, 5,269 groups in all. They were split by a hash of their names into 36 for calibration and 40 held out. Because GL1X and PROPKA use different model pKa values (PROPKA uses, for example, 3.80 for Asp and 12.5 for Arg), the comparison is of shifts: each program's pKa minus its own model value. Arginines were left out of the scores, since GL1X's model value of 13.8 puts almost every arginine above pH 14. Parameters were set by coordinate descent on the calibration structures, minimising the root-mean-square difference of the shifts; the weight of hydrogen bonds to bases came out near zero and was kept at a small 0.2 rather than dropped, because a hydrogen bond does stabilise an ammonium.

Results on held-out structures #

On the 40 held-out structures (2,505 groups other than arginines) GL1X's shifts differ from PROPKA's by 0.69 pK units root-mean-square, and by 0.45 on average in absolute value; 71% of groups agree within 0.5 units, 92% within 1 and 98% within 2. Taking no shift at all (every group at its model pKa) would differ by 1.53 units root-mean-square.

GltPh's file in the held-out set comes from OPM, whose dummy membrane atoms a program would read as part of the structure: GL1X leaves them out, and PROPKA was run on the file without them.

The largest disagreements are for buried tyrosines and free cysteines, for which PROPKA predicts larger desolvation shifts (pKa values of 12–13) than GL1X does.

Membrane placement check #

The membrane placement was checked on GltPh against the OPM database (above): normal within 1.3°, mid-plane within 0.6 Å, half-thickness 13 Å against 13.4 Å. Among the 92 readable structure files of the benchmark folder (the 93rd holds no polymer chain), the slab passed the thresholds for the two membrane proteins (GltPh and the porin OmpK36, 1OSM) and for no soluble protein; the soluble structure that came closest (chain A of 1LPB) had a membrane-like score but was rejected because it does not reach beyond both faces. A unit test checks the placement and the orientation on a model with an apolar band between polar ends.

What the comparison does and does not show #

Agreement with PROPKA shows that GL1X captures the same structural effects with a different model: it is not a measure of accuracy against experiment. Against measured values, Olsson et al. (2011) reported root-mean-square errors for PROPKA 3 of 0.79 for Asp and Glu, 0.75 for Tyr, 0.65 for Lys and 1.00 for His, and errors of different programs need not add or cancel. GL1X's own errors against measured values are below.

Against measured pKa values #

Since 3.4.0 GL1X is also checked against experiment. Data: the wild-type table of PKAD, a curated database of measured protein pKa values (Pahari et al. 2019), read from its copy in the MCCE benchmarking repository; every value was cross-checked against PKAD-3 (the values of the KaML split; 331 of the 368 primary values appear there, identical but for two the original abstracts correct) and, where it could be read, PKAD-R, and each was checked against the abstract of its paper. One value per group was kept, from the reference giving the most values for its protein; values that are not a single well-defined pKa (multi-step or biphasic titrations, approximate or estimated values, values the sources disagree on) and values that do not belong to the deposited structure (a variant, another construct, protonation on denaturation) were left out — 87 in all, each with its reason. That leaves 368 measured values in 47 proteins, one PDB entry each, analysed with ph_environment at the ionic strength of the measurement (0.15 M when the source states none usable). Three of the groups coordinate a metal in the crystal (barnase Glu60 a zinc, ribonuclease T1 Asp15 a calcium, the N-terminus of thioredoxin a copper); GL1X holds such groups fixed and gives them no pKa, so 365 measured values are compared. The null model leaves every group at GL1X's model pKa (model_pka): a predictor of shifts must beat it to add anything.

Overall GL1X is no more accurate than the null model: 0.88 against 0.91 pK units root-mean-square, mean absolute errors 0.63 and 0.62, 81% of both within one unit; the difference, −0.03, has a 95% interval of −0.19 to +0.15 (2,000 bootstrap resamples of the proteins). GL1X is closer to the measurement for 179 of the 365 groups. Where it helps and where it hurts:

• Large shifts and their direction. GL1X earns its keep where groups are strongly shifted (1.82 against 2.64) and tracks the direction of shifts of 0.5 units or more in 120 of 168 groups; its shifts correlate with the measured ones (r 0.48), but a GL1X shift of one unit goes with a measured shift of about half a unit on average. Aspartates gain most (0.75 against 1.12; r 0.77).
• Histidines and buried groups. GL1X sets histidines too low on average (−0.60), where measured histidines sit slightly above their model value; its burial term pushes buried groups towards their neutral form, while in active sites the protein often stabilizes the charged form: 8 of the 10 largest errors are such groups (a buried cysteine of methionine sulfoxide reductase computed 5.3 units too high, a histidine of a phospholipase C 3.5 too low).
• Unshifted groups. Where the measured value is within half a unit of the model value, GL1X's shifts are noise (0.59 against 0.26).
• NMR models. On the first model of an NMR ensemble GL1X does worse than no shift at all (0.95 against 0.54).
• Salt. For 27 groups measured twice at different ionic strengths the measured values moved by 0.42 units root-mean-square, GL1X's by 0.77, with no correlation (r −0.22): its response to salt looks too strong. Two independent measurements of the same group differ by 0.35 — the floor any predictor faces on these data.
• Hen lysozyme. Glu35, measured at 6.2, comes out at 5.30–5.46 on three crystal structures (right direction, not far enough); Asp52, measured at 3.68, at 4.27–4.90 (wrong direction).

Limits of the benchmark: GL1X computes at 25 °C while the measurements range from 10 to 42 °C (errors inside and outside 20–30 °C are alike); a salt concentration stands for the ionic strength; some structures are not exactly the molecule measured (another species' ovomucoid, full-length L9 for its N-terminal domain, calcium-loaded calbindin for the apo protein — on the apo NMR structure calbindin's eight groups come out better, 0.44 against 0.56); the values are mostly NMR titrations of small soluble proteins, dominated by surface carboxylates with small shifts, which favours the null model, and buried groups, cysteines and lysines are few, so their numbers are uncertain. The values are PKAD's transcriptions, checked against the abstracts, not the full texts. The measured values with their sources and the reason for each one left out, GL1X's answers, the scripts that compute the tables and the full report come with the release (gl1x-3.4.0-pka-benchmark.tar.gz), and the benchmark is written up as a paper: GL1X pKa values against experiment.

Reading GL1X's values. Take a single pKa as an estimate within about a unit — more for histidines, cysteines, buried groups and NMR models — and look at its partners to see why it is shifted; the strongly shifted groups and the groups that switch charge between environments are where the analysis is most informative, as candidates to check, not facts. The model's parameters are as calibrated against PROPKA; recalibrating them against measured values (with a held-out split of proteins) is future work.

In the explorer #

The pH environments tool of the research workbench shows the record's structure on the left and the analysis beside it.

• Environments. Pick up to six, grouped as digestive tract, blood and tissue, and cell compartments; the experiment's pH appears once the file has been read, and any pH can be added. The first environment chosen is the reference (marked ref); click a chosen one to remove it.
• Settings. Ionic strength (0.15 M by default), the membrane mode (off, auto, on), the inside pH and Swap sides in membrane mode; Analyse runs the analysis (it also runs when the tool opens).
• Summary. The conditions of the experiment, the isoelectric point folded and unfolded, the most stable pH, and a tile per environment with its net charge (and, in membrane mode, the charge of each side).
• Comparisons. One line per environment against the reference: the change of net charge and stability and the switching groups, each a button that shows the group in 3D.
• Charts. Net charge over pH for the folded and unfolded protein, the stability profile, and in membrane mode the charge of each side as the outer pH varies; each chart marks the environments, and hovering reads the values.
• Colour. The structure is coloured by each residue's charge in an environment (blue positive, red negative, as electrostatic surfaces are drawn), shown as a surface, or by the change from the reference to an environment (blue for a residue that gained protons, red for one that lost them), shown as a cartoon; switching groups, when there are no more than 60, are drawn as sticks.
• Table. Every group with its pKa (the model value in small type), its shifts as four small bars (desolvation, hydrogen bonds, fixed charges, coupling; 12 pixels per pK unit, lighter when they lower the pKa, with a key above the table), its burial and its charge in each environment; filters keep the switching groups, the groups whose pKa moved by 1 or more, the buried ones, or all. A click on a group shows it in the structure with its partners and their shifts.

API and MCP #

curl -s 'https://YOUR-DOMAIN/api/bio/ph/environments'
curl -s 'https://YOUR-DOMAIN/api/bio/records/0xde3723766bc32dcacd03c17baa400a7b36837eba/1/ph?env=blood,stomach,experiment'
curl -s 'https://YOUR-DOMAIN/api/bio/records/0xde3723766bc32dcacd03c17baa400a7b36837eba/1/ph?env=7.4,6.1,4.7&ionic_strength=0.1&curves=1'
curl -s 'https://YOUR-DOMAIN/api/bio/records/{contract}/{tokenId}/ph?env=blood,brush-border&membrane=auto&inside_ph=7.2'

GET /api/bio/ph/environments lists the named environments (id, name, group, pH, spread, note, source), the defaults (inside pH 7.2, ionic strength 0.15 M) and the model's rules.

GET /api/bio/records/{contract}/{tokenId}/ph takes:

The answer holds conditions (the experiment's pH and what it was, temperatures, resolution, crystallisation method and details), model (version, temperature, ionic strength), environments (id, name, pH, net_charge, protein_charge, and stability, or in membrane mode inside_ph, outside_charge, inside_charge and membrane_charge; with the source of a named environment), sites (for each group: index, kind, label such as A:GLU35, chain, number, ins, residue, model_pka, pka or pka_below/pka_above, hill, intrinsic_pka, shift with desolvation, hydrogen_bonds, fixed_charges and coupling, burial, burial_atoms, partners, fixed and coordinated for metal-bound groups, irregular, compartment, height and lipid_facing in membrane mode, and charge and protonated with one value per environment), comparisons (from, to, charge_change, stability_change, switches, switch_count), summary (groups by kind, groups shifted by 1 or more, the groups not titrated by reason, metal ions, phosphates, pi, pi_unfolded, stability_optimum, the number of clusters and the largest), curves (ph, charge_folded, charge_unfolded, and stability or outside_scan), membrane and membrane_used when a membrane was placed, notes, and the model's rules.

Errors say why: bad_request (400: ionic_strength or inside_ph is not a number), invalid_contract, invalid_token, invalid_environment, too_many_environments, no_experiment_ph (the file states no pH), invalid_ionic_strength (outside 0 to 2), invalid_membrane, invalid_inside_ph (outside 0 to 14), record_not_found (404), structure_not_ready (409: the structure has not been read yet), structure_changed (409), structure_unavailable (502), structures_unavailable (503: the server cannot read structure files), structure_unreadable (422), no_protein (422), structure_too_large (422: more than 15,000 ionizable groups), ph_busy (503: two analyses run at once and a third waits up to 20 seconds), analysis_rate_limited (429: 20 at once and one more a second for each client). Analyses are kept for 30 minutes (the last 64); answers are cached for five minutes. On one core an analysis of lysozyme takes about 6 ms and one of GltPh with its membrane placed about 0.3 s.

MCP: ph_environment takes contract, token_id, environments, ionic_strength, membrane, inside_ph, flip, chain and max_groups. Its answer is the API's without the curves and the rules, with the groups that switch charge first and then those whose pKa moved most, at most max_groups of them (40 by default), each with up to four partners; it says how many groups it shows and how many there are. → MCP

Limits #

• Fixed coordinates. The structure is analysed as deposited: side chains do not move, loops do not open and the protein does not change shape with pH. Real proteins do all of these, and a group's pKa can shift when its surroundings reorganise; the strongly shifted groups are the most affected.
• One structure, first model. Crystal contacts with neighbouring molecules in the lattice are not considered, an NMR ensemble is represented by its first model, and alternate conformations by the first.
• Temperature and salt. The model works at 25 °C, where its model pKa values were measured, and with one ionic strength for every environment.
• What carries charge. Modified residues are not titrated (a phosphoserine's phosphate, which titrates near neutral pH, is missing), nor are nucleotide bases, ligands carry no charge except single-atom metal ions, and metals inside clusters and cofactors carry none. Waters in the file are ignored.
• Membranes. The membrane is a flat slab of fixed composition; curvature, lipid headgroup charges and the membrane potential (about −70 mV, inside negative, in many cells) are not modelled, and the positive-inside rule can orient a protein wrongly.
• What the file holds. GL1X analyses what the file holds. Ligands that a file writes as atom records (as some processed files do) are read as short chains and titrated; the dummy atoms with which OPM marks membrane planes (DUM) are left out; waters are ignored.
• Predicted models. A computed model's side chains are predictions too; the pKa values of its buried groups inherit their uncertainty.
• Accuracy. Against 365 measured values GL1X's error is 0.88 pK units root-mean-square, no better overall than leaving every group at its model value (0.91); it is better for aspartates and strongly shifted groups, worse for histidines, cysteines, barely shifted groups and NMR models (see above). Individual values can be off by more than a unit, by several for buried active-site groups.

What a physiological structure would need #

A crystal structure cannot be recalculated in other conditions: its electron density was measured in the crystal, and the coordinates describe that density. What can be asked is how the molecule would behave in a physiological environment, and that has levels:

1. Protonation and charge on the deposited coordinates, which is what this tool does: which groups change charge, how the net charge and the electrostatic stability change, and which groups would feel a change of environment most.
2. Conformational response, which needs simulation: molecular dynamics of the protein in explicit water or a membrane at 310 K and 150 mM salt, with protonation states set for the target pH (or allowed to change during the simulation, as constant-pH molecular dynamics does). This takes hours to days of GPU time per structure and is outside what a single explorer server does; GL1X's answer gives such a simulation its starting protonation states and tells it which groups to watch.

References #

• Boudker O, Ryan RM, Yernool D, Shimamoto K, Gouaux E (2007) Coupling substrate and ion binding to extracellular gate of a sodium-dependent aspartate transporter. Nature 445:387–393.
• Casey JR, Grinstein S, Orlowski J (2010) Sensors and regulators of intracellular pH. Nat Rev Mol Cell Biol 11:50–61.
• Debye P, Hückel E (1923) Zur Theorie der Elektrolyte. I. Gefrierpunktserniedrigung und verwandte Erscheinungen. Phys Z 24:185–206.
• Evans DF, Pye G, Bramley R, Clark AG, Dyson TJ, Hardcastle JD (1988) Measurement of gastrointestinal pH profiles in normal ambulant human subjects. Gut 29:1035–1041.
• Fitch CA, Platzer G, Okon M, García-Moreno E B, McIntosh LP (2015) Arginine: its pKa value revisited. Protein Sci 24:752–761.
• Gilson MK (1993) Multiple-site titration and molecular modeling: two rapid methods for computing energies and forces for ionizable groups in proteins. Proteins 15:266–282.
• Grimsley GR, Scholtz JM, Pace CN (2009) A summary of the measured pK values of the ionizable groups in folded proteins. Protein Sci 18:247–251.
• Hamm LL, Nakhoul N, Hering-Smith KS (2015) Acid-base homeostasis. Clin J Am Soc Nephrol 10:2232–2242.
• Kabsch W, Sander C (1983) Dictionary of protein secondary structure: pattern recognition of hydrogen-bonded and geometrical features. Biopolymers 22:2577–2637.
• Li H, Robertson AD, Jensen JH (2005) Very fast empirical prediction and rationalization of protein pKa values. Proteins 61:704–721.
• Lomize MA, Lomize AL, Pogozheva ID, Mosberg HI (2006) OPM: orientations of proteins in membranes database. Bioinformatics 22:623–625.
• Lucas M (1983) Determination of acid surface pH in vivo in rat proximal jejunum. Gut 24:734–739.
• Olsson MHM, Søndergaard CR, Rostkowski M, Jensen JH (2011) PROPKA3: consistent treatment of internal and surface residues in empirical pKa predictions. J Chem Theory Comput 7:525–537.
• Pahari S, Sun L, Alexov E (2019) PKAD: a database of experimentally measured pKa values of ionizable groups in proteins. Database 2019:baz024. doi:10.1093/database/baz024
• Shrake A, Rupley JA (1973) Environment and exposure to solvent of protein atoms. Lysozyme and insulin. J Mol Biol 79:351–371.
• Thurlkill RL, Grimsley GR, Scholtz JM, Pace CN (2006) pK values of the ionizable groups of proteins. Protein Sci 15:1214–1218.
• Thwaites DT, Anderson CMH (2007) H⁺-coupled nutrient, micronutrient and drug transporters in the mammalian small intestine. Exp Physiol 92:603–619.
• Tusnády GE, Dosztányi Z, Simon I (2004) Transmembrane proteins in the Protein Data Bank: identification and classification. Bioinformatics 20:2964–2972.
• von Heijne G (1992) Membrane protein structure prediction. Hydrophobicity analysis and the positive-inside rule. J Mol Biol 225:487–494.
• Yang AS, Honig B (1993) On the pH dependence of protein stability. J Mol Biol 231:459–474.