Project1-Q1-Handout
Read the paper, and prove (show your steps) that the model described by equations (2) indeed reduce to the model defined by equations (3) for the initial condition: $x_2(0) = x_6(0)$ and $\dot{x}_2(0) = \dot{x}_6(0)$.
Simulate the model defined by the set of equations (3) and (4). Use $\alpha_1=-0.024, \alpha_2=0.0216, \alpha_3=-0.0012, \alpha_4=0.12, C=1.35, \beta=4, H=3$ and initial conditions $x_1=0.0, x_2=0.0, x_3=0.1, x_4=0.0$, along with a timestep of $dt=10^{-4}s$.
Compute the BPM from your simulation to verify that the formula for the multiplicative constant $\Gamma_t$ given by (5) is valid.
Choose $\Gamma_t$ so that you get 100 BPM from your model. Use Winfree’s method to estimate the Phase Response Curve of the ECG signal with respect to (i) $x_2$ (SA node), and (ii) $x_4$ (AV node).
[Bonus]: Compare your result from Q.4(i) against the SA node plot given in Figure 4(a) of this paper. Note that you will have to report the PRC in a slightly different way to make this comparison (refer to the paper for details).
[Further reading]: We can also incorporate breathing and how it couples with the heart’s rhythm into an ECG model, and non-invasively estimate the joint phase response manifold and phase-coupling functions for this model. Details can be found in this paper.
[Your answer for Task 1]¶
# import libraries
import numpy as np
import matplotlib.pyplot as plt
from collections import OrderedDict
from compneuro.base_model import BaseModel
from compneuro.utils.signal import spike_detect
# define ECG model
## you can inherit from BaseModel if you wish to use CompNeuro to simulate your model
# instantiate model and simulation parameters
dt = 1e-4 # set time step
t = np.arange(0, 10, dt) # set time course
## Your model here
## Simulate your model
## Plot resultant ECG signal
## Plot estimated BPM for various Gamma_t
## On the same plot, compare BPM based on the formula (5) to verify its validity
## Use (5) to choose Gamma_t for 100 BPM for the rest of the assignment
## Compute and plot the PRC of the ECG signal with respect to x_2 (SA node)
winfree_steps = 100
perturbation_amp = -0.1
## Compute and plot the PRC of the ECG signal with respect to x_4 (AV node)
winfree_steps = 100
perturbation_amp = 0.1