Model

Estimation of GZMK Concentration in Mucus

Overview: A Multi-Scale Model for a Clinical Question

Our DOCTOR project aims to develop a rapid diagnostic for Chronic Rhinosinusitis (CRS) by detecting Granzyme K (GZMK) in nasal mucus. A critical unknown, however, is the actual concentration of GZMK available for detection. To bridge this knowledge gap and establish a theoretical foundation for our device, we constructed a predictive mathematical model.

The core of our approach is a multi-scale framework that connects the large-scale dynamics of the sinus cavity with the microscopic properties of the cellular barrier. At the macroscopic level, we model the overall sinus environment, establishing a dynamic equilibrium where the diffusion of GZMK from inflamed tissue is balanced by its removal through natural mucus clearance. This equilibrium determines the final steady-state concentration of the biomarker.

To accurately quantify this process, the model zooms into the microscopic level, treating the nasal epithelium as a porous barrier. Here, we analyze how the physical properties of both the GZMK molecule and the intercellular pores hinder transport from tissue to mucus. Crucially, our model incorporates how CRS-induced inflammation damages this barrier, altering its permeability. By linking these two scales, our model provides a comprehensive and theoretically grounded estimate of the GZMK concentration, offering essential guidance for the sensitivity requirements of our diagnostic test.

Model Design and Assumptions

Our model integrates two distinct scales: a macroscopic view of the sinus environment and a microscopic view of the epithelial barrier. It is built on the following key assumptions:

1. Constant GZMK concentration in tissues: In the state of rhinosinusitis, the rates of GZMK production and consumption in the sinus mucosal tissue are relatively stable. Thus, the tissue concentration can be considered constant over a certain period.
2. Dynamic equilibrium of the system: Mucus is continuously renewed and cleared. After GZMK enters mucus from tissues, it is continuously carried away. When the transport rate of GZMK from tissues to mucus equals its clearance rate with mucus, the concentration in mucus reaches a steady state.
3. Passive, Paracellular Transport of GZMK: We assume that GZMK, as a large protein, primarily crosses the epithelial barrier through the water-filled channels of the cell-cell tight junctions (the paracellular pathway) via passive diffusion. Transport directly through the cells is considered negligible.
4. Multi-Pore Medium Model for Epithelium: The nasal epithelial layer is modeled as a porous barrier, where the paracellular pathways are treated as equivalent cylindrical pores. The transport of GZMK is therefore subject to steric and hydrodynamic hindrance within these pores.

Model Establishment

Macroscopic Equilibrium

According to Fick's first law, the diffusion rate of GZMK from tissues to mucus (J_diff) is proportional to the concentration difference across the mucosal barrier:

$$ J_{\text{diff}} = P \cdot A \cdot (C_t - C_m) $$

The clearance rate of GZMK (J_clear) is determined by the rate of mucus renewal:

$$ J_{\text{clear}} = Q \cdot C_m $$

At equilibrium:

$$ P \cdot A \cdot (C_t - C_m) = Q \cdot C_m $$

Solving for steady-state mucus concentration:

$$ C_m = \frac{P \cdot A}{P \cdot A + Q} \cdot C_t $$

Microscopic Permeability

The permeability coefficient (P) is determined by molecular transport through epithelial pores:

$$ P = \left( \frac{A_p}{A_{total}} \right) \cdot \frac{D_{pore}}{h_{epithelium}} $$

Diffusion inside pores is:

$$ D_{pore} = K(\lambda) \cdot D_w $$

Free diffusion coefficient:

$$ D_w = \frac{k_B T}{6 \pi \eta r_s} $$

Hindrance function:

$$ K(\lambda) = (1 - \lambda)^2 \cdot \left(1 - 2.104\lambda + 2.09\lambda^3 - 0.95\lambda^5\right) $$

Barrier damage coefficient from TEER:

$$ K_{damage} = \frac{\text{TEER}_{\text{normal}}}{\text{TEER}_{\text{CRS}}} $$

Pore size and area scaling under CRS:

$$ r_p = K_{damage}^{1/2} \cdot r_{p0} $$

$$ \frac{A_p}{A_{total}} = K_{damage} \cdot \left(\frac{A_p}{A_{total}}\right)_0 $$

Model Parameter

Parameter	Symbol	Value/Range	Justification / Reference
Polyp Surface Area	A	3–100 cm² (Median: 12.5 cm²)	Testa et al. (2020)
Mucus Clearance Rate	Q	2 mL/h	Gizurarson (2015)
GZMK Tissue Conc.	Ct	0.5–12.5 ng/mL (Avg: 2.5 ng/mL)	Lan et al. (2025)
GZMK Hydrodynamic Radius	rs	2.05 nm	Estimated from MW
Epithelial Thickness	hepithelium	50 μm	Typical airway epithelium
Baseline Pore Radius	rp0	10 nm	Deen (1987)
Baseline Pore Area Fraction	(Ap/Atotal)0	5 × 10⁻⁵	Renkin (1954)
Barrier Damage Coefficient	Kdamage	3–9	Ramezanpour et al. (2025)

Calculation and Sensitivity Analysis

Using median parameter values, our model predicts a baseline GZMK concentration in mucus of:

C_m = 2.64 × 10⁻³ ng/mL

The maximum concentration is:

C_m = 1.2084 ng/mL

A sensitivity analysis was conducted using a tornado plot:

Figure 1: Tornado plot illustrating sensitivity of predicted GZMK concentration in mucus. The black dashed line marks baseline 2.64e−03 ng/mL.

Analysis of Key Parameters

• High-Impact Parameters: A, C_t, and (A_p/A_total)₀. Larger polyp surface area can increase GZMK concentration above 0.02 ng/mL.
• Low-Impact Parameters: h_epithelium and r_p0.

Conclusion and Theoretical Feasibility

Our model predicts a baseline GZMK concentration of 2.64 pg/mL.

Typical colloidal gold LFAs detect in low ng/mL range, but optimized systems can reach pg/mL sensitivity. Thus, detection of GZMK in mucus is theoretically feasible.

Reference

Deen, W. M. (1987). “Hindered transport of large molecules in liquid-filled pores”. AIChE Journal.

Renkin, E. M. (1954). “Filtration, diffusion, and molecular sieving…”. J Gen Physiol.

Ramezanpour et al. (2025). Int Forum Allergy Rhinol.

Testa et al. (2020). Am J Case Rep.

Gizurarson S. (2015). Biol Pharm Bull.

Lan et al. (2025). Nature.

GZMK Active Structure Prediction and Refinement

Since accurate structural information is indispensable for inhibitor design and virtual screening, our first step was to reconstruct the active conformation of GZMK. The only crystal structure currently available is a catalytically inactive S195A pro-enzyme mutant, which was engineered to avoid autocleavage but does not represent the native functional state. To overcome this limitation, we employed a combination of AlphaFold-based prediction and molecular dynamics (MD) refinement to obtain a reliable structural model under near-physiological conditions.

We began by submitting the full-length GZMK amino acid sequence to the AlphaFold server, which generated five candidate models. The best-ranked structure achieved a predicted TM-score:

$$ \text{pTM} \approx 0.92 $$

indicating high confidence in its global fold.

Structural alignment against the pro-enzyme crystal structure confirmed that the predicted fold retained the expected serine protease backbone architecture. Alignment quality was quantified by the root-mean-square deviation (RMSD) between Cα atoms:

$$ \text{RMSD} = \sqrt{\frac{1}{N}\sum_{i=1}^N \left| \mathbf{r}_i^{\text{pred}} - \mathbf{r}_i^{\text{exp}} \right|^2} $$

where $N$ is the number of aligned atoms, $\mathbf{r}_i^{\text{pred}}$ are coordinates from the AlphaFold model, and $\mathbf{r}_i^{\text{exp}}$ are from the crystal structure. RMSD values below $\sim 2$ Å indicated strong topological consistency. This phase established a high-fidelity initial model.

A potential concern was whether post-translational modifications, particularly N-glycosylation, might alter the protein’s surface chemistry and confound hotspot selection. To evaluate this, we applied the Re-Glyco tool from the GlycoShape platform, which computationally screens for canonical N–X–S/T motifs with sufficient solvent accessibility. No reliable glycosylation sites were detected on the mature protein, suggesting that surface interaction pockets would not be significantly affected by glycan modifications. This outcome allowed us to focus our subsequent efforts on structural and energetic hotspot identification, free from the uncertainty of unknown PTMs.

While AlphaFold provides reliable backbone topology, it is less accurate for side-chain orientations and flexible loops, which are essential for ligand docking. To refine the model, we performed a 100 ns MD simulation using the Desmond module under near-physiological conditions (310 K, 0.15 M NaCl, S-OPLS force field). The time evolution of the structure was monitored by RMSD($t$):

$$ \text{RMSD}(t) = \sqrt{\frac{1}{N}\sum_{i=1}^N \left| \mathbf{r}_i(t) - \mathbf{r}_i(0) \right|^2} $$

where $\mathbf{r}_i(t)$ is the position of atom $i$ at time $t$, and $\mathbf{r}_i(0)$ is its reference position.

Throughout the trajectory, RMSD fluctuations remained below 3.8 Å, demonstrating the global stability of the protein structure and the absence of large-scale structural drift. Conformational clustering identified a dominant equilibrium state, and the centroid of this cluster was selected as the final optimized model.

Through this integrative approach—AlphaFold prediction, glycosylation site screening, and MD relaxation—we obtained a structurally and dynamically stable model of wild-type GZMK. This optimized structure provides not only a robust framework for defining druggable pockets but also a realistic basis for subsequent energetic analyses and inhibitor screening.