Dissecting Cleavage Mechanisms
in Tensed Collagen Fibrils
using Reactive Molecular Dynamics Simulations

Jannik Buhr

Introduction

Concerning hobbits collagen

An orange a day keeps the scurvy away

Citrus fruits contain Vitamin C ( 🟠 )

1747, HMS Salisbury, James Lind performs the first clinical trial with a control group [1]

Collagen keeps us together

• single peptide

• triple helix

• collagen fibril

The Gly-X-Y motif allows tight turns and packing

Glycine (Gly)

Proline (Pro)

Hydroxyproline (Hyp)

X,Y: proline or hydroxyproline

Prolyl hydroxylase ( 🟠 )

Collagen I is cross-linked

 

 

 

Lysine

Lysyl-hydroxylase ( 🟠 )
→

Hydroxylysine

Lysyl oxidase
→

Center of a pyridinoline crosslink

 

A rather radical discovery: mechanoradicals in collagen

“Mechanoradicals in tensed tendon collagen as a source of oxidative stress”, Zapp et al. [3]

Goals

for this talk, and by extension, my thesis

Dissecting cleavage mechanisms in tensed collagen fibrils using reactive molecular dynamics simulations

“Mechanical Activation Drastically Accelerates Amide Bond Hydrolysis,
Matching Enzyme Activity”, Pill et al., 2019 [4]

Hydrolysis
vs.
Homolysis

What factors influence their competition?

What is the influence of the collagen structure?

Why do we still see radicals?

Dissecting cleavage mechanisms in tensed collagen fibrils using reactive molecular dynamics simulations

QM/MM simulations of
base-catalyzed hydrolysis

 

Kinetic Monte Carlo Molecular Dynamics:
KIMMDY

Methods and Theory

Molecular dynamics / Molecular Mechanics

How do we make the atoms move?

\[ \begin{array}{rcl} F &=& ma \\ -\frac{dV}{dr} &=& m\frac{d^2r}{dt^2}, \end{array} \tag{1}\]

\(V\): potential
\(r\): position
\(F\): force
\(a\): acceleration
\(t\): time
\(m\): mass

Force Fields and Topologies

Bond \(\mathrm{r}\), angle \(\mathrm{\alpha}\), dihedral \(\mathrm{\phi}\), non-bonded terms \(\mathrm{nb}\), \(\mathrm{1-4}\).

Morse potential.

Harmonic potential.

Lennard-Jones potential.

Coulomb potential.

Connectivity and parameters pre-defined!

So how do we do chemistry?

Quantum mechanics/molecular mechanics

Electrostatic embedding scheme for QM/MM. Figure by Dmitry Morozov [5]

DFT: Density Functional Theory

Collagen Hydrolysis by Quantum Mechanics/Molecular Mechanics

Base-catalyzed hydrolysis

TS: transition state
TI: Tetrahedral Intermediate
ZI: Zwitterionic Intermediate

Previous work on base-catalyzed hydrolysis

Figure from [4]. From top to bottom: educt, TS1, TI.

Data from Pill et al. [4] visualized. Atomic Force Microscopy (experimental) and ab initio QM (theoretical) energies

System Setup

[]: input parameters.

QM region selection and setup

\(BD\): Bürgi-Dunitz angle, 107°

\(FL\): Flippin-Lodge angle, 0°

Umbrella sampling the reaction coordinate

QM/MM simulations reveal mechanistic details of the tetrahedral intermediate formation

QM/MM simulations reveal mechanistic details of the tetrahedral intermediate formation

The triple helix marginally impedes hydroxide attack thermodynamically

\(\Delta\Delta G\) 9.02 ±4.41 kJ/mol (±SEM)

Rate reduction due to triple helix:

30 to 40 times.

Proton mobility in explicit QM solvent molecules allows breaking of the peptide bond

Discussion

• Effect size?
• Survivorship bias?
• Cause?
• Computational cost

Kinetic Monte Carlo Molecular Dynamics: KIMMDY

Kinetic Monte Carlo Molecular Dynamics: KIMMDY

with Eric Hartmann and Kai Riedmiller,
inspired by Benedikt Rennekamp [6]

graeter-group.github.io/kimmdy/

Read now on bioRxiv:
KIMMDY: A biomolecular reaction emulator [7],
or soon in Nature Communications.

KIMMDY is easy to use

KIMMDY is easy to extend to new chemistries

KIMMDY takes care of the complicated topology modifications.

Parametrization interface to obtain parameters after a reaction

Graph Attentional Protein Parametrization by Leif Seute [8]

also a good excuse to show the mol* plugin I wrote for our documentation

Slow-growth to execute even complex reactions like hydrolysis

Does not have to follow the chemically
accurate reaction coordinate

Applying KIMMDY to collagen
hydrolysis and homolysis

Applying KIMMDY to collagen hydrolysis and homolysis

Applying KIMMDY to collagen hydrolysis and homolysis

Arrhenius equation

\[ k = A\,\mathrm{e}^{\left(\frac{-\mathbf{E_\mathrm{a}}}{RT}\right)} \]

\(A\): pre-exponential factor (attempt frequency)
\(R\): gas constant
\(T\): temperature

Morse potential

---

Homolysis

Physical (Bell-Evans [9, 10]) model
by Benedikt Rennekamp [6]

Hydrolysis

Hybrid physical + experimental model
based on AFM data by Pill et al. [4]

Hydrolysis rate incorporates SASA

SASA

Hydrolysis rate incorporates SASA

Szasza SASA = Solvent Accessible Surface Area

Hydrolysis rate incorporates pH, SASA and force

\[ \mathrm{Rate}(bond, t) = \ \frac{\mathrm{SASA}(bond, t)}{\mathrm{SASA}_{max}} \cdot \ \frac{c_{OH^-}}{c_{OH^-_{\mathrm{exp}}}} \ \cdot \mathrm{rate}_{\mathrm{exp}}(F_{bond,t}) \tag{2}\]

Solvent accessibility in the fibril is heterogeneous

 

 

 

Peptide bonds near the center are still accessible

 

 

Homolysis becomes competitive with hydrolysis by force concentration

 

 

Discussion and outlook

• pH just a factor, could be e.g. constant pH simulations [11]
• SASA influence likely to be underestimated
• Need a link to real tissue experiments

Concluding Thoughts

Concluding Thoughts

• KIMMDY is a versatile tool for reactive MD simulations:
  graeter-group.github.io/kimmdy/

• Triple helix marginally impedes hydrolysis.
  Homolysis becomes competitive to hydrolysis due to force concentration

• Thesis and slides available at
  jmbuhr.de/phd-thesis

• Eat your oranges! 🟠

Thank You!

Thank You!

Questions?

• jmbuhr.de
• jmbuhr
• jmbuhr
• mpip-mainz.mpg.de/en/graeter
• graeter-group

This work was supported by the Klaus Tschira Foundation and has received funding from the European Research Council (ERC).

References

1.

Lind, J. A Treatise on the Scurvy: In Three Parts. Containing an Inquiry Into the Nature, Causes, and Cure, of That Disease. Together with a Critical and Chronological View of What Has Been Published on the Subject. S. Crowder. (S. Crowder, 1772).

2.

The University of Wales Bioimaging laboratory, Institute of Biological Sciences, The University of Wales, Aberystwyth, Wales, UK, Hughes, L., Archer, C. & Ap Gwynn, I. The ultrastructure of mouse articular cartilage: Collagen orientation and implications for tissue functionality. A polarised light and scanning electron microscope study and review. European Cells and Materials 9, 68–84 (2005).

3.

Zapp, C., Obarska-Kosinska, A., Rennekamp, B., et al. Mechanoradicals in tensed tendon collagen as a source of oxidative stress. Nature Communications 11, 2315 (2020).

4.

Pill, M. F., East, A. L. L., Marx, D., Beyer, M. K. & Clausen‐Schaumann, H. Mechanical Activation Drastically Accelerates Amide Bond Hydrolysis, Matching Enzyme Activity. Angewandte Chemie International Edition 58, 9787–9790 (2019).

5.

Morozov, D. Webinar: Multiscale QM/MM simulations: Exploring chemical reactions using novel GROMACS/CP2K interface (2020-12-08). (2020).

6.

Rennekamp, B., Kutzki, F., Obarska-Kosinska, A., Zapp, C. & Gräter, F. Hybrid Kinetic Monte Carlo/Molecular Dynamics Simulations of Bond Scissions in Proteins. Journal of Chemical Theory and Computation 16, 553–563 (2020).

7.

Buhr*, J., Hartmann*, E., Riedmiller*, K., et al. KIMMDY: A biomolecular reaction emulator. bioRxiv. (2025) doi:10.1101/2025.07.02.662624 *These authors contributed equally to this work.

8.

Seute, L., Hartmann, E., Stühmer, J. & Gräter, F. Grappa – a machine learned molecular mechanics force field. Chemical Science 16, 2907–2930 (2025).

9.

Bell, G. I. Models for the Specific Adhesion of Cells to Cells. Science 200, 618–627 (1978).

10.

Evans, E. & Ritchie, K. Dynamic strength of molecular adhesion bonds. Biophysical Journal 72, 1541–1555 (1997).

11.

Aho, N., Buslaev, P., Jansen, A., et al. Scalable Constant pH Molecular Dynamics in GROMACS. Journal of Chemical Theory and Computation 18, 6148–6160 (2022).

12.

Procida, D. Diátaxis documentation framework.

Appendix

Force Fields and Topologies

Bond \(\mathrm{r}\), angle \(\mathrm{\alpha}\), dihedral \(\mathrm{\phi}\), non-bonded terms \(\mathrm{nb}\), \(\mathrm{1-4}\).

Listing 1: Select lines from ffnonbonded.itp

[ atomtypes ]
; name at.num  mass    charge ptype  sigma      epsilon
C       6      12.01   0.0000  A   3.39967e-01  3.59824e-01
CA      6      12.01   0.0000  A   3.39967e-01  3.59824e-01
H       1       1.008  0.0000  A   1.06908e-01  6.56888e-02

Listing 2: Select simplified lines from a .top file of a capped glycine dipeptide

[ moleculetype ]
; Name            nrexcl
Protein             3

[ atoms ]
; nr  type resnr residue  atom  cgnr   charge  mass  typeB ..
  1    CT     1    ACE    CH3    1  -0.3662   12.01
  2    HC     1    ACE   HH31    2   0.1123   1.008

[ bonds ]
;  ai    aj funct    c0    c1      c2    c3
    1     2     1
    1     3     1

[ angles ]
;  ai    aj    ak funct   c0   c1   c2    c3
    2     1     3     1

More SASA

KIMMDY finds a new radical scavanging candidate in collagen

DOPA

Pyridinoline

Proton mobility in explicit QM solvent molecules allows breaking of the peptide bond

 

KIMMDY is easy to install and well documented

graeter-group.github.io/kimmdy/

Running KIMMDY: Other cool features

• checkpoints and restarting
• high performance computing support
• extensive testing and logging
• easy to install, use and extend

Direct comparison of hydrolysis and homolysis force response

 

The unique effect of prolines

 

 

 

Supplementary Information

Small multiples plot of the approach of the hydroxide in terms of distances of the carbonyl \(\ce{O}\) of the TI to nearest proton of either a MM solvent molecule (cyan), the protein (black) or the QM water (magenta) across all umbrella windows. Distances are average within each window and the windows belonging to the same approach simulation are connected. The thick lines are a LOESS curve across all windows.

Supplementary Information

Figure 1: The protonation states are show as the average number of protons assigned to either the hydroxyl \(\ce{O}\) or the carbonyl \(\ce{O}\) with all sampling windows averaged for each system. The region for detection of TS1 is shaded cyan, while the region for detection of the TI is shaded gray.