Exploring Caspase Mutations and Post-Translational Modification by Molecular Modeling Approaches

Summary

The present protocol uses a biomolecular simulation package and describes the molecular dynamics (MD) approach for modeling the wild-type caspase and its mutant forms. The MD method allows for assessing the dynamic evolution of the caspase structure and the potential effect of mutations or post-translational modifications.

Abstract

Apoptosis is a type of programmed cell death that eliminates damaged cells and controls the development and tissue homeostasis of multicellular organisms. Caspases, a family of cysteine proteases, play a key role in apoptosis initiation and execution. The maturation of caspases and their activity is fine-tuned by post-translational modifications in a highly dynamic fashion. To assess the effect of post-translational changes, potential sites are routinely mutated with residues persistent to any modifications. For example, the serine residue is replaced with alanine or aspartic acid. However, such substitutions could alter the caspase active site's conformation, leading to disturbances in catalytic activity and cellular functions. Moreover, mutations of other amino acid residues located in critical positions could also break the structure and functions of caspases and lead to apoptosis perturbation. To avoid the difficulties of employing mutated residues, molecular modeling approaches can be readily applied to estimate the potential effect of amino acid substitutions on caspase structure. The present protocol allows the modeling of both the wild-type caspase and its mutant forms with the biomolecular simulation package (Amber) and supercomputer facilities to test the effect of mutations on the protein structure and function.

Introduction

Apoptosis is one of the most widely-studied cellular processes that regulate the morphogenesis and tissue homeostasis of multicellular organisms. Apoptosis can be initiated by a wide range of external or internal stimuli, such as the activation of death receptors, disturbance in the cell cycle signals, DNA damage, endoplasmic reticulum (ER) stress, and various bacterial and viral infections¹. The caspases - key apoptotic players - are conventionally classified into two groups: initiators (caspase-2, caspase-8, caspase-9, and caspase-10) and effectors (caspase-3, caspase-6, and caspase-7), depending on their domain structure and the place in the caspase cascade^2,3. Upon cell death signals, the initiator caspases interact with adaptor molecules that facilitate proximity-induced dimerization and autoprocessing to form an active enzyme. The effector caspases are activated through cleavage by initiator caspases and perform downstream execution steps by cleaving multiple cellular substrates⁴.

The maturation and function of the initiator and effector caspases are regulated by a large number of different intracellular mechanisms, among which the post-translational modification plays an indispensable role in cell death modulation⁵. The addition of modifying groups (phosphorylation, nitrosylation, methylation, or acetylation) or proteins (ubiquitination or SUMOylation) changes the enzymatic activity of caspases or the protein conformation and stability that regulate apoptosis. Site-directed mutagenesis is widely applied to investigate the potential post-translational modification sites and discern their role. A putative modification site is usually replaced by another amino acid, which cannot be further modified. Thus, potentially phosphorylated serine and threonine are mutated to alanine, and lysine ubiquitination sites are replaced with arginine. Another strategy includes substituting an amino acid that particularly mimics post-translational modification (e.g., glutamate and aspartate have been used to mimic phosphorylated serine or threonine)⁶. However, some of these substitutions located in the high vicinity of an active site or in critical positions could change caspase structure, disturb catalytic activity, and suppress apoptotic cell death⁷. Similar effects could be observed in cases of tumor-associated missense mutations in caspase genes. For example, the tumor-associated mutation of caspase-6 - R259H - resulted in conformational changes of loops in the substrate-binding pocket, reducing the efficient catalytic turnover of substrates⁸. The G325A amino acid substitution in caspase-8 identified in head and neck squamous cell carcinoma could hamper caspase-8 activity, which led to the modulation of nuclear factor-kB (NF-kB) signaling and promoted tumorigenesis⁹.

To assess the potential effect of amino acid substitutions on caspase structure and function, molecular modeling can be applied. The molecular dynamics (MD) approach is described in this work for modeling the wild-type caspase and its mutant forms using the biomolecular simulation package (Amber). The MD method gives a view of the dynamic evolution of the protein structure following the introduction of mutations. Originally developed by Peter Kollman's group, the Amber package became one of the most popular software tools for biomolecular simulations^10,11,12,13. This software is divided into two parts: (1) AmberTools, a collection of programs routinely used for system preparation (atom type assignment, adding hydrogens and explicit-water molecules, etc.) and trajectory analysis; and (2) Amber, which is centered around the pmemd simulation program. AmberTools is a free package (and a prerequisite for installing Amber itself), while Amber is distributed with a separate license and fee structure. Parallel simulations on a supercomputer and/or using graphics processing units (GPUs) can substantially improve the performance for the scientific research of protein structure dynamics¹⁴. The latest available software versions are AmberTools21 and Amber20, but the described protocols can also be used with the former versions.

Protocol

1. System preparation

NOTE: The molecular models of the native and mutant protein forms are built based on an appropriate crystal structure obtained from the Protein Data Bank^15,16.

1. To retrieve the selected PDB structure, use the Download Files drop-down list and click on PDB Format. Remove remarks and connectivity data, and insert a TER card between separate protein chains in the PDB file. Optionally, set the protonation state of histidine residues by replacing the residue name HIS with HIE, HID, or HIP (Supplementary File 1).
   NOTE: It is reasonable not to remove crystallographically resolved water molecules (if present in the PDB file).
2. To prepare the starting model, launch the tleap program (AmberTools package) and then enter commands that manipulate tleap objects.
   1. Launch the program: tleap. Load the ff14SB force field to describe the protein with molecular mechanics: > source leaprc.protein.ff14SB.
   2. Load force fields for water molecules and atomic ions (Na⁺, Cl^-): > source leaprc.water.tip3p.
   3. Load the PDB file and build coordinates for hydrogens, creating an object named mol: > mol = loadpdb protein.pdb.
   4. Check for internal inconsistencies that could cause problems: > check mol.
      NOTE: Disulfide bridges, if present, must be specified manually. Replace the names of covalently bound cysteines CYS with CYX and type the following command in tleap to create a bond between SG atoms: bond mol.X.SG mol.Y.SG (X and Y are cysteine residue numbers).
   5. Create a solvent box around the protein (TIP3P water, 12 Å distance): > solvatebox mol TIP3PBOX 12.0.
   6. Check the total charge: > charge mol, and add counterions (Na⁺ or Cl^-) to neutralize the system:
      > addions mol Na⁺ 0.
   7. Create the topology prmtop file and the coordinate inpcrd file (inputs for running a simulation): > saveamberparm mol protein.prmtop protein.inpcrd. Quit the tleap program: > quit.
      ​NOTE: Amber-format coordinate files can be easily converted to PDB format using the ambpdb tool (refer to the Users' Manual, see Table of Materials).

2. Energy minimization

NOTE: The energy minimization is necessary to remove any bad contacts and overlaps between atoms in the starting system that lead to instability when running MD.

1. Carry out the first stage of the energy minimization (2,500 steps of the steepest descent algorithm + 2,500 steps of the conjugate gradient algorithm) to optimize the positions of added hydrogen atoms and water molecules while keeping the protein coordinates fixed by positional restraints on heavy atoms.
   1. Run the pmemd program as follows:
      pmemd -O -i min1.in -p protein.prmtop -c protein.inpcrd -o protein_min1.out
      -r protein_min1.rst -ref protein.inpcrd.
   2. Follow the required arguments: -i FILE control data; -p FILE molecular topology, force field parameters, atom names; -c FILE initial coordinates; -o FILE user-readable log output; -r FILE final coordinates; -ref FILE reference coordinates for position restraints.
   3. Specify options in the input file min1.in:
      &cntrl
      imin=1, maxcyc=5000, ncyc=2500,
      cut=10.0, ntb=1,
      ntc=1, ntf=1,
      ntpr=10,
      ntr=1,
      restraintmask=':1-517 & !@H=',
      restraint_wt=2.0
      /
      NOTE: Input options: imin=1 perform minimization; maxcyc=5000 maximum number of minimization cycles; ncyc=2500 minimization method is switched from steepest descent to conjugate gradient after ncyc cycles; cut=10.0 non-bonded cutoff (Å); ntb=1 periodic boundaries are imposed, constant volume; ntc=1 no constraints are applied to bond lengths (for better energy convergence); ntf=1 all interactions are calculated; ntpr=10 energy information is printed to the out file every ntpr steps; ntr=1 flag for restraining specified atoms using a harmonic potential; restraintmask=':1-517 & !@H=' specifies restrained atoms; restraint_wt=2.0 weight (kcal/mol∙Ų) for positional restraints.
2. Carry out the second stage of the energy minimization (5,000 steepest descent steps + 5,000 conjugate gradient steps) without restraints to optimize the whole system.
   1. Use min2.in and protein_min1.rst as inputs: pmemd -O -i min2.in -p protein.prmtop
      -c protein_min1.rst -o protein_min2.out -r protein_min2.rst.
   2. Specify options in the input file min2.in:
      &cntrl
      imin=1, maxcyc=10000, ncyc=5000,
      cut=10.0, ntb=1,
      ntc=1, ntf=1,
      ntpr=10
      /

3. Heating

NOTE: This stage aims to heat the system from 0 K to 300 K. Initial velocities are assigned to atoms since the starting model based on the PDB file does not contain velocity information.

1. Carry out the heating process with positional restraints on the protein atoms (50 ps, constant volume).
   1. Use heat.in and protein_min2.rst as inputs: pmemd -O -i heat.in -p protein.prmtop
      -c protein_min2.rst -o protein_heat.out
      -r protein_heat.rst -x protein_heat.mdcrd -ref protein_min2.rst
   2. Follow the required arguments: -x FILE coordinate sets saved over MD trajectory.
   3. Specify options in the input file heat.in:
      &cntrl
      imin=0, irest=0, ntx=1,
      nstlim=25000, dt=0.002,
      ntc=2, ntf=2,
      cut=10.0, ntb=1,
      ntpr=500, ntwx=500,
      ntt=3, gamma_ln=2.0,
      tempi=0.0, temp0=300.0,
      ntr=1, restraintmask=':1-517',
      restraint_wt=1.0,
      nmropt=1
      /
      &wt TYPE='TEMP0',
      istep1=0, istep2=25000,
      value1=0.1, value2=300.0 /
      &wt TYPE='END' /
      NOTE: Input options: imin=0 run MD; irest=0 run a new simulation; ntx=1 no initial velocities are read from the rst file; nstlim=25000 number of MD steps; dt=0.002 time step (ps); ntc=2 bonds involving hydrogen are constrained using the SHAKE algorithm; ntf=2 forces for the constrained bonds are not calculated; ntwx=500 coordinates are written to the mdcrd file every ntwx steps; ntt=3 Langevin temperature regulation; gamma_ln=2.0 collision frequency for Langevin dynamics (ps⁻¹); tempi=0.0 initial temperature; temp0=300.0 reference temperature; nmropt=1 parameters specified in an &wt namelist are read; TYPE='TEMP0' vary the target temperature.

4. Equilibration

NOTE: This stage is necessary to adjust the density of water and obtain the equilibrium state of the protein.

1. Perform the equilibration at 300 K without any restraints (500 ps, constant pressure).
   1. Use equil.in and protein_heat.rst as inputs: pmemd -O -i equil.in -p protein.prmtop
      -c protein_heat.rst -o protein_equil.out -r protein_equil.rst -x protein_equil.mdcrd.
   2. Specify options in the input file equil.in:
      &cntrl
      imin=0, irest=1, ntx=5,
      nstlim=250000,
      dt=0.002,
      ntc=2, ntf=2,
      cut=10.0, ntb=2, ntp=1, taup=2.0,
      ntpr=1000, ntwx=1000, ntwr=50000,
      ntt=3, gamma_ln=2.0,
      temp0=300.0
      /
      NOTE: Input options: irest=1 restart the simulation; ntx=5 coordinates and velocities are read from a previously generated rst file; ntb=2 periodic boundaries are imposed, constant pressure; ntp=1 isotropic pressure scaling; taup=2.0 pressure relaxation time (ps); ntwr=50000 every ntwr steps, the rst file is written, ensuring recovery from a crash.

5. Production dynamics

1. After equilibrium has successfully been achieved, carry out a production MD simulation (10 ns or longer) at constant pressure, and generate the trajectory file for subsequent analysis of the protein structure.
   1. Use prod.in and protein_equil.rst as inputs: pmemd -O -i prod.in -p protein.prmtop
      -c protein_equil.rst -o protein_prod.out -r protein_prod.rst -x protein_prod.mdcrd.
   2. Specify options in the input file prod.in:
      &cntrl
      imin=0, irest=1, ntx=5,
      nstlim=5000000,
      dt=0.002,
      ntc=2, ntf=2,
      cut=10.0, ntb=2, ntp=1, taup=2.0,
      ntpr=1000, ntwx=1000, ntwr=50000,
      ntt=3, gamma_ln=2.0,
      temp0=300.0, ig=-1
      /
      NOTE: The parallel version of pmemd (pmemd.MPI) or GPU-accelerated version (pmemd.cuda) can be used on computer clusters and supercomputers. A long MD simulation may be broken into several segments and performed sequentially. It is strongly recommended to set ig=-1 (random seed option) for every simulation restart. The ptraj or cppptraj programs can be used for the analysis and processing of coordinate trajectories. For more information, see the software Users' Manual.

Results

The present protocol can be readily applied in studies of post-translational modification of caspases or pathogenic mutations. In this section, the MD modeling workflow is illustrated (Figure 1), which has been successfully used in the study of caspase-2⁷. Using in vitro site-directed mutagenesis of potential phosphorylation sites (Ser/Thr to Ala) and biochemical approaches, it was demonstrated that the Ser384Ala mutation prevented caspase-2 processing and bl...

Discussion

The described MD approach allows for modeling both the wild-type and mutant forms of caspase using the biomolecular simulation packages. Several important issues of the methodology are discussed here. First, a representative crystal structure of caspase needs to be selected from the Protein Data Bank. Importantly, both monomeric and dimeric forms of caspase are acceptable. Choosing high-resolution structures with a minimum number of missing residues is a good idea. The protonation state of some residues can be set manual...

Materials

Name	Company	Catalog Number	Comments
Amber20	University of California, San Francisco		Software for molecular dynamics simulation http://ambermd.org
AmberTools21	University of California, San Francisco		Software for molecular modeling and analysis http://ambermd.org