Location: Computational analysis of the human sinus node action potential @ 7f1518641c18 / Figure6.py

Author:

nima <nafs080@aucklanduni.ac.nz>

Date:

2021-06-22 14:51:01+12:00

Desc:

Merge branch 'master' of https://models.physiomeproject.org/workspace/648

Permanent Source URI:

https://staging.physiomeproject.org/workspace/648/rawfile/7f1518641c18a9570b48e21e74c4836a960d5a07/Figure6.py

# To reproduce Figure 3 in the associated Physiome paper,
# execute this script from the command line:
#
#   cd [PathToThisFile]
#   [PathToOpenCOR]/pythonshell Figure6.py

import matplotlib

matplotlib.use('agg')

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

import opencor as opencor

# different values for gf to decrease the If
K_NaCa = [1.6715, 0.3343]

# load the reference model
simulation = opencor.open_simulation("HumanSAN_Fabbri_Fantini_Wilders_Severi_2017.sedml")
data = simulation.data()
data.set_ending_point(1.8)
data.set_point_interval(0.001)

simulation.reset(True)

results = np.zeros((13, 1801))

for value in range(len(K_NaCa)):
    simulation.reset(True)
    simulation.clear_results()

    data.constants()["i_NaCa/K_NaCa"] = K_NaCa[value]

    for i in range(8):
        simulation.run()
        simulation.clear_results()

    simulation.run()

    ds = simulation.results().data_store()

    results[value] = ds.voi_and_variables()["Membrane/V"].values()

simulation.reset(True)
simulation.clear_results()
data.constants()["i_NaCa/K_NaCa"] = 0.83575

for i in range(3):
    simulation.run()
    simulation.clear_results()

simulation.run()

ds = simulation.results().data_store()

results[3] = ds.voi_and_variables()["Membrane/V"].values()

simulation.reset(True)
simulation.clear_results()
data.constants()["i_NaCa/K_NaCa"] = 3.343

for i in range(42):
    simulation.run()
    simulation.clear_results()

simulation.run()

ds = simulation.results().data_store()

results[2] = ds.voi_and_variables()["Membrane/V"].values()

for i in range(0, 1):
    for i in range(3):
        simulation.run()
        simulation.clear_results()
    simulation.run()

    ds = simulation.results().data_store()
    results[4] = ds.voi_and_variables()["environment/time"].values()

for value in range(len(K_NaCa)):
    simulation.reset(True)
    simulation.clear_results()

    data.constants()["i_NaCa/K_NaCa"] = K_NaCa[value]
    for i in range(8):
        simulation.run()
        simulation.clear_results()
    simulation.run()

    ds = simulation.results().data_store()

    results[value + 5] = ds.voi_and_variables()["Ca_dynamics/Cai"].values()

simulation.reset(True)
simulation.clear_results()
data.constants()["i_NaCa/K_NaCa"] = 0.83575

for i in range(3):
    simulation.run()
    simulation.clear_results()

simulation.run()

ds = simulation.results().data_store()

results[8] = ds.voi_and_variables()["Ca_dynamics/Cai"].values()

simulation.reset(True)
simulation.clear_results()
data.constants()["i_NaCa/K_NaCa"] = 3.343

for i in range(42):
    simulation.run()
    simulation.clear_results()

simulation.run()

ds = simulation.results().data_store()

results[7] = ds.voi_and_variables()["Ca_dynamics/Cai"].values()

for value in range(len(K_NaCa)):
    simulation.reset(True)
    simulation.clear_results()

    data.constants()["i_NaCa/K_NaCa"] = K_NaCa[value]
    for i in range(8):
        simulation.run()
        simulation.clear_results()
    simulation.run()

    ds = simulation.results().data_store()

    results[value + 9] = ds.voi_and_variables()["i_NaCa/i_NaCa"].values()

simulation.reset(True)
simulation.clear_results()
data.constants()["i_NaCa/K_NaCa"] = 0.83575

for i in range(3):
    simulation.run()
    simulation.clear_results()

simulation.run()

ds = simulation.results().data_store()

results[12] = ds.voi_and_variables()["i_NaCa/i_NaCa"].values()

simulation.reset(True)
simulation.clear_results()
data.constants()["i_NaCa/K_NaCa"] = 3.343

for i in range(42):
    simulation.run()
    simulation.clear_results()

simulation.run()

ds = simulation.results().data_store()

results[11] = ds.voi_and_variables()["i_NaCa/i_NaCa"].values()

# define the x and y axis and match the units
X = results[4] * 1000
Y1 = results[2]
Y2 = results[0]
Y3 = results[3]
Y4 = results[1]
Y5 = results[7] * 1e6
Y6 = results[5] * 1e6
Y7 = results[8] * 1e6
Y8 = results[6] * 1e6
Y9 = results[11] * 1000 / 57
Y10 = results[9] * 1000 / 57
Y11 = results[12] * 1000 / 57
Y12 = results[10] * 1000 / 57

plt.figure(figsize=(16, 14))
plt.subplot(2, 2, 1)

plt.plot(X, Y1, 'navy', linestyle='-', label='CTRL', linewidth=3)
plt.plot(X, Y2, 'red', linestyle='-', label='Block 50%', linewidth=3)
plt.plot(X, Y3, 'green', linestyle='-', label='Block 75%', linewidth=3)
plt.plot(X, Y4, 'purple', linestyle='-', label='Block 90%', linewidth=3)

plt.gca().add_patch(Rectangle((440, -65), 695, 30, linewidth=2, edgecolor='black', linestyle='--', facecolor='none'))
plt.grid()
plt.xlim(0, 1800)
plt.ylim(-70, 30)
plt.xticks(np.arange(0, 1800, 500))

plt.tick_params(axis='both', labelsize=14)
plt.ylabel('V$_m$ (mV)', fontsize=16)
plt.title('A', loc='left', y=1.05, x=-0.06, fontsize='20')

plt.subplot(2, 2, 2)
plt.plot(X, Y1, 'navy', linestyle='-', label='CTRL', linewidth=3)
plt.plot(X, Y2, 'red', linestyle='-', label='Block 50%', linewidth=3)
plt.plot(X, Y3, 'green', linestyle='-', label='Block 75%', linewidth=3)
plt.plot(X, Y4, 'purple', linestyle='-', label='Block 90%', linewidth=3)

plt.grid()
plt.xlim(445, 1120)
plt.ylim(-63, -35)

plt.tick_params(axis='both', labelsize=14)
plt.ylabel('V$_m$ (mV)', fontsize=16)
plt.title('B', loc='left', y=1.05, x=-0.06, fontsize='20')

plt.subplot(2, 2, 3)
plt.plot(X, Y5, 'navy', linestyle='-', label='CTRL', linewidth=3)
plt.plot(X, Y6, 'red', linestyle='-', label='Block 50%', linewidth=3)
plt.plot(X, Y7, 'green', linestyle='-', label='Block 75%', linewidth=3)
plt.plot(X, Y8, 'purple', linestyle='-', label='Block 90%', linewidth=3)

plt.grid()
plt.xlim(0, 1800)
plt.xticks(np.arange(0, 1800, 500))
plt.ylim(50, 450)
plt.xlabel('Time (ms)', fontsize=16)
plt.tick_params(axis='both', labelsize=14)
plt.ylabel('[Ca$^{2+}]_i$ (nM)', fontsize=16)
plt.title('C', loc='left', y=1.05, x=-0.06, fontsize='20')

plt.subplot(2, 2, 4)
plt.plot(X, Y9, 'navy', linestyle='-', label='CTRL', linewidth=3)
plt.plot(X, Y10, 'red', linestyle='-', label='Block 50%', linewidth=3)
plt.plot(X, Y11, 'green', linestyle='-', label='Block 75%', linewidth=3)
plt.plot(X, Y12, 'purple', linestyle='-', label='Block 90%', linewidth=3)

plt.grid()
plt.xlim(445, 1130)
plt.ylim(-0.16, 0.05)
plt.yticks(np.arange(-0.15, 0.01, 0.05))
plt.xlabel('Time (ms)', fontsize=16)
plt.tick_params(axis='both', labelsize=14)
plt.ylabel('I$_{NaCa}$ (pA/pF)', fontsize=16)
plt.title('D', loc='left', y=1.05, x=-0.06, fontsize='20')
plt.legend(bbox_to_anchor=(0., 0.85, 1, 0.1), loc='best', fontsize=14,
           ncol=10, mode="expand")

plt.tight_layout(pad=0.5, w_pad=3.5, h_pad=3)

plt.savefig('Figure6.png')