Skip to content

API Reference

deuterium_uptake(pf, k_int, time)

Compute time-dependent deuterium uptake using EX2 kinetics. D(t) = 1 - exp(-k_obs * t) where k_obs = k_int / PF (Hvidt & Nielsen, 1966).

Parameters:

Name Type Description Default
pf ndarray

(N,) protection factors.

 required
k_int ndarray

(N,) intrinsic exchange rates.

 required
time float

Exposure time in minutes.

 required

Returns:

Type Description
ndarray

D(t) (N,) fractional deuterium uptake.
Source code in diff_hdx/kernels.py 
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
def deuterium_uptake(
    pf: jnp.ndarray,
    k_int: jnp.ndarray,
    time: float,
) -> jnp.ndarray:
    """
    Compute time-dependent deuterium uptake using EX2 kinetics.
    D(t) = 1 - exp(-k_obs * t)
    where k_obs = k_int / PF (Hvidt & Nielsen, 1966).

    Args:
        pf: (N,) protection factors.
        k_int: (N,) intrinsic exchange rates.
        time: Exposure time in minutes.

    Returns:
        D(t) (N,) fractional deuterium uptake.
    """
    k_obs = k_int / pf
    return 1.0 - jnp.exp(-k_obs * time)

h_bond_energy(donor_coords, acceptor_coords, cutoff=3.5, sigma=0.5)

Compute a differentiable approximation of H-bond energy/count. Uses a sigmoid-like distance cutoff.

Parameters:

Name Type Description Default
donor_coords ndarray

(N, 3) coordinates of donors.

 required
acceptor_coords ndarray

(M, 3) coordinates of acceptors.

 required
cutoff float

Distance cutoff in Angstroms.

 3.5
sigma float

Smoothing parameter for the transition.

 0.5

Returns:

Type Description
ndarray

Approximate H-bond energy/count for each donor (N,).
Source code in diff_hdx/kernels.py 
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
def h_bond_energy(
    donor_coords: jnp.ndarray,
    acceptor_coords: jnp.ndarray,
    cutoff: float = 3.5,
    sigma: float = 0.5,
) -> jnp.ndarray:
    """
    Compute a differentiable approximation of H-bond energy/count.
    Uses a sigmoid-like distance cutoff.

    Args:
        donor_coords: (N, 3) coordinates of donors.
        acceptor_coords: (M, 3) coordinates of acceptors.
        cutoff: Distance cutoff in Angstroms.
        sigma: Smoothing parameter for the transition.

    Returns:
        Approximate H-bond energy/count for each donor (N,).
    """
    # Compute pairwise distances (N, M)
    diff = donor_coords[:, None, :] - acceptor_coords[None, :, :]
    dist_sq = jnp.sum(diff**2, axis=-1)
    # Safe distance for gradients
    dist = jnp.sqrt(jnp.where(dist_sq > 0, dist_sq, 1.0))
    dist = jnp.where(dist_sq > 0, dist, 0.0)

    # Soft-cutoff: 1 / (1 + exp((r - r_cutoff) / sigma))
    # Sum over all potential acceptors for each donor
    hb_counts = jnp.sum(jax.nn.sigmoid((cutoff - dist) / sigma), axis=-1)
    return hb_counts

intrinsic_rates(sequence, ph=7.0, temperature=293.15)

Compute intrinsic exchange rates (k_int) using the Bai et al. (1993) model. Includes full side-chain correction factors for all 20 standard amino acids.

Per Bai et al. (1993) the correction for residue i uses: - the left neighbour (residue i-1) via the "al" / "bl" factors, and - the right neighbour (residue i+1) via the "ar" / "br" factors. Boundary residues (N-terminus, C-terminus) use Ala as a placeholder.

This implementation is fully vectorised and compatible with JAX JIT.

Parameters:

Name Type Description Default
sequence str

Protein sequence string (one-letter amino-acid codes).

 required
ph float

pH value.

 7.0
temperature float

Temperature in Kelvin.

 293.15

Returns:

Type Description
Array

k_int array of shape (N,), rates in min⁻¹.
Source code in diff_hdx/kernels.py 
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
def intrinsic_rates(
    sequence: str,
    ph: float = 7.0,
    temperature: float = 293.15,
) -> Array:
    """
    Compute intrinsic exchange rates (k_int) using the Bai et al. (1993) model.
    Includes full side-chain correction factors for all 20 standard amino acids.

    Per Bai et al. (1993) the correction for residue *i* uses:
      - the **left** neighbour  (residue i-1) via the "al" / "bl" factors, and
      - the **right** neighbour (residue i+1) via the "ar" / "br" factors.
    Boundary residues (N-terminus, C-terminus) use Ala as a placeholder.

    This implementation is fully vectorised and compatible with JAX JIT.

    Args:
        sequence: Protein sequence string (one-letter amino-acid codes).
        ph: pH value.
        temperature: Temperature in Kelvin.

    Returns:
        k_int array of shape (N,), rates in min⁻¹.
    """
    n = len(sequence)  # noqa: F841 -- kept for readability; not used in vectorised ops

    if not sequence:
        return jnp.zeros(0)

    # Encode sequence as integer indices (unknown residues → Ala)
    seq_idx = [_AA_IDX.get(aa, _ALA_IDX) for aa in sequence]

    # Left-neighbour indices: residue i-1; N-terminal boundary → Ala
    left_idx = [_ALA_IDX] + seq_idx[:-1]
    # Right-neighbour indices: residue i+1; C-terminal boundary → Ala
    right_idx = seq_idx[1:] + [_ALA_IDX]

    # Look up correction arrays — pure Python lists, converted to JAX once
    corr = jnp.array(_CORRECTIONS)  # (20, 4)
    left_corr = corr[jnp.array(left_idx)]  # (N, 4)
    right_corr = corr[jnp.array(right_idx)]  # (N, 4)

    # Reference rates for NH in H₂O at 20 °C (293.15 K)
    k_a_ref = 10.0**1.39
    k_b_ref = 10.0**10.08
    k_w_ref = 10.0**-1.50  # estimated

    # [H⁺] and [OH⁻]; pKw at 20 °C ≈ 14.17
    h_plus = 10.0 ** (-ph)
    oh_minus = 10.0 ** (ph - 14.17)

    # Arrhenius temperature corrections (activation energies in kcal/mol)
    e_a, e_b, e_w = 14.0, 17.0, 19.0
    r_gas = 1.987e-3  # kcal / (mol·K)

    def temp_corr(k_ref: float, e_act: float) -> jnp.ndarray:
        return k_ref * jnp.exp(-e_act / r_gas * (1.0 / temperature - 1.0 / 293.15))  # type: ignore[no-any-return]

    ka_ref_t = temp_corr(k_a_ref, e_a)
    kb_ref_t = temp_corr(k_b_ref, e_b)
    kw_ref_t = temp_corr(k_w_ref, e_w)

    # Log-additive corrections — vectorised over all residues simultaneously
    # Columns: [al=0, ar=1, bl=2, br=3]
    ka = ka_ref_t * 10.0 ** (left_corr[:, 0] + right_corr[:, 1])  # al + ar
    kb = kb_ref_t * 10.0 ** (left_corr[:, 2] + right_corr[:, 3])  # bl + br
    kw = kw_ref_t * 10.0 ** (left_corr[:, 2] + right_corr[:, 3])  # same as kb

    rates = ka * h_plus + kb * oh_minus + kw

    # Proline has no amide proton, so its intrinsic rate is physically 0.0.
    # Note: we use 0.0 here as the rate; downstream tools should handle this
    # appropriately if a sentinel value or NaN is preferred instead.
    pro_idx = _AA_IDX.get("P", -1)
    pro_mask = jnp.where(jnp.array(seq_idx) == pro_idx, 0.0, 1.0)

    return jnp.asarray(rates * pro_mask)  # explicit Array, satisfies mypy

protection_factors(coords, h_bond_energies, beta_c=1.0, beta_asa=1.0, probe_radius=1.4)

Compute HDX protection factors (PF). PF = k_int / k_obs

Uses the Linderstrøm-Lang model with separate scaling coefficients for H-bond and burial contributions:

ln(PF) = beta_c * N_HB + beta_asa * (1 − SASA)

Both coefficients default to 1.0, matching the original single-beta formulation for backward compatibility. When fitting against experimental protection factors, beta_c and beta_asa should be treated as independent free parameters.

Parameters:

Name Type Description Default
coords ndarray

(N, 3) coordinates.

 required
h_bond_energies ndarray

(N,) hydrogen bond energies (or counts).

 required
beta_c float

Scaling coefficient for the H-bond contribution.

 1.0
beta_asa float

Scaling coefficient for the burial (1 − SASA) contribution.

 1.0
probe_radius float

Solvent probe radius passed to sasa_approx (Å).

 1.4

Returns:

Type Description
ndarray

PF (N,) protection factors.
Source code in diff_hdx/kernels.py 
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
def protection_factors(
    coords: jnp.ndarray,
    h_bond_energies: jnp.ndarray,
    beta_c: float = 1.0,
    beta_asa: float = 1.0,
    probe_radius: float = 1.4,
) -> jnp.ndarray:
    """
    Compute HDX protection factors (PF).
    PF = k_int / k_obs

    Uses the Linderstrøm-Lang model with separate scaling coefficients for
    H-bond and burial contributions:

        ln(PF) = beta_c * N_HB + beta_asa * (1 − SASA)

    Both coefficients default to 1.0, matching the original single-beta
    formulation for backward compatibility.  When fitting against experimental
    protection factors, beta_c and beta_asa should be treated as independent
    free parameters.

    Args:
        coords: (N, 3) coordinates.
        h_bond_energies: (N,) hydrogen bond energies (or counts).
        beta_c: Scaling coefficient for the H-bond contribution.
        beta_asa: Scaling coefficient for the burial (1 − SASA) contribution.
        probe_radius: Solvent probe radius passed to sasa_approx (Å).

    Returns:
        PF (N,) protection factors.
    """
    sasa = sasa_approx(coords, probe_radius=probe_radius)
    # ln PF = beta_asa*(1 − SASA) + beta_c*N_HB
    ln_pf = beta_asa * (1.0 - sasa) + beta_c * h_bond_energies
    return jnp.exp(ln_pf)

sasa_approx(coords, probe_radius=1.4, sigma=2.0)

Differentiable approximation of Solvent Accessible Surface Area (SASA). Uses a Gaussian occlusion model.

The probe radius is incorporated as an additive contribution to the effective Gaussian width (effective_sigma = sigma + probe_radius), so a larger probe widens the occlusion shell around each atom, reducing the accessible surface — consistent with standard SASA intuition.

Note: this is a differentiable surrogate, not a true Shrake–Rupley SASA. It lacks per-atom van-der-Waals radii and returns dimensionless values in (0, 1]. It is suitable as a smooth proxy for gradient-based refinement.

Parameters:

Name Type Description Default
coords ndarray

(N, 3) atomic coordinates in Angstroms.

 required
probe_radius float

Radius of the solvent probe in Angstroms (default 1.4 Å).

 1.4
sigma float

Base Gaussian width for the occlusion kernel in Angstroms.

 2.0

Returns:

Type Description
ndarray

Approximate accessibility values (N,) in (0, 1]; 1 = fully exposed.
Source code in diff_hdx/kernels.py 
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
def sasa_approx(
    coords: jnp.ndarray,
    probe_radius: float = 1.4,
    sigma: float = 2.0,
) -> jnp.ndarray:
    """
    Differentiable approximation of Solvent Accessible Surface Area (SASA).
    Uses a Gaussian occlusion model.

    The probe radius is incorporated as an additive contribution to the
    effective Gaussian width (effective_sigma = sigma + probe_radius), so
    a larger probe widens the occlusion shell around each atom, reducing the
    accessible surface — consistent with standard SASA intuition.

    Note: this is a differentiable *surrogate*, not a true Shrake–Rupley SASA.
    It lacks per-atom van-der-Waals radii and returns dimensionless values in
    (0, 1].  It is suitable as a smooth proxy for gradient-based refinement.

    Args:
        coords: (N, 3) atomic coordinates in Angstroms.
        probe_radius: Radius of the solvent probe in Angstroms (default 1.4 Å).
        sigma: Base Gaussian width for the occlusion kernel in Angstroms.

    Returns:
        Approximate accessibility values (N,) in (0, 1]; 1 = fully exposed.
    """
    # Effective width combines atom-atom smoothing and the probe size
    effective_sigma = sigma + probe_radius

    # Pairwise squared distances
    diff = coords[:, None, :] - coords[None, :, :]
    dist_sq = jnp.sum(diff**2, axis=-1)

    # Occlusion kernel: nearby atoms reduce accessibility.
    # Subtract the self-contribution (exp(0) = 1) from each row.
    occlusion = jnp.sum(jnp.exp(-dist_sq / (2 * effective_sigma**2)), axis=-1) - 1.0
    accessibility = 1.0 / (1.0 + occlusion)

    return accessibility