How paracetamol, dipyrone metabolites and arachidonic acid interact with the two active sites of cyclooxygenase-2
MRes research project · all-atom MD of the COX-2 homodimer
Overview • Systems • Workflow • Results • Quick start • Layout • Reading
Cyclooxygenase-2 (COX-2) is a homodimer, and each monomer has two catalytic sites:
|
The cyclooxygenase channel, where the substrate arachidonic acid binds and Tyr385 starts catalysis. Classical NSAIDs and coxibs such as rofecoxib block this site. |
The peroxidase site at the heme. One proposed mechanism is that paracetamol and dipyrone metabolites act here as reducing co-substrates rather than blocking the channel. |
This project uses all-atom molecular dynamics to test where these drugs sit and what they do to the enzyme. For every system it compares:
| 📏 Binding stability | 🌀 Protein dynamics | 🔗 Allostery | ⚡ Energetics |
|---|---|---|---|
| distance to Tyr385 or the heme iron | RMSD · RMSF · radius of gyration | residue correlation × contact networks between monomers | MM/PBSA and MM/GBSA binding energies |
| Ligands | 6: two references and four analgesics/metabolites |
| Binding sites | COX site, POX site, and both at once (1C_2P) |
| Sampling | 300 ns per replica, 2–3 replicas per system |
| Force field | ff19SB · GAFF + RESP · MCPB.py heme/Fe model |
| Solvent | OPC or TIP3P water, Na⁺/Cl⁻ |
| Time step | 4 fs with hydrogen mass repartitioning |
| Ligand | What it is | Role in this study | |
|---|---|---|---|
| 🧪 | AA: arachidonic acid | Natural COX-2 substrate | Reference. Productive (tail-up) and non-productive (head-up) poses |
| 🛑 | Vioxx: rofecoxib | Selective COX-2 inhibitor | Reference for a COX-site blocker |
| 💊 | APAP: paracetamol | Analgesic / antipyretic | COX site and POX site |
| 🔁 | AM404: N-arachidonoylphenolamine | Paracetamol metabolite | COX site and POX site |
| 💉 | MAA: 4-methylaminoantipyrine | Active metabolite of dipyrone | COX site and POX site |
| 🔁 | 4-AA: 4-aminoantipyrine | Metabolite of dipyrone | COX site and POX site |
Tip
Naming convention: MAA_POX means MAA at the POX site. 1C_2P systems have two copies of the ligand, one in a COX site and one at a POX site. rep2/rep3 are independent replicas.
flowchart LR
A["🧱 COX-2 dimer<br/>crystal structure"] --> B["⚛️ Heme / Fe site<br/>MCPB.py"]
L["💊 Ligand"] --> Q["🧮 QM: B3LYP opt<br/>HF/6-31G* ESP"]
Q --> R["🔋 RESP charges<br/>GAFF params"]
B --> S["🏗️ tleap<br/>solvate + ions"]
R --> S
S --> H["⚖️ ParmEd<br/>H-mass repartition"]
H --> M["🔥 min → heat →<br/>10 ns equilibration"]
M --> P["🚀 Production<br/>6 × 50 ns"]
P --> T["📈 cpptraj + pytraj<br/>RMSD · RMSF · distances"]
P --> E["⚡ MMPBSA.py<br/>binding energy"]
T --> F["📊 Summary figures"]
E --> F
| Step | Stage | Details | Folder |
|---|---|---|---|
| 1 | Ligand parameters | Gaussian optimisation → RESP charges → GAFF (antechamber, parmchk2) |
01_ligand_parameterisation |
| 2 | System setup | ff19SB protein + MCPB.py heme/Fe + docked ligand, solvated in tleap, HMR with ParmEd | 02_system_setup |
| 3 | Minimisation | 10 000 steps restrained, then 16 000 unrestrained | 03_md_protocol |
| 4 | Heating | 0 → 100 K (NVT), then 100 → 320 → 300 K (NPT), restrained | 03_md_protocol |
| 5 | Equilibration | 10 ns, 300 K, NPT, Monte Carlo barostat | 03_md_protocol |
| 6 | Production | 50 ns segments (300 ns per replica), frames every 50 ps | 03_md_protocol |
| 7 | Analysis | cpptraj stripping, pytraj notebooks, MMPBSA.py array jobs | 04_analysis · 05_project_summary |
Two figures from the project's summary notebooks, as saved from the original run (legends kept as originally produced).
|
Arachidonic acid reference poses 
The productive (tail-up) pose stays closer to Tyr385 than the non-productive (head-up) pose. |
Drugs and metabolites in the COX site 
Ligand–Tyr385 distance for APAP, 4-MAA, 4-AA and AM404, 1 ns averages. |
More figures (RMSD/Rg joint plots, RMSF profiles, MM/PBSA traces) are saved in the notebooks under 04_analysis and 05_project_summary.
Warning
Corrections to the original project. Check these before reusing figures from the original work:
- 🔥 Heating protocol:
Heat_2.innever ran its final 320 → 300 K ramp because of an extra&wt type='END'(details). Now fixed. - 📉 POX-site figure: in
05_project_summary/POX_site.ipynb, the two AM404 lines of the ligand–heme distance plot were drawn from 4-AA data. The code is fixed; re-run the notebook to regenerate the figure. - ⚡ MM/PBSA units: two notebooks labelled the energies kJ/mol; MMPBSA.py reports kcal/mol.
- 📐 RMSF subplots: the summary notebooks had the x and y axis labels swapped.
# 1. Environment (AmberTools, pytraj, Jupyter, ...)
conda env create -f environment.yml
conda activate cox2-md
# 2. Build a system (example: paracetamol at the POX site)
cd 02_system_setup/APAP_POX
tleap -f APAP_POX.in
parmed -p APAP_POX.parm7 -i ../hmr.parmed # set outparm to APAP_POX_HMR.parm7
# 3. Simulate on a local GPU: NAME, last solute residue, GPU id
../../03_md_protocol/run_equilibration.sh APAP_POX 1107 0
../../03_md_protocol/run_production.sh APAP_POX 6 0 # 6 × 50 ns = 300 ns🖥️ Running on a PBS cluster instead
# Equilibration
qsub -v NAME=APAP_POX,RESNUM=1107,PROTOCOL_DIR=/path/to/03_md_protocol \
/path/to/03_md_protocol/hpc/equilibration.pbs
# Production: a self-resubmitting chain, one 50 ns segment per job
python /path/to/03_md_protocol/hpc/next_job.py APAP_POX 1 0 --init --nseg 6
# MM/PBSA over 6000 frames split into 200 array tasks
qsub -J 1-200 -v JOB=APAP_POX,CHUNK=30 04_analysis/mmpbsa/mmpbsa_array.pbsSee 03_md_protocol/README.md for details.
Note
pmemd.cuda is licensed separately from AmberTools. Amber changes between releases, so check the current manual and tutorials before reusing the inputs.
The folders follow the order of the workflow. Each one has its own README.
📦 MRes_Project_COX-2
├── 📂 01_ligand_parameterisation QM optimisation, RESP charges, GAFF parameters per ligand
├── 📂 02_system_setup tleap scripts: protein + heme (MCPB.py) + ligand → solvated system
│ ├── common/ shared COX-2 model and heme/Fe parameters
│ └── AA_COX/ APAP_COX/ APAP_POX/ MAA_POX/
├── 📂 03_md_protocol pmemd inputs + run scripts (local GPU and PBS)
├── 📂 04_analysis
│ ├── trajectory/ per-system pytraj notebooks + cpptraj/parmed inputs
│ └── mmpbsa/ MMPBSA.py input, PBS array script, result notebooks
├── 📂 05_project_summary cross-system comparison notebooks (main figures)
├── 📂 docs/figures figures shown in this README
└── 📄 environment.yml conda environment for the analysis
Important
Trajectories and topologies (.nc, .parm7, …) are too large for git and are not included. The notebooks are kept with their outputs so the results can still be viewed. Each notebook opens with a header describing the system and any corrections made since its outputs were produced.
| Tool | Used for |
|---|---|
| Amber / AmberTools | tleap, antechamber, parmchk2, MCPB.py, ParmEd, pmemd.cuda, cpptraj, MMPBSA.py |
| Gaussian 16 / GaussView | Ligand geometry optimisation and electrostatic potential |
| Python | pytraj, NumPy, SciPy, pandas, Matplotlib, seaborn, Jupyter |
🧬 COX-2 review
- Chem. Rev., doi:10.1021/acs.chemrev.0c00215: a comprehensive starting point