Computer Simulation of Biological Processes

Project Leader
André H. Juffer, Ph.D., Biocomputing Coordinator
Biocenter Oulu and Faculty of Biochemistry and Molecular Medicine, University of Oulu

Background and Significance

Computer simulation techniques have become important tools to study the properties of biomolecules and their interactions with other molecules. These techniques are employed with the general objective of predicting and explaining physical quantities and phenomena, such as protein-ligand associations, reaction rate constants of enzyme-catalyzed reactions, acid dissociation constants of proteins, protein domain motion, protein folding and stability, and the like. The underlying computational model provides direct insight into the workings of these molecules, such that statements about their functional properties can be offered. Such knowledge is frequently relied upon, for instance, in the drug design industry to speed up the process of drug finding and drug target analysis.

Computer simulation and modeling are not limited to the study of proteins, but can be and are being applied to many different domains of science. This project is specifically concerned with the study of biological processes. Whereas bioinformatics is largely concerned with the generation and analysis of biological data, computer simulation and modeling are applied with the objective of explaining the behavior of a given biological process. Each simulation model typically deals with different time- and length scales. Chemical reactions at the active site of an enzyme occur in very short timescales (fs), protein dynamics generally take place in the range of ps to ns or even longer timescales, interfacial processes (e.g. protein diffusion at membrane surfaces) cover the ms to ms timescale, protein translation falls in the ms to s range, and the growth of tumors occurs over months to years. Each model must therefore involve different simulation and modeling methods.

Simulation and modeling methods rely on the application of techniques that are based in physics (theoretical biophysics), mathematics and chemistry. This project is predominantly concerned with the development and application of such techniques. These range from very detailed QM/MM (quantum mechanics/molecular mechanics) and dynamic models of proteins to multiple-scale approaches that combine, for instance, diffusion models with cellular automata (CA) in the case of brain tumor simulation.

Recent Progress

The antifungal protein chitosanase/glucanase, recently isolated from the soil bacterium Paenibacillus sp. IK-5, displays two CBM32 chitosan-binding modules commonly referred to as DD1 and DD2. They are linked in tandem at the C-terminus1. In order to identify the strength of interaction between ligands and these modules, a series of ligand dockings and free energy calculations (affinity predictions) have been completed. A chitosan trimer (GCS) served as a test ligand. The affinity calculations are based on an application of the so-called double-decoupling method in which the ligand is removed (or decoupled) from a given system, meaning that all interactions between the ligand and all other atoms in the system are slowly switched off in the course of multiple simulations. The system is either the ligand free in solution or the ligand in complex with a binding module (either DD1 or DD2). A typical outcome of a ligand-docking calculation is provided in Figure 1, which shows the complex immediately after docking (left) and after a 50 ns follow-up molecular dynamics simulation (right). Figure 2 displays the change in the free energy as a function of a decoupling parameter. (In collaboration with Tamo Fukamizo.1)

Figure 1: Left: result of a ligand docking with DD1. The docking was completed with AutoDock Vina. Right: the same complex after a 50 ns molecular dynamics simulation. Distances in Å between chemical groups in the ligand and residues in DD1 known to be involved in binding (from NMR data1) are highlighted. These results indicate that the ligand remains bound to DD1 over the course of the simulation, but it tends to drift away from its original binding mode in the docked complex (left). For free energy predictions, it is essential to accurately sample this motion to arrive at a reasonable estimate of the affinity.

Figure 2: An example of a decoupling “reaction” during which the ligand (GCS) is decoupled from the rest of the system. The y-axis is the Gibbs (G) free energy change in units of kT associated with a change in the coupling parameter. The x-axis is an index identifying the current value of the coupling parameter λ, which itself ranges from 0 to 1. For i=20 (corresponding to λ=1), the ligand is fully decoupled from all other atoms in the system. Grey: ligand free in solution; red: ligand in complex with DD1; blue: ligand in complex with DD2. Similar profiles can be found for the reverse reaction.

1 Shinya et al, Biochem. J. (2016) 473 ,1085–1095.

Water networks in thiolase
This subproject is concerned with the importance of hydrogen bond networks (HBNs) for the efficient catalysis of the Zoogloea ramigera thiolase active site. This is investigated by in silico modeling and simulation techniques (molecular dynamics in particular) and structural enzymological approaches. These studies will benefit very much from the presence of a range of key structures of this thiolase, complexed with various substrates, substrate analogs and intermediates. Two HBNs will be studied, being related to the properties of oxyanion hole 1 (OAH1, stabilizing an enolate intermediate) and oxyanion hole 2 (OAH2, stabilizing a tetrahedral intermediate) of thiolase. The hydrogen bond donors (HBDs) of these oxyanion holes are from main chain, side chain or water molecules. In the HBN of OAH1 water molecules play a key role, whereas in OAH2 a peptide-mediated HBN is involved. This study capitalizes on the enormous amount of structural data already available on the thiolase system. In addition, the computational studies will be complemented and validated by structural enzymological studies (collaboration with Prof. Wierenga’s group, See also Figure 3).

Figure 3: Evolution of the amount of buried water (black line) and the length of the water trail (red line) as a function of time. The top left inset displays the catalytic site of the thiolase, with S-acetylated cysteine and bound acetyl-CoA (green) in the oxyanion hole OAH2 formed by an histidine and a water molecule (red). The latter is part of a water trail.The top right inset displays an example of a buried water analysis. The alpha shape surface (grey) defines an outer surface of the protein. Inside the surface, buried water molecules (blue cross) and the water trail (red) are detected by an in-house developed algorithm that analyses the trajectory generated from molecular dynamics simulations. The properties of the water trails may be affected by mutations altering the catalytic efficiency of thiolase. For this particular analysis, a whole set of new software tools were developed.

Proton transfer events in coarse-grained molecular dynamics
Any theoretical molecular model for describing biological processes involving proteins should account for possible variations in the protein's charge distribution due to proton (dis)association by ionizable groups in the protein. In this project we propose a novel model for charge fluctuations in coarse-grained molecular dynamics (CGMD), a simulation methodology that has become one of the most important techniques for the study of proteins and membranes at large length and long time scales. Conventional CGMD is based on Newtonian mechanics. Like standard atomistic molecular dynamics (MD), it computes forces on CG particles or “beads” to generate a sequence of states (particle coordinates and velocities) as a function of time, from which numerous molecular, thermodynamic and time-dependent properties of the system can be calculated. However, charge  (and associated mass) fluctuations in CGMD are completely ignored. This subproject is aimed at correcting this serious omission (See also Figure 4).

Figure 4: Simple simulation of a stochastic current I(t) (in blue) between two coarse-grained  particles. The x-axis represents time t. The time scale is not representative of the proton transfer process. The current is a consequence of a proton transfer or jump event, the occurrence of which is controlled by the state of the system, and results in a new current or a change of direction of still ongoing current (highlighted by vertical lines at the jump time). The current affects the charges Q(t) of the particles. Note that the change in Q(t) is continues so that force calculations can be carried out as in standard CGMD.

Future Goals

  1. Development of a novel method for the inclusion of proton binding in coarse-grained simulation of proteins.
  2. Deriving scientifically sound and probabilistic optimal principles for coarse-grain molecular systems for simulating protein adsorption at interfaces.

Doctoral Theses 2017

Leonardo Garma, Structural bioinformatics tools for the comparison and classification of protein interactions, Acta Univ. Oul. D 1422, 2017 (

Research Group Members

Project Leader:
André H. Juffer, Ph.D. (Biocenter Oulu)

Ph.D. Students:
Pierre Leprovost, MSc (Biocenter Oulu)
Ioannis Beis, M.Sc.
Outi Lampela, M.Sc.
Leonardo Garma, M.Sc.

Florian Borse, M.Sc. (Biocenter Oulu, 3 months)

Undergraduate Students:
Ella Rukajärvi

Foreign Scientists, 4

Last updated: 5.7.2018