v=-140:0.01:50;
na = 0:0.01:100;
ca = 0:0.01:2;

thetaM=-35; sigmaM=-5;
thetaA=-20; sigmaA=-10;
thetaS=-20; sigmaS=-1;
thetaMP=-40; sigmaMP=-6;
thetaaT=-65; sigmaaT=-7.8;
thetaRF=-105; sigmaRF=5;
thetaN=-30; sigmaN=-5;
thetaRS=-105; sigmaRS=25;
thetaHP=-48; sigmaHP=6;
thetaE=-60; sigmaE=5;
thetarT = -67; sigmarT = 2; 
thetabT = 0.4; sigmabT = -0.1;
pmax=0.37;
ec=38.7;
nh=3.5;
ks=0.5;


minf = 1./(1+exp((v-thetaM)/sigmaM));
ainf = 1./(1+exp((v-thetaA)/sigmaA));
sinf = 1./(1+exp((v-thetaS)/sigmaS));
mpinf = 1./(1+exp((v-thetaMP)/sigmaMP));
aTinf = 1./(1+exp((v-thetaaT)/sigmaaT));
rinff = 1./(1+exp((v-thetaRF)/sigmaRF));
ninf = 1./(1+exp((v-thetaN)/sigmaN));
rinfs = 1./(1+exp((v-thetaRS)/sigmaRS));
hpinf = 1./(1+exp((v-thetaHP)/sigmaHP));
einf = 1./(1+exp((v-thetaE)/sigmaE));
alphah=0.128*exp(-(v+50)/18);
betah=4./(1+exp(-(v+27)/5));
hinf=alphah./(alphah + betah);
rTinf = 1./(1+exp((v-thetarT)/sigmarT));
bTinf = 1./(1+exp((rTinf-thetabT)/sigmabT))-1./(1+exp(-thetabT/sigmabT));
winf = pmax./(1+(ec./na).^nh);
kinf = (ca.^2)./(ca.^2+(ks^2));

% Panel A
figure();
plot(v, minf, 'k-', 'LineWidth', 2.5);
hold on
plot(v, ainf, 'b-', 'LineWidth', 2.5);
plot(v, sinf, 'c-', 'LineWidth', 2.5);
plot(v, aTinf, 'm-', 'LineWidth', 2.5);
plot(v, rinff, 'r-', 'LineWidth', 2);
plot(v, mpinf, 'g-', 'LineWidth', 2.5);
legend('m_{\infty}', ...
    'a_{\infty}', ...
    's_{\infty}', ...
     'a_{T_{\infty}}', ...
    'rf_{\infty}', ...
    'mp_{\infty}');
xlabel('Voltage(mV)');
ylabel('Gating');



% Panel B
figure();
plot(v, ninf, 'm-', 'LineWidth', 2.5);
hold on
plot(v, rinfs, 'r-', 'LineWidth', 2.5);
plot(v, einf, 'b--', 'LineWidth', 2.5);
plot(v, hinf, 'k--', 'LineWidth', 2.5);
plot(v, bTinf, 'c--', 'LineWidth', 2.5);
plot(v, hpinf, 'g--', 'LineWidth', 2.5);
legend('n_{\infty}', ...
    'rs_{\infty}', ...
    'e_{\infty}', ...
    'h_{\infty}', ...
    'b_{T_{\infty}}', ...
    'hp_{\infty}');
xlabel('Voltage(mV)');
ylabel('Gating');


% Panel C
figure();
plot(ca, kinf, 'k-', 'LineWidth', 2);
xlabel('[Ca^{++}] (\muM)');
ylabel('Steady-state activation function');


% Panel D
figure();
plot(na, winf, 'k-', 'LineWidth', 2);
xlabel('[Na^+] (mM)');
ylabel('Steady-state activation function');
