Journal of Computational Molecular and Cell Biology Subject Code: IMBI802
Journal of Computational Molecular and Cell Biology Subject Code: IMBI802
Journal of Computational Molecular and Cell Biology Subject Code: IMBI802
BY
Miss. Kasturi Sonawane
Enrolment No: MITU17IMBI0012
Year: 4th year Integrated M. Tech
SUBMITTED TO
2
Practical 01
The Van Der Pol oscillator is an oscillator with nonlinear damping governed by the second-order
differential equation. It is particularly a simple differential equation that can produce a stable limit
cycle. The damping term dictates the dynamics of the oscillator. The equation reduces to the
harmonic oscillator when ǫ = 0. For positive values of ǫ, the system is non conservative, and a limit
cycle exists. The existence of the limit cycle is verified by noticing that when |y| > 1, the damping
is positive (the oscillator dissipates energy) and the amplitude of y decreases. When |y| < 1, the
damping is negative (the system receives energy), and the amplitude of y increases. The limit cycle
is the trajectory for which the average energy balance is null.
First Order Equation: A first order differential equation is an equation of the form F(t,y,˙y)=0. A
solution of a first order differential equation is a function f(t) that makes F(t,f(t),f′(t))=0 for every
value of t.
Lorenz Equation: The Lorenz system is a system of ordinary differential equations first studied by
Edward Lorenz. It is notable for having chaotic solutions for certain parameter values and initial
conditions. In particular, the Lorenz attractor is a set of chaotic solutions of the Lorenz system. In
popular media the 'butterfly effect' stems from the real-world implications of the Lorenz attractor,
i.e. that in any physical system, in the absence of perfect knowledge of the initial conditions (even
the minuscule disturbance of the air due to a butterfly flapping its wings), our ability to predict its
future course will always fail. This underscores that physical systems can be completely
deterministic and yet still be inherently unpredictable even in the absence of quantum effects. The
shape of the Lorenz attractor itself, when plotted graphically, may also be seen to resemble a
butterfly.
3
Procedure (second order differential equation of Van der pol oscillator):
1. Open OpenCOR.
2. From file select new CellML file.
3. Type the code for Van der pol oscillator ode.
4. Save the file.
5. Go to simulation tab on right hand side of the panel.
6. Add three graphs by clicking on the ‘+’ sign on top.
7. Change the ending point to 100 and point interval to 0.1.
8. Select the first graph; from the parameter column right click on x to plot against the
variable of integration.
9. Select the second graph and execute step 8, instead plot y against variable of integration.
10. Select the third graph and plot y against x by right clicking on y from parameter column
and select main to plot against x.
11. Run the simulation by clicking on the run option from the bar on the panel.
12. Analyze the graphs generated through simulation.
1. Open OpenCOR.
2. From file select new CellMl file.
3. Type the code for First order ode.
4. Save the file.
5. Go to simulation tab on right hand side of the panel.
6. Change the ending point to 10 and point interval to 0.5.
7. Select the graph.
8. In the parameter list right click on y to plot against variable of integration.
9. Observe the graph.
10. Exchange the values of y and b.
11. Refresh the graph and plot y against variable of integration.
12. Analyze the graphs generated through simulation.
1. Open OpenCOR.
2. From file select new CellMl file.
3. Type the code for Lorenz attractor ode.
4. Save the file.
5. Go to simulation tab on right hand side of the panel.
6. Add three graphs by clicking on the ‘+’ sign on top.
7. Change the ending point to 50 and point interval to 0.001.
8. Select the first graph; from the parameter column right click on x to plot against the
variable of integration.
4
9. Select the second graph; from the parameter column right click on x, select main to plot
y against x.
10. Select the third graph and plot z against x by right clicking on z from parameter column
and select main to plot against x.
11. Run the simulation by clicking on the run option from the bar on the panel.
12. Observe the graphs generated.
Results:
Fig 2: graphs generated after simulation of Van der Pol Oscillator equation
5
Fig 3: first order differential equation code
6
Fig 5: Graph observed showed exponential growth after exchanging y and b values in first
ODE.
7
Fig 7: graph observed after simulating the Lorenz equation.
Conclusion: A simple second order differential equation of Van der pol oscillator was simulated,
a first order differential equation was simulated to show exponential decay and growth and lastly,
Lorenz attractor ODE was simulated.
8
Practical No. 02: Simulation of neuronal spiking models and bifurcation studies.
Introduction: A model is presented that reproduces spiking and bursting behaviour of known
types of cortical neurons. The model combines the biologically plausibility of Hodgkin-Huxley-
type dynamics and the computational efficiency of integrate-and-fire neurons. Using this model,
one can simulate tens of thousands of spiking cortical neurons in real time (1 mS resolution). The
model for a single neuron must be: 1) computationally simple, yet 2) capable of producing rich
firing patterns exhibited by real biological neurons. Bifurcation methodologies enable to reduce
many biophysically accurate Hodgkin–Huxley-type neuronal models to a two-dimensional (2-D)
system of ordinary differential equations of the form v ‘ = 0.04 v2 + 5v + 140 – u + I (equation
1) and u ‘ = a(bv – u) ---(equation 2) with the auxiliary after-spike setting if v ≥ 30 mV, then v
c, u u + d (equation 3). Here, v and u are dimensionless variables, and a, b, c, and d are
dimensionless parameters, and ‘= d/dt, where t is the time. The variable v represents the
membrane potential of the neuron and u represents a membrane recovery variable, which
accounts for the activation of K+ ionic currents and inactivation of Na+ ionic currents, and it
provides negative feedback to v. After the spike reaches its apex (+30 mV), the membrane
voltage and the recovery variable are reset according to the (equation 3). Synaptic currents or
injected dc-currents are delivered via the variable I. The part 0.04 v2 + 5v + 140 was obtained by
fitting the spike initiation dynamics of a cortical neuron so that the membrane potential v has mV
scale, and the time t has mS scale. The resting potential in the model is between -70 and -60 mV
depending on the value of b. As most real neurons, the model does not have a fixed threshold;
depending on the history of the membrane potential prior to the spike, the threshold potential can
be as low as -55 mV or as high as -40 mV.
The parameter ‘a’ describes the time scale of the recovery variable u. Smaller values result in
slower recovery. A typical value is a = 0.02.
The parameter ‘b’ describes the sensitivity of the recovery variable u to the subthreshold
fluctuations of the membrane potential v. Greater values couple v and u more strongly resulting
in possible subthreshold oscillations and low-threshold spiking dynamics. A typical value is b =
0.2.
The parameter ‘c’ describes the after-spike reset value of the membrane potential v caused by the
fast high-threshold K+ conductances. A typical value is c = =65 mV.
The parameter ‘d’ describes after-spike reset of the recovery variable u caused by slow high-
threshold Na+ and K+ conductances. A typical value is d = 2.
Various choices of the parameters result in various intrinsic firing patterns, including those
exhibited by the known types of neocortical and thalamic neurons. The model is canonical in the
sense that the difference between it and a whole class of biophysically detailed and accurate
9
Hodgkin–Huxley-type models, including those consisting of enormous number of equations and
considering all possible information about ionic currents, is just a matter of coordinate change.
This model was simulated in the CellDesigner software.
Procedure:
4. For simulation, click on ‘Simulation’ on the menu bar and click on ‘Control Panel’.
6. Now, check for the parameters once. Keep ‘Parameter Scan’ and ‘Interactive Simulation’ as it
is.
Results:
Fig. 2 a.) v spike reaching its apex (-30 mV) from initial
point -70 mV.
Conclusion: Hence, the spiking neuron model was successfully simulated in the CellDesigner
software.
11
Practical 03: Simulation and study of Gene regulatory model to understand different
regulatory mechanisms.
Aim: To understand the mechanism involved in a gene regulatory model using in silico tool
TinkerCell.
Introduction: Gene regulatory networks involve a collection of various molecules that interact
with each other and other substances in the cell to control the expression levels of genes and
consequently the mRNA and protein expression levels, ultimately determining a certain function
of the cell.
Lac operon was the first example of a prokaryotic gene regulation studied and understood and has
since been studied extensively and used widely to outline the basic characteristics of a gene
regulatory model.
Lac operon is used by E. Coli to control the transport and metabolism of lactose in situations
where glucose is not available, and lactose needs to be utilised as an energy source. Thus, the lac
operon genes are expressed only when lactose is present and glucose is absent in the cells and is
regulated by two regulators turning the operon on and off, the lac repressor and the CAP
activator, respectively.
Lac repressor here acts as a sensor and blocks transcription of the operon when no lactose is
present but stops inhibiting when lactose is present. This is due to some molecules of lactose
changing to the isomer allolactose that bids to the lac repressor. As a result, the lac repressor can
no longer bind to the DNA at the operator region and so transcription of the operon is carried out
and the lac genes get expressed. Consequently, the lacZ gene is expressed which breaks down
lactose into smaller molecules.
TinkerCell is a computational design software tool used in Synthetic biology that uses a visual
interface combined with an API of Python allowing users to share their code with each other via a
central repository. TinkerCell can be used to draw biological systems with different components
and simulate the system while reiterating the design step multiple times after analysis of the
system so s to get the system to perform a certain task. In this practical we will be using Tinker
Cell to create a bistable lac operon. It is called bistable because the final value of lactose and lacZ
depends on the amount of lactose present inside the cell in the beginning. So, if there is
insufficient lactose in the beginning, the amount of lactose stays low while if there is enough in
the beginning, then this lactose will upregulate its own uptake.
13
Fig 2. Concentration time graph when concentration of lactose is 0.
14
Fig 4. graph showing effect on model due to varying promoter strength.
- Search for glycolysis model in EMBL and download the particular file.
- For simulation, click on ‘Simulation’ on the menu bar and click on ‘Control Panel’.
- Now, check for the parameters once. Keep ‘Parameter Scan’ and ‘Interactive Simulation’ as it
is.
- You can try and change the input concentration values from 5-10 to see whether there is any
change observed in simulation.
15
Fig 1. Imported model of lac operon from EMBL.
Conclusion: Thus, TinkerCell was used to design the Lac Repressor system and the initial
Lactose concentration was varied to see the effect on LacI and LacZ concentrations in the system
with time. Further, promoter strength was altered to see the effect on the model built and, the
model was successfully studied in the CellDesigner software too.
18
Practical No 04: Introduction to CellDesigner.
Aim: To understand the CellDesigner software basics by model building and connecting to
databases, importing the models and performing simulations.
I] Importing model.
- In the menu bar, click on ‘Database’. Click on ‘Import models from BioModels.net.
- Click on Import.
19
Fig 1. The imported model of glycolysis from BioModels.net
20
- Unselecting all the species involved in the reaction. Checking only ‘Glucose’.
- Observe the peak.
- Drag and place the component in the draw area and name it as ‘A’.
- From reaction toolbar, select ‘Catalysis’. Drag it to the state transition reaction between A and
B.
- Add ‘E’ component. Select ‘Inhibition’ from the reaction toolbar and drag it to the State
Transition reaction between B and D.
21
Fig 4. A simple reaction network
Here, protein A undergoes State Transition to protein B under the catalysis of protein C. Protein
B in turn undergoes State transition to protein D and the transition is inhibited by protein E.
- Search for glycolysis model in EMBL and download the particular file.
- For simulation, click on ‘Simulation’ on the menu bar and click on ‘Control Panel’.
- Now, check for the parameters once. Keep ‘Parameter Scan’ and ‘Interactive Simulation’ as it
is.
-You can try and change the input parameter values to see whether there is any change observed
in simulation.
22
Fig 5. The model imported from EMBL database.
23
Fig 7. Graph highlighting important components of the pathway (E4P, GAP, NADH)
Conclusion: Hence, the CellDesigner software was successfully studied using the Glycolytic
pathway and analysing the graphs.
24