Code repository for MSc Physics thesis: "Coarse-Grained Molecular Dynamics of MUC2 Mucin Rheology, Pore Characteristics, and Chain Length Dependencies"
Author: Jo Lavoie Institution: Concordia University, Department of Physics Supervisor: Dr. Re Mansbach Date: 2024-2026
This repository contains simulation setup, analysis, and visualization code for studying the rheological properties of MUC2 mucin using coarse-grained molecular dynamics. Simulations were performed using LAMMPS on the Narval cluster (Compute Canada/Alliance).
The project investigates concentration scaling, chain-length dependence, and pore-size distributions of MUC2 mucin solutions using a bead-spring model based on Ford et al. (2021).
muc2_simulations/
├── README.md
├── .gitignore
├── rheology_common.py # Shared rheological analysis module (GK method)
│
├── parameter_sweep/ # Phase 1 parameter sweep (4,410 simulations)
│ ├── parameter_sweep_generator.py
│ ├── generate_slurm_arrays.py
│ ├── phase2_corrected_generator.py
│ ├── analyze_sweep_results.py
│ ├── analyze_phase2_corrected.py
│ ├── in.muc2_production
│ └── README.md
│
├── phase1_corrected/ # Phase 1 with CORRECTED full LJ interactions
│ ├── generate_phase1_narval.py
│ ├── analyze_phase1_corrected.py
│ ├── analyze_phase1_thesis.py
│ ├── compare_wca_vs_lj.py
│ ├── create_phase1_scripts.py
│ └── INSTRUCTIONS.md
│
├── phase2_analysis/ # Phase 2 P1Winners optimized parameter study
│ ├── generate_phase2_p1winners.py
│ ├── analyze_phase2_rheology.py
│ ├── analyze_phase2_pores.py
│ ├── analyze_winners.py
│ ├── create_phase2_figures.py
│ └── create_pore_thesis_figure.py
│
├── chain_length_study/ # Chain length dependence studies
│ ├── generate_chain_length_sims.py # Pre-correction generator (eps_TT=0.4, T=1.0)
│ ├── generate_chain_length_p1winners.py # Set D generator (short timescale)
│ ├── generate_chain_length_p1winners_med.py
│ ├── generate_chain_length_p1winners_long.py
│ ├── analyze_rheology.py # Pre-correction rheology (40 sims)
│ ├── analyze_rheology_p1winners.py # P1Winners Set D short (~14 μs)
│ ├── analyze_rheology_p1winners_med.py # P1Winners Set D medium (~120-210 μs)
│ ├── analyze_rheology_p1winners_long.py # P1Winners Set D long (~250-420 μs)
│ ├── analyze_pores.py # Pre-correction pore analysis
│ ├── analyze_pores_p1winners.py # P1Winners pore analysis (short, KDTree)
│ ├── analyze_pores_p1winners_med.py # P1Winners pore analysis (medium, KDTree)
│ ├── analyze_pores_p1winners_long.py # P1Winners pore analysis (long, KDTree)
│ ├── analyze_rg.py
│ ├── create_chainlength_rheology_figure.py
│ └── run_*.sh # SLURM job scripts
│
├── muc2-pore-char/ # Bhattacharya-Gubbins pore characterization
│ ├── pore_analysis.py
│ └── batch_pore_analysis.py
│
├── thesis_figures/ # Thesis figure generation scripts
│ ├── Figure_4_1_flory_scaling/
│ ├── Figure_5_1_concentration_scaling/
│ ├── Figure_5_2_cutoff_independence/
│ ├── Figure_5_3_pore_size_distribution/
│ ├── Figure_5_4_pore_distributions/
│ ├── Figure_S1_temperature_dependence/
│ ├── Figure_S2_parameter_comparison/
│ ├── Figure_chain_length_Rg/
│ ├── Figure_chain_length_pores/
│ ├── Figure_chain_length_scaling/
│ └── Figure_stress_comparison/
│
├── figure_scripts/ # Volume-corrected thesis figure generators
│ ├── generate_chainlength_figures.py # Chain length G', tan(δ) (hardcoded corrected data)
│ ├── generate_chainlength_setD_figures.py # Set D scaling and pore sizes (hardcoded corrected data)
│ ├── generate_phase2_corrected_figures.py # Phase 2 G', G'', tan(δ) panels
│ ├── generate_pore_distributions_setD.py # Set D pore CDFs with SE bands
│ ├── generate_pore_histograms_setD.py # Set D pore histograms
│ └── generate_pore_validation.py # N=52 pore validation histogram
│
├── docs/ # Documentation
│ └── LAMMPS_TIPS.md
│
└── appendix_code/ # Supplementary/appendix code
├── analyze_rheology.py
├── compute_full_rheology.py
├── create_scaling_figures.py
├── generate_replicate_table.py
├── data_tables.tex
└── examples/
Based on Ford et al. (2021), MUC2 dimers are represented as bead-spring chains:
| Model | Beads/dimer | Pattern | Description |
|---|---|---|---|
| Pattern-Preserving | 6 | T-G-T x 2 | Minimal domain structure |
| Ratio-Preserving | 18 | T2-G6-T1 x 2 | Domain ratio-preserving |
| Full (half-length) | 52 | T5-G18-T3 x 2 | Half-length physiological dimer |
| Extended | 104 | T11-G35-T6 x 2 | Full MUC2 dimer |
Bead types:
- Type 1 (T): Hydrophobic terminal/cysteine-knot domains
- Type 3 (G): Glycosylated central domains
Force field: FENE bonds (K=30, R0=1.5) with Lennard-Jones interactions (cutoff 2.5 sigma).
Physical mapping: sigma = 9 nm, epsilon = k_BT at 298 K, tau = 1403.85 ps.
- Systematic variation of epsilon_TT, epsilon_TH, epsilon_HH, temperature, LJ cutoff, and concentration
- Short equilibration (50k steps) + production (600k steps, ~4.2 μs)
- Note: Original runs used WCA (repulsive-only) interactions due to a pair_coeff cutoff bug. Corrected runs in
phase1_corrected/use full LJ. Seedocs/LAMMPS_TIPS.mdfor details. - Median tan(delta) = 3.1, 87.9% sol-like; strongest elastic response at epsilon_TT=0.7, T=0.9
- Four parameter sets (epsilon_TT = 0.7, 0.8, 0.9, 1.0) x 11 concentrations (10-60 mg/mL) x 3 replicates
- Fixed: T=0.9, epsilon_TH=0.4, epsilon_HH=0.4, rc=2.0 sigma
- All conditions liquid-like (tan(delta) >> 1); G'_GK = 40-175 Pa
- Set D (epsilon_TT=1.0) selected for chain length study
- N = 6, 18, 52, 104 beads/dimer at 10 and 30 mg/mL, 5 replicates each
- Uses pre-correction parameters (epsilon_TT=0.4, T=1.0) identified before WCA bug fix, with correct full LJ
- G' ~ N^2.3 scaling, pore sizes 70-200 nm (consistent with experiment)
- Same chain lengths and concentrations as pre-correction study
- Uses optimized Set D parameters (epsilon_TT=1.0, T=0.9, rc=2.0)
- Three timescales: short (~14 μs, dt=0.005), medium (~120-210 μs, dt=0.012), long (~250-420 μs, dt=0.012)
- Medium and long runs partially completed due to Narval 4-day/7-day walltime limits; actual production time varies by system size (larger chains complete fewer steps)
- All conditions liquid-like at all timescales
Each figure has a dedicated generation script in thesis_figures/:
| Thesis Figure | Script | Data Source |
|---|---|---|
| Fig 4.1 (Flory scaling) | thesis_figures/Figure_4_1_flory_scaling/create_figure_4_1_improved.py |
Single-chain validation sims |
| Fig 5.1 (Concentration scaling) | thesis_figures/Figure_5_1_concentration_scaling/create_figure_5_1.py |
Phase 1 parameter sweep |
| Fig 5.2 (Cutoff independence) | thesis_figures/Figure_5_2_cutoff_independence/create_figure_5_2_updated.py |
Phase 1 parameter sweep |
| Fig 5.3 (Pore validation) | thesis_figures/Figure_5_3_pore_size_distribution/create_figure_5_3_chainlength_data.py |
Chain length pore analysis |
| Fig 5.4 (Pore distributions) | thesis_figures/Figure_5_4_pore_distributions/create_figure_pore_distributions_with_zoom.py |
Chain length pore analysis |
| Fig S1 (Temperature) | thesis_figures/Figure_S1_temperature_dependence/create_figure_temperature.py |
Phase 1 parameter sweep |
| Fig S2a (SE vs SD) | thesis_figures/Figure_S2_parameter_comparison/create_S2b_SE_vs_SD_comparison.py |
Phase 1 parameter sweep |
| Fig S2b (Param effects) | thesis_figures/Figure_S2_parameter_comparison/create_parameter_effects.py |
Phase 1 parameter sweep |
| Fig S2c (Phase 2 comparison) | thesis_figures/Figure_S2_parameter_comparison/create_figure_phase2_3panel.py |
Phase 2 results |
| Fig S2d (Sensitivity) | thesis_figures/Figure_S2_parameter_comparison/create_sensitivity_analysis.py |
Phase 1 parameter sweep |
| P1Winners concentration | phase2_analysis/create_phase2_figures.py |
Phase 2 P1Winners |
| P1Winners pore analysis | phase2_analysis/create_pore_thesis_figure.py |
Phase 2 P1Winners |
| Chain length rheology | chain_length_study/create_chainlength_rheology_figure.py |
Chain length study |
| Chain length Rg | thesis_figures/Figure_chain_length_Rg/create_figure_chain_length_Rg_updated.py |
Chain length study |
| Chain length pores | thesis_figures/Figure_chain_length_pores/create_figure_chain_length_pores.py |
Chain length study |
| Chain length scaling | thesis_figures/Figure_chain_length_scaling/create_figure_chain_length_scaling_updated.py |
Chain length study |
| Stress comparison | thesis_figures/Figure_stress_comparison/create_figure_stress_comparison.py |
Phase 1 vs Phase 2 |
| Chain length G' vs N (corrected) | figure_scripts/generate_chainlength_figures.py |
Volume-corrected Set D data |
| Set D scaling + pores | figure_scripts/generate_chainlength_setD_figures.py |
Volume-corrected Set D data |
| Phase 2 corrected panels | figure_scripts/generate_phase2_corrected_figures.py |
Phase 2 rheology JSON |
| Set D pore CDFs | figure_scripts/generate_pore_distributions_setD.py |
Chain length pore data |
| Set D pore histograms | figure_scripts/generate_pore_histograms_setD.py |
Chain length pore data |
| Pore validation (N=52) | figure_scripts/generate_pore_validation.py |
Chain length pore data |
numpy >= 1.20
scipy >= 1.7
matplotlib >= 3.4
pandas >= 1.3
MDAnalysis >= 2.0 # For pore characterization only
See requirements.txt for pip installation, or appendix_code/requirements.txt for the same list.
# Chain length study (corrected, full LJ)
cd chain_length_study
python generate_chain_length_sims.py
# Phase 1 corrected
cd phase1_corrected
python generate_phase1_narval.py --batch-id 0 --count-only # Check total count
# Phase 2 corrected
cd parameter_sweep
python phase2_corrected_generator.py# Transfer to cluster
rsync -avz <sim_directory>/ user@narval.computecanada.ca:~/scratch/<sim_directory>/
# Submit jobs
cd ~/scratch/<sim_directory>
bash submit_all.sh# Rheological analysis
cd chain_length_study
python analyze_rheology.py
# Pore characterization
cd muc2-pore-char
python pore_analysis.py <trajectory_file> --lammps-data <data_file>
# Phase 2 corrected analysis
cd parameter_sweep
python analyze_phase2_corrected.pycd thesis_figures/Figure_5_1_concentration_scaling
python create_figure_5_1.pyLAMMPS LJ units to SI (sigma = 9.0 nm):
- Flory Scaling: nu = 0.475 +/- 0.001 (near-ideal chain behavior in implicit solvent)
- WCA Bug Discovery: pair_coeff cutoff 1.1225 sigma overrides pair_style cutoff 2.5 sigma, eliminating all attractive interactions. Fix: explicitly set cutoff in each pair_coeff line. See
docs/LAMMPS_TIPS.md. - Liquid-Like Behavior: All tested conditions (Phase 1, Phase 2 P1Winners, chain length study) are sol-like with tan(delta) >> 1 from consistent Green-Kubo analysis
- Hydrophobic Strength: Increasing epsilon_TT from 0.7 to 1.0 increases G' but system remains liquid-like (tan(delta) >> 1 at all concentrations)
- Pore Structure: Chain length study pores (70-200 nm) consistent with AFM/SEM measurements of MUC2 mucus
- Chain Length Scaling: G' ~ N^2.3, Rg ~ N^0.5 (polymer physics validated across N = 6 to 104)
- Timescale Independence: P1Winners Set D results consistent across short (~14 μs), medium (~120-210 μs), and long (~250-420 μs) production runs — liquid-like behavior is not a short-timescale artifact
Compute Canada - Narval Cluster
- Account: def-rmansbac
- LAMMPS: lammps-omp/20230802
- Typical job: 4 cores, 16 GB RAM, 6-24 hours
- Ford, A. G., et al. (2021). Molecular dynamics simulations to explore the structure and rheological properties of normal and hyperconcentrated airway mucus. Studies in Applied Mathematics, 147(4), 1369-1387.
- Kremer, K. & Grest, G. S. (1990). Dynamics of entangled linear polymer melts: A molecular-dynamics simulation. J. Chem. Phys., 92(8), 5057-5086.
- Bhattacharya, S. & Gubbins, K. E. (2006). Fast method for computing pore size distributions of model materials. Langmuir, 22, 7726-7731.
- Parsons, J. M. (2021). A Computational Model of Mucus. M.S. thesis, San Diego State University.
- Round, A. N., et al. (2012). Lamellar structures of MUC2-rich mucin: a potential role in governing the barrier and lubricating functions of intestinal mucus. Biomacromolecules, 13(10), 3253-3261.
- Lai, S. K., et al. (2009). Micro- and macrorheology of mucus. Advanced Drug Delivery Reviews, 61(2), 86-100.
Academic use only. Contact author for permissions.