Script for Enzyme QM/MM Studies

The structure file is now built. Firstly, the substrate is generated and designated as segment NEU5. The crystallographic determined water molecules are then set up and identified as WATR; the calcium ion bound to the protein is then generated with identifier CION. Finally, the co-ordinate file for the ES complex is then read.

 

READ SEQUENCE CARD                                                                UNIT 16

GENERATE neu5 SETUP

 

READ SEQUENCE CARD                                                                UNIT 15

GENERATE prot SETUP

 

READ SEQUENCE TIP3 116

GENERATE watr SETUP NOANG NODIHE

 

READ SEQUENCE CARD                                                                UNIT 17

GENERATE cion SETUP

 

READ COOR PDB                                                                              UNIT 13

 

Because there are disulfide links in the NA structure, CHARMM must be instructed to delete the ─SH hydrogen atoms and form the S─S links between the cysteine residues.

 

PATCH DISU prot 94  prot 112 SETUP

PATCH DISU prot 197 prot 210 SETUP

PATCH DISU prot 199 prot 208 SETUP

 

A 17Ĺ water sphere is then added to the system; the 694 added solvent molecules are designated SOLV.

 

READ SEQUENCE TIP3 694

GENERATE solv SETUP NOANG NODIHE

   READ COOR PDB APPEND                                                         UNIT 14

The QM link atom is then added to the Asp-151 side chain. To incorporate it into the “start-up” system, it is placed half way between the Ca─Cb bond axis. The QM link atom is known as HQA.

 

COOR AXIS SELECT ( (SEGID prot .AND. RESID 70 ) .AND. TYPE CA) END -

SELECT ((SEGID prot .AND. RESID 70 ) .AND. TYPE CB) END

 

COOR SET xdir ?xcen ydir ?ycen zdir ?zcen SELECT (SEGID prot .AND. –

RESID 70 .AND. TYPE HQA) END

 

The next stage is to remove unwanted atoms from the system. Protein and WATR residues beyond 15Ĺ of the substrate are deleted. A SOLV molecule whose oxygen is within 2.5Ĺ of any heavy atom in the ES complex is also deleted.

 

DELETE ATOM SELECT .NOT. (.BYRES. (SEGID neu5 .AROUND. 15.00) .OR. -

(SEGID solv)) END

 

DELETE ATOM SELECT (.BYRES. ( (SEGID solv .AND. TYPE OH2 ) .AND.-

 ( (.NOT. SEGID solv .AND. .NOT. HYDROGEN) .AROUND. 2.5) ) ) END

 

Atom-based non-bonded cut-offs are set to very large (40Ĺ) and a 19Ĺ boundary potential is applied to the SOLV and WATR molecules. For heavy protein and crystallographic water atoms beyond 13Ĺ of the origin, a harmonic force constraint is applied. The force constant is mass-weighted with a value of 1.0.

 

nbonds atom switch cutnb 40.00 ctofnb 38 ctonnb 36

 

 

OPEN NAME  ./wat19.pot UNIT 4 FORM READ

sbound READ UNIT 4

sboundary set xref 0.0 yref 0.0 zref 0.0 -

          assign 1 selection (segid watr .and. type O*) .OR.-

(segid solv .and. type O*) end

 

CONSTRAINT HARMONIC FORCE 1.0 SELECT ((SEGID prot .OR. SEGID watr) -

 .AND.  (.NOT. TYPE H*) .AND. (.NOT. POINT 0.0 0.0 0.0 CUT 13.0))-

 END  MASS EXPONENT=2

 

After setting up the system, it is now possible to restart a simulation if required. This is done by reading the restart file into the comparison co-ordinate array and swapping it with the co-ordinates of the current system. This ensures that all the atoms in the restart file are set up as in the original simulation. The QM link atom in the restart file will have already been added; initially such an atom is placed 1.16Ĺ from the Cb atom.

 

 open name fluptr1s7mk2pt3.pdb unit 77 form read

 read coor comp pdb unit 77

 close unit 77

 coor swap

 

It is often useful to define a sub-set of atoms. The QM link atom is included in the QuantumASP definition.

 

 

DEFINE QuantumASP SELECT -

(RESID 70 .AND. (TYPE CB .OR. TYPE HB1 .OR. TYPE HB2 .OR. TYPE CG -

.OR. TYPE OD2 .OR. TYPE OD1 .OR. TYPE HD2 .OR. TYPE HQA) ) END

 

DEFINE MolMechSolv SELECT ( (SEGID solv) .OR. (SEGID watr .AND. -.NOT. (RESID 146 .OR. RESID 148) ) ) END

 

DEFINE pnp SELECT ( (SEGID neu5) .AND. (TYPE C19 .OR. TYPE C20.OR. -

TYPE C21 .OR. TYPE C22 .OR. TYPE N23 .OR. TYPE C24 .OR. TYPE C25 -.OR. TYPE H45 .OR. TYPE H46 .OR. TYPE O47 .OR. TYPE O48 .OR. TYPE-

 H49 .OR. TYPE H50) ) END

 

The initial adiabatic mapping process can be achieved in at least two ways. The first way is to define a vector for a transferring atom. This was done in the enzyme study when investigating the proton transfer and attack of water to the zwitterion. The translation vector for the proton transfer was calculated from the proton’s position in the reactant and its position in the intermediate. The vector is then divided into roughly ten sections (in Cartesian co-ordinates) and the atom is translated by this amount to produce the energy profile. The Quick command can give the distance between two atoms (in this case 38 and 40) and thus serves as a check to see if the translation has correctly occurred. The second approach was used in the a-lactone/zwitterion solvation study. This involves extending the bond to a specific value and then applying an internal co-ordinate constraint. This approach is included here but is commented out.

 

SET a -0.06

SET b 0.065

SET c -0.0783

COOR TRANS SELECT (RESID 70 .AND. TYPE HD2) END xdir @a ydir @b –

zdir @c

 

QUICK 38 40

 

!set dist1 ?DIST

!incr dist1 by 0.05

!IC DELETE SELECT .NOT. (BYNUM 38 .OR. BYNUM 40) END

!IC EDIT

!BOND BYNU 38 BYNU 40 @dist1

!END

!CONS IC BOND 1e7 EXPO 2

 

The substrate and Asp(H)-151 side chain are treated with the AM1 Hamiltonian. Because the Asp residue is protonated, the overall charge of the QM system is –1 (associated with the carboxylate group of the substrate).

 

QUANTUM SELECT (SEGID neu5 .OR. QuantumASP) END AM1 CHARGE -1

 

The script up to this point is defined as the CHARMM QM/MM system set up.

 

Rigid constraints are now applied. For the enzyme system, this will include the QM link atom. To generate the reaction profile, the proton on the Asp(H)-151 carboxylate group and the corresponding carboxylate oxygen are also fixed.

 

CONS FIX SELECT (TYPE HQA .OR. (RESID 70 .AND. (TYPE OD2 .OR. TYPE –

 HD2))) END

 

A file called arguments is generated containing lines within a title (such lines start with a “*”). These report the type of optimiser used, the gradient tolerance (in kcal mol-1 Ĺ-1), the maximum allowed optimisation cycles, the restart file name, the number of optimisation cycles already done and the number of cycles allowed before each that restart file is generated.

 

OPEN NAME arguments                                            UNIT  77 FORM WRITE

WRITE TITLE UNIT 77

* set 1 abnr                   ! select ABNR minimiser

* set 2 0.005                 ! select GRMS exit criterion

* set 5 600                    ! select the maximum number of cycles

* set 6 fluptr2s8.res                  ! select the name of the restart file

* set 7 0                                                 ! number of cycles already done

* set 8 100                                            ! Set the required number of cycles

* return

*

CLOSE UNIT 77

 

The variable,  r, is used to control the loops that occur later in the script. The goto command goes to the minimisation part of the script.

 

set r min_home_2

goto MINIMISE

label min_home_2

 

The command, STOP, is used to finish the simulation and to terminate the running of CHARMM.

 

STOP

 

 

The label sets the position of the script recognised by the goto command; “goto MINIMISE” will thus initiate a minimisation of the system. The first stage of minimisation using this script is to read in the arguments.

 

label MINIMISE

 

OPEN NAME arguments                                UNIT 77 FORM READ

STREAM UNIT 77

 

Minimise the system using the ABNR optimiser for 600 steps with an initial step size of 0.02Ĺ; system energies are printed after 20 cycles. The 600 steps are fully carried out before the total-number-of-steps-done counter is then increased by 600. The restart file title gives information on the total energy and RMS gradient of the system, and the number of ABNR cycles already done. It also includes the current length of the bond being elongated. The title may also contain a history of previous optimisations done to arrive at the restart structure. Unit 77 is then closed off for future use.

 

label minimise_loop

 

mini @1 nstep @8 nprint 20 step 0.02 tolgrd 0.0

incr 7 by @8

 

 

OPEN NAME @6                                                                   UNIT 77 FORM WRITE

WRITE COOR PDB UNIT 77

* PREVIOUS OPTIMISATIONS:

*   As For SIA and QM/MM Quantised (flusi0qmmm) and

*   QM/MM 800 cycles ABNR SIP QM ELSE MM; H20 FIXED

*   QM/MM 3000 cycles ABNR SIP/ASP151 QM ELSE MM; H20 FIXED

* THIS OPTIMISATION:

* @7 cycles @1 

*       ASP-151 H20 DISTANCE IS @9

* created by STUART FIRTH-CLARK, UNIVERSITY OF BATH

*

CLOSE UNIT 77

 

Now the gradient tolerance is tested. If it is below the specified value (0.005 kcal mol-1 Ĺ-1) then the script moves forward to minimiser_gohome which sends the program back to the main label (just before the STOP command) in the script. However, if the gradient is too high, the program is sent back to the minimise_loop label for a further 100 optimisations cycles. Applying the gradient tolerance test in this way (rather than specifying it when minimising the system) reduces the chances of premature terminations in the minimisation procedure. “Premature terminations” can sometimes occur because the RMS gradient of the (large) system “bounces”; such an effect can sometimes drop the RMS gradient below the tolerance before the system is fully minimised. The script also moves forward if the number of optimisation cycles has reached the maximum value (given by variable 5)

 

if 7 eq @5 goto minimiser_gohome

if 2 gt ?GRMS goto minimiser_gohome

goto minimise_loop

 

label minimiser_gohome

 

goto @r