Model is modified from original Moldakarimov model seen in this paper, “Competitive Dynamics in Cortical Responses” by Moldakarimov et al.

The below code should produce an image like this in XPP:

#Moldakarimov model with two excitatory, two inhibitory neurons

#Parameters

par Jee = 0.0000771185, Jii = 0.0008748672024132885, Jei = 1.4090255443147903, Jiei = 0.013533455878780463, Jiec = 0.1269213277585255, AA = 34.259728192954, psi = 1.5371052948343185, fe = 0.006579634656629687, fi = 0.0021210156666630615

par te = 10, ti = 8, tge = 100, tgi = 80

par I1 = 0.8, I2 = 0.2, NC = 2, I1i=1, I2i=1

par cs1 = 9.114343343789908, cs2 = 0.00010877536894466077, cs3 = 0.00149053109891304, cs4 = -0.48447454404577495, cs0 = -0.3316691998368452

par VL = -65, VNa = 55, VK = -80, Vei = -80, Vee = 0, Vie = 0, Vii = -80

par gL = 0.05, gNa = 100, gK = 40, gCa = 0.975

#Equations

dVe1/dt = Iexte1 – (gL * (Ve1 – VL) + gK * ne1^4 * (Ve1 – VK) + gNa * minf(Ve1)^3 * he1 * (Ve1 – VNa)) – (Jee * ge1 * (Ve1 – Vee) / NC + Jei * gi1 * (Ve1 – Vei) / NC)

dVe2/dt = Iexte2 – (gL * (Ve2 – VL) + gK * ne2^4 * (Ve2 – VK) + gNa * minf(Ve2)^3 * he2 * (Ve2 – VNa)) – (Jee * ge2 * (Ve2 – Vee) / NC + Jei * gi2 * (Ve2 – Vei) / NC)

dhe1/dt = psi * (alphah(Ve1) * (1 – he1) – betah(Ve1) * he1)

dhe2/dt = psi * (alphah(Ve2) * (1 – he2) – betah(Ve2) * he2)

dne1/dt = psi * (alphan(Ve1) * (1 – ne1) – betan(Ve1) * ne1)

dne2/dt = psi * (alphan(Ve2) * (1 – ne2) – betan(Ve2) * ne2)

dse1/dt = (AA * sigma(Ve1) * (1 – se1) – se1) / te

dse2/dt = (AA * sigma(Ve2) * (1 – se2) – se2) / te

dphie1/dt = (-fe * sigma(Ve1) * phie1 + (1 – phie1)) / tge

dphie2/dt = (-fe * sigma(Ve2) * phie2 + (1 – phie2)) / tge

ge1 = se1 * phie1

ge2 = se2 * phie2

dVi1/dt = I1i-(gL * (Vi1 – VL) + gK * ni1^4 * (Vi1 – VK) + gNa * minf(Vi1)^3 * hi1 * (Vi1 – VNa)) – (Jii * gi1 * (Vi1 – Vii) / NC + Jiei * ge1 * (Vi1 – Vie) / NC + Jiec * ge2 * (Vi1 – Vie) / NC)

dVi2/dt = I2i-(gL * (Vi2 – VL) + gK * ni2^4 * (Vi2 – VK) + gNa * minf(Vi2)^3 * hi2 * (Vi2 – VNa)) – (Jii * gi2 * (Vi2 – Vii) / NC + Jiei * ge2 * (Vi2 – Vie) / NC + Jiec * ge1 * (Vi2 – Vie) / NC)

dhi1/dt = psi * (alphah(Vi1) * (1 – hi1) – betah(Vi1) * hi1)

dhi2/dt = psi * (alphah(Vi2) * (1 – hi2) – betah(Vi2) * hi2)

dni1/dt = psi * (alphan(Vi1) * (1 – ni1) – betan(Vi1) * ni1)

dni2/dt = psi * (alphan(Vi2) * (1 – ni2) – betan(Vi2) * ni2)

dsi1/dt = (AA * sigma(Vi1) * (1 – si1) – si1) / ti

dsi2/dt = (AA * sigma(Vi2) * (1 – si2) – si2) / ti

dphii1/dt = (-fi * sigma(Vi1) * phii1 + (1 – phii1)) / tgi

dphii2/dt = (-fi * sigma(Vi2) * phii2 + (1 – phii2)) / tgi

gi1 = si1 * phii1

gi2 = si2 * phii2

#Where

alpham(x) = 0.1 * (x + 30) / (1 – exp(-0.1 * (x + 30)))

betam(x) = 4 * exp((-x – 55) / 18)

alphan(x) = 0.01 * (x + 34) / (1 – exp(-0.1 * (x + 34)))

betan(x) = 0.125 * exp((-x – 44) / 80)

alphah(x) = 0.07 * exp((-x – 44) / 20)

betah(x) = 1 / (1 + exp(-0.01 * (x + 14)))

sigma(x) = 1 / (1 + exp((-x + 20) / 4))

minf(x) = alpham(x) / (alpham(x) + betam(x))

Iexte1 = cs0 + cs1 * I1 + cs2 * I1 ^ 2 + cs3 * I1 ^ 3 + cs4 * I1 ^ 4

Iexte2 = cs0 + cs1 * I2 + cs2 * I2 ^ 2 + cs3 * I2 ^ 3 + cs4 * I2 ^ 4

#Initial conditions

init Ve1 = -80

init Ve2 = -80

init Vi1 = -80

init Vi2 = -80