Zinc-Methylimidazole Speciation & Kinetics Model

A computational chemistry tool that predicts pH-dependent speciation of Zn(II) with N-methylimidazole and derives second-order rate constants for CO2 hydration from the resulting active-species distribution. Built to model the catalytic behavior of carbonic anhydrase biomimetic complexes.

Scientific Context

Carbonic anhydrase mimics based on Zn-imidazole coordination reproduce the enzyme's inner-sphere hydroxide mechanism for CO2 hydration. The catalytic rate depends on which zinc-imidazole species are present in solution, and these populations shift dramatically with the ligand-to-metal ratio.

This model computes the full speciation across 10 zinc species using a partition-function approach, then predicts observable rate constants by weighting each catalytically active species by its quantum-chemically derived intrinsic rate.

Model Architecture

Speciation (Partition Function)

The model tracks 10 species in the Zn-methylimidazole system:

Index	Species	Description
0	M	Free Zn(II)
1-6	MA1 - MA6	Successive ligand complexes
7	MH⁻¹A3	Tris-hydroxo (catalytic)
8	MH⁻¹A4	Tetrakis-hydroxo (catalytic)
9	MH⁻²A4	Tetrakis-dihydroxo

Species fractions are computed from the binding polynomial:

Q = 1 + Σ β(1,0,n)[L]^n + β(MH⁻¹A3)[L]³[H⁺]⁻¹ + β(MH⁻¹A4)[L]⁴[H⁺]⁻¹ + β(MH⁻²A4)[L]⁴[H⁺]⁻²

Free ligand concentration [L] is solved via bisection on the ligand mass balance.

Kinetics Pipeline

The predicted rate constant is built in progressive stages:

1. Speciation-only -- weighted sum of intrinsic rates over active hydroxo species
2. Bicarbonate inhibition -- product inhibition from buffer-derived HCO3⁻
3. Logistic damping -- tetrakis channel attenuation at high ligand excess
4. Mass-action damping -- alternative competitive-displacement model

Thermodynamic Constants

All equilibrium constants from Appleton & Sarkar (1974), measured at 0.16 M KNO3, 25 °C:

• Stepwise log K: 2.380, 2.544, 1.676, 2.614, 0.791, 1.042
• pKa(N-MeIm): 7.209
• pKa(M-OH2): 9.12
• log β(MH⁻¹A4): 0.157
• log β(MH⁻²A4): -10.615

Project Structure

.
├── speciation/                    # Core library
│   ├── config.py                  # Model configuration and constants
│   ├── thermodynamics.py          # Partition function and species fractions
│   ├── solver.py                  # Free-ligand mass-balance solver
│   ├── kinetics.py                # Rate predictions, inhibition, fitting
│   ├── plotting.py                # Figure generation
│   └── runner.py                  # Pipeline orchestrator
├── tests/                         # Automated test suite (pytest)
│   ├── test_thermodynamics.py     # Species fraction properties
│   ├── test_solver.py             # Mass-balance convergence
│   └── test_kinetics.py           # Rate predictions + regression checks
├── run_model.py                   # Entry point: default configuration
├── run_corrected.py               # Entry point: ionic-strength corrected
├── make_publication_figure.py     # Publication-quality two-panel figure
├── data/
│   └── observed.xlsx              # Experimental rate constants
├── docs/
│   ├── model_primer.md            # Pedagogical guide to the model
│   ├── supplementary_derivation.md  # Full mathematical derivation
│   └── breakdown.md               # Non-technical project walkthrough
├── .github/workflows/tests.yml    # CI: runs pytest on push/PR
├── requirements.txt
└── LICENSE

Installation

python -m venv .venv
source .venv/bin/activate        # Windows: .venv\Scripts\Activate
pip install -r requirements.txt

Requires Python 3.9+.

Usage

Run the speciation model

python run_model.py

Outputs written to outputs/:

• full_speciation_all_compositions.csv -- complete speciation matrix across all compositions
• prediction_summary_by_ratio.csv -- averaged predictions by ligand:Zn ratio
• phi_correction.csv -- multiplicative correction factors for high-throughput screening
• overlay_stage*.png -- progressive overlay plots comparing predicted vs observed rates

Run the ionic-strength-corrected variant

python run_corrected.py

Applies Davies equation corrections from I = 0.16 M (literature) to I = 0.125 M (experimental). Outputs written to outputs_corrected/.

Generate publication figure

python make_publication_figure.py

Produces a two-panel figure with observed vs predicted rate constants and a stacked-area speciation breakdown.

Run tests

python -m pytest tests/ -v

The test suite covers thermodynamic identities (species fractions sum to 1), solver convergence, physical constraints (non-negative rates), and regression checks against known-good outputs.

Configuration

All model parameters are defined in speciation/config.py under DEFAULT_CONFIG. Key settings:

DEFAULT_CONFIG = {
    "pH": 9.00,
    "ratios_L_to_Zn": [3, 4, 10, 15, 25, 37.5, 50, 75, 100, 175, 250],
    "Zn_totals_mM": [1.00, 0.50, 0.25],
    "k_MA3_OH": 511.0,       # intrinsic rate for tris-hydroxo (M⁻¹ s⁻¹)
    "k_MA4_OH": 1494.0,      # intrinsic rate for tetrakis-hydroxo (M⁻¹ s⁻¹)
    ...
}

The corrected variant (run_corrected.py) overrides only three values from the default, demonstrating how the shared codebase eliminates duplication.

References

• Appleton, D. W. & Sarkar, B. (1974). The activity-related ionization in carbonic anhydrase model systems. J. Biol. Chem., 249(14), 4780-4786.
• Rains, G. et al. (2019). Bicarbonate inhibition of carbonic anhydrase mimics hinders catalytic efficiency. Dalton Trans., 48, 5045-5056.

License

MIT License. See LICENSE for details.