%% Example code for BENG 260 Homework 2, Fall 2019
% Questions to mwagner@ucsd.edu or jmford@ucsd.edu

%% Morris-Lecar (barnacle muscle fiber) Model
% Constants
C_m  = 1.0; % membrane capacitance, in uF/cm^2
g_Ca = 1.1; % maximum conducances, in mS/cm^2
g_K  = 2.0;
g_L  = 0.5;
E_Ca = 100; % Nernst reversal potentials, in mV
E_K  = -70;
E_L  = -50;

% Channel gating kinetics
% Functions of membrane voltage
m_infty = @(V) (1 + tanh((V + 1) ./ 15)) ./ 2;
w_infty = @(V) (1 + tanh(V ./ 30)) ./ 2;
tau_w   = @(V) 5 ./ cosh(V ./ 60);  % in ms

% Membrane currents (in uA/cm^2)
I_Ca = @(V)    g_Ca .* m_infty(V) .* (V - E_Ca);
I_K  = @(V, w) g_K  .* w          .* (V - E_K);
I_L  = @(V)    g_L                .* (V - E_L);

% External current
% step up 10 uA/cm^2 every 100ms
I_ext = @(t) 10 .* floor(t ./ 100);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Problem 2.1: Channel gating kinetics
Vsweep = -20 : 0.25 : 120;
figure();
title('2.1: Channel Gating Kinetics');
% recommended style: use same colors for same state variable (m, n, h) and solid/dashed for alpha/beta
plot(Vsweep, m_infty(Vsweep), 'r-');
hold on;
plot(Vsweep, w_infty(Vsweep), 'g-');
plot(Vsweep, tau_w(Vsweep), 'g--');
hold off;
xlabel('V (mV)');
ylabel('Kinetics Value');
xlim([Vsweep(1) Vsweep(end)]);
legend('m_\infty', 'w_\infty', 'tau_w', 'Location', 'SouthEast');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Problem 2.2: leaky passive neuron
dVdt_leak = @(t, V) (I_ext(t) - I_L(V)) ./ C_m;
[t_leak, V_leak] = ode45(dVdt_leak, [0 400], E_L);

figure();

subplot(2,1,1);
plot(t_leak, V_leak, 'k');
title('2.2: Leaky Passive Neuron');
ylabel('V (mV)');

subplot(2,1,2);
plot(t_leak, I_ext(t_leak), 'k');
xlabel('t (ms)');
ylabel('I_{ext} (\mu{A}/cm^2)');
ylim([-1 max(I_ext(t_leak))+1]);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Problem 2.3: Morris-Lecar neuron
% vector coding of state variables: X = [V, w]
dXdt = @(t, X) ...
    [(I_ext(t) - I_Ca(X(1)) - I_K(X(1), X(2)) - I_L(X(1))) ./ C_m; ...
    (w_infty(X(1)) - X(2)) ./ tau_w(X(1))];

[t, X] = ode23(dXdt, [0 400], [-44 0.05]);

V = X(:,1); % the first column is the V values
w = X(:,2); % the second column is the w values

figure;

subplot(3,1,1);
plot(t, V, 'k');
title('2.3: Morris-Lecar Neuron');
ylabel('V (mV)');

subplot(3,1,2);
% recommended style: use the came colors as you did in part 2.1
plot(t, w, 'g');
ylabel('Gating Value');
legend('w');

subplot(3,1,3);
plot(t, I_ext(t), 'k');
xlabel('t (ms)');
ylabel('I_{ext} (\mu{A}/cm^2)');
ylim([-1 max(I_ext(t))+1]);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Problem 4: Bonus: Morris-Lecar model with expanded channel gating dynamics  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constants
C_m  = 1.0; % membrane capacitance, in uF/cm^2
g_Ca = 1.1; % maximum conducances, in mS/cm^2
g_K  = 2.0;
g_L  = 0.5;
E_Ca = 100; % Nernst reversal potentials, in mV
E_K  = -70;
E_L  = -50;

% Channel gating kinetics
% Functions of membrane voltage
m_infty = @(V) (1 + tanh((V + 1) ./ 15)) ./ 2;
w_infty = @(V) (1 + tanh(V ./ 30)) ./ 2;
tau_w   = @(V) 5 ./ cosh(V ./ 60);  % in ms

% Membrane currents (in uA/cm^2)
I_Ca = @(V)     g_Ca .* m_infty(V) .* (V - E_Ca);
I_K  = @(V, w)  g_K  .* w          .* (V - E_K);
I_L  = @(V)     g_L                .* (V - E_L);
I_ext = @(t)    40+floor(t/20);

t_start = 0;
t_stop = 200;
t_step = 0.1;
t = [t_start:t_step:t_stop];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4.1 Numerical solution (ODE solver)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Differential equation for X(1)=Vm and X(2)=w
dXdt = @(t, X) ...
    [(I_ext(t) - I_K(X(1),X(2)) - I_Ca(X(1)) - I_L(X(1))) ./ C_m; ...
    (w_infty(X(1))-X(2)) ./ tau_w(X(1))];

[t, X] = ode23(dXdt, t, [-70 0]);

Vm_numerical = X(:,1); % the first column is the V values
w_numerical = X(:,2); % the second column is the w values

figure();
plot(t, Vm_numerical, 'k');
figure();
plot(t, w_numerical);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4.2 Stochastic solution (Markov)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha_w = @(Vm) w_infty(Vm) ./ tau_w(Vm);
beta_w  = @(Vm) (1./tau_w(Vm)) - alpha_w(Vm);

N = 1e3;                        % number of gates
gates = zeros(1,N);             % all gates start in state 0
Vm_stoch = zeros(size(t));      % store Vm
w_stoch = zeros(size(t));       % store the fraction of gates open

Vm_current = -70;
w_current = 0;

for t_index=1:1:length(t)
    
    % ODE solver
    dVdt = @(t, Vm) ...
    (I_ext(t) - I_K(Vm,w_current) - I_Ca(Vm) - I_L(Vm)) ./ C_m;

    [t_ode, Vm] = ode23(dVdt, [t(t_index):t_step/10:t(t_index)+t_step], Vm_current);
    
    Vm_current = Vm(end);
    Vm_stoch(t_index) = Vm_current;
    
    for gate_index=1:length(gates)
        prob_sample = rand;             %probabilistic sample
        if gates(gate_index)==0
            % transition from state 0->1
            gates(gate_index) = prob_sample < t_step*alpha_w(Vm_current);
        else
            % remain in state 1
            gates(gate_index) = prob_sample < (1-t_step*beta_w(Vm_current));
        end
    end
    
    w_current = sum(gates)/N;
    w_stoch(t_index) = w_current;
    
end

figure;
plot(t,w_stoch);
figure;
plot(t,Vm_stoch);
figure;
plot(t,I_ext(t));
