Circadian Clock ODE Model (MATLAB)

A MATLAB re-implementation of a cell-autonomous mammalian circadian clock, using 21 ODEs and 132 parameters to reproduce key figures and knockout behaviors described in the literature. The code simulates paired mRNA/protein dynamics (PER, CRY, REV-ERBα, CLK, BMAL1, RORC) and exports publication quality figures.

This project began as a University of Maryland College Park BIOE 241 final project and has since been lightly refactored for clarity, reproducibility, and modern MATLAB (R2024b). The original modeling goals and figures are summarized in the included report.

---

Brief Background

The mammalian circadian clock is a biochemical oscillator that regulates sleep, metabolism, hormone release, and temperature Disruptions are linked to conditions like insomnia, myocardial infarction, obesity, diabetes, and immune dysfunction.

Core clock genes (Per1/2/3, Cry1/2) are activated by the CLK/BMAL1 complex and form delayed negative feedback loops, producing ~24-hour rhythms at the cellular level; the master pacemaker resides in the SCN (~20k neurons) with peripheral clocks elsewhere.

In the class project, we reproduced two key figures from the literature: (i) eight paired mRNA/protein traces (Rev-erbα, Rorc, Per1/2, Cry1/2, Bmal1, Clk), and (ii) a Per2 comparison in wild-type vs Cry2−/− (knockout shows larger amplitude/period). The model implements the published 21 ODEs with 132 parameters to capture these oscillations.

Contents

circadian_clock/
├─ run_simulation.m         % runs ODEs + saves figures
├─ ode_circadian.m          % RHS of the 21-dimensional ODE system
├─ model_params.m           % 132 model parameters as a struct (p.*)
├─ plot_timeseries.m        % 8-panel mRNA/protein grid (Figure 2-style)
├─ figures/                 % auto-created; exported .pdf/.png/.fig written here
└─ report/
   └─ BIOE241_Final_Report.pdf  % project report (context + objectives)

---

What it does

• Implements the 21 ODEs from a standard mammalian circadian clock model.

• Uses Michaelis–Menten for transcription and mass action elsewhere.

• Enforces nonnegativity on all 21 states during integration.

• Reproduces:

  • Figure 2-style 8 subplot grid (mRNA in blue, protein in red; paired axes).
  • Figure 3-style comparison of Per2 mRNA (wild type vs Cry2−/−) as a separate figure.

• Saves vector PDF/SVG (optional), PNG, and MATLAB .fig to a figures/ folder.

---

Requirements

• MATLAB R2024b or newer recommended (orginally tested with R2019a+).
• No toolboxes beyond base MATLAB.
• OS: Windows/macOS/Linux.

---

Quick start

1. Open the folder in MATLAB.

2. Edit the output path in run_simulation.m:

   % Absolute path (so exports don’t depend on current folder)
   output_directory = 'D:\User Data\Downloads\circadian_clock\figures';

3. Run:

4. Figures are written into figures/ with a timestamped prefix:

   • ..._timeseries.pdf/.png/.fig (8-panel grid)
   • ..._per2_compare.pdf/.png/.fig (Per2 vs Cry2−/−)

---

How it’s organized

model_params.m

All parameter scalars are placed on a struct p (ex: p.v0per1, p.KA1cry1, ...).

ode_circadian.m

Right-hand side of the ODEs. Signature:

function dydt = ode_circadian(t, y, p)

• y is a 21×1 column vector with this order:

  [per1; per2; cry1; cry2; reverba; clk; bmal1; rorc; PER1; PER2; CRY1; CRY2; REVERBA; CLK; BMAL1; RORC; PER1CRY1; PER2CRY1; PER1CRY2; PER2CRY2; CLKBMAL1]
  

• All parameter references are p. names.

run_simulation.m

• Defines initial conditions and tspan.

• Sets solver options with nonnegativity constraints:

  opts = odeset('NonNegative',1:21,'RelTol',1e-6,'AbsTol',1e-9);

• Calls:

  rhs = @(t,y) ode_circadian(t,y,p);
  [t,y] = ode45(rhs, [0 50], y0, opts);

• Plots:

  • 8-panel grid via plot_timeseries.
  • Separate Per2 vs Cry2−/− line plot.

• Exports figures via exportgraphics and savefig.

plot_timeseries.m

Builds a tiled layout (4×2) and reproduces the dual-axis styling: mRNA on left (blue), protein on right (red), with axis colors set accordingly.

---

Reproducing the figures

• 8-panel (Figure 2-style) uses:

  • Initial mRNA: per1=.25, per2=.25, cry1=.1, cry2=1, reverba=.05, clk=1.15, bmal1=.712, rorc=.1
  • Initial proteins: PER1=.78, PER2=12, CRY1=5, CRY2=4.5, REVERBA=.29, CLK=8.8, BMAL1=7.6, RORC=.05
  • Heterodimers & CLKBMAL1 start at 0.

• Per2 vs Cry2−/− (Figure 3-style):

  • The script generates a separate comparison figure after integrating; the dashed line corresponds to the Cry2 knockout condition as in the class report.

See the report in report/ for background, objectives, and figure descriptions.

---

Output

• figures/circadian_YYYYMMDD_HHMMSS_timeseries.pdf|png|fig
• figures/circadian_YYYYMMDD_HHMMSS_per2_compare.pdf|png|fig

---

Citation

If you use or extend this code, please cite the class report and the original model paper referenced within it:

• BIOE241 Final Project Report (included in this repository).
• Mirskya, H.P., Liu, A.C., Welsh, D.K., Kay, S.A., and Doyle, F.J. 2009. A model of the cell-autonomous mammalian circadian clock. PNAS Vol. 106, No. 27 (2009) 11107-11112.