[[model]]
name: Simon-Chica-2021
desc: """
Murine resident cardiac maacrophage model [1].

** Type 2 **


[1] Simon-Chica et al. (2021) Circ Res
"""


# Initial values
iKur.R = 3.3840e-4
iKur.S = 0.9652

cell.V = -39.55 [mV]

[mac]
Cm = 18.32e-6 in [uF]

[Environment]
time = 0 bind time
pace = 0 bind pace


#MACROPHAGE MODEL
# Background current
[ib]
use cell.V as V 
use mac.Cm as Cm
gb = 2.2e-7 in [mS] #mS
Eb = 0 in [mV]
Ib = (gb/Cm * (V - Eb)) in [pA/pF] # Background current


[ik1]
use cell.V as V
F = 96487 [C/mol]
R = 8314 [mJ/K/mol]
T = 275+22 [K]
K_i = 5.4 # is actually mM but should be dimensionless here because of the sqrt below
GKir = 7.3e-7 [mS]
aKir = 0.8
bKir = 0.7 
EK = -79.362 [mV]
IK1=((GKir/mac.Cm * sqrt(K_i) * (V - (EK))) / (aKir + exp((bKir * (V - (EK)) * F)/(R * T)))) in [pA/pF]

[iKur]
use cell.V as V
EK = -79.362 [mV]
gKur = 3e-7 [mS]
a = 86 
b = 18.8
c = 1.2
d = 1
IKrinf = 1/(1 + exp((V + 6)/a))
IKsinf = 1/(1 + exp((V + 7.5)/b)) 
IKrtau = (0.0066/(1 + exp((V + 5)/c)) + 0.00036)
IKstau = (0.43/(1 + exp((V + 60)/d)) + 2.2)
dot(R) = (IKrinf - R)/IKrtau
dot(S) = (IKsinf - S)/IKstau

IKur= gKur/mac.Cm * R * S * (V - EK) in [pA/pF]


[cell]
i_diff = 0 bind diffusion_current 
     
dot(V) = -(iKur.IKur + ib.Ib + ik1.IK1 + (i_diff)/mac.Cm) in [mV] label membrane_potential

[[protocol]]
# Level  Start    Length   Period   Multiplier
1.0      0.0      2        1000.0   10

[[script]]
import matplotlib.pyplot as plt
import myokit
import numpy as np
from mpl_toolkits.mplot3d import axes3d
import math
from scipy.signal import find_peaks, peak_prominences
from scipy import signal


numstim = 0 #5 #50 #for prestimulation 
bcl = 200

# Get model and protocol, create simulation
   
dt = 0.005
sim_time = bcl*5
n_stim = 1
#g_vec= [(0,1,2,8)]

n_end = 1
n_mid = 1
n_epi = 0
n = n_end + n_mid + n_epi
m = get_model()
m.resolve_interdependent_components()
p = myokit.Protocol()
p.schedule(level = 1, start = 0, duration = 2, period=bcl, multiplier=10)
s = myokit.SimulationOpenCL(m, p, ncells=n)
s.set_step_size(dt)
s.set_paced_cells(n_stim)
s.pre(numstim*bcl)
    
#CONNECTIONS
g= 167e-6 #mS/uF
g_NM = 0 #5e-6 #mS/uF
#connections = [(i, i+1, g_NM if n_end-1<= i <=n_end+n_mid-1 else g) for i in range(n-1)]
connections = [(0, i+1, g_NM) for i in range(n-1)]
s.set_connections(connections)
  
#DEFINE CABLE
types = [0] * n_end + [2] * n_mid + [1] * n_epi
s.set_field('cell2.type', types)
    
d = s.run(sim_time, log_interval=0.01, log=['Environment.time', 'remaining.V','ecc.na_m','ecc.na_h','ecc.na_j'],report_nan=False) #maybe sim_time can be smaller for checking amp?
#d.save_csv('mouse_type0_NOinward_5Cm_5nS_1Macs_5g.txt')

numts=(d['remaining.V',0])
vm=np.zeros((len(numts),n))
t = d['Environment.time']   
#voltage vector         
for j in range(0, n):
    vm[:,j]=d['remaining.V', j]

plt.figure()
plt.plot(t,vm[:,-1])
plt.title(' sim_time={}, t={}'.format(sim_time, t[-1]))
plt.plot(t,vm[:,0],'r')