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

%% Setup constants
V_t = 26.7; % thermal voltage, in mV
C_m = 1;    % membrane capacitance, in uF/cm^2

% Free ion concentrations of Na+, K+, and Cl-
C_Na_o = 140;  % in mM
C_Na_i = 11;
C_K_o  = 6;
C_K_i  = 145;
C_Cl_o = 111;
C_Cl_i = 5;

% Conductances
g_Na = 13; % in mS/cm^2
g_K  = 5;
g_Cl = 0.15;

% Nernst reversal potentials
E_Na = V_t * log(C_Na_o / C_Na_i); % in mV
E_K  = V_t * log(C_K_o  / C_K_i );
E_Cl = V_t * log(C_Cl_i / C_Cl_o); % reversed since Cl is negatively charged

% Permeabilities
P_Na = g_Na * V_t / E_Na * (1 / C_Na_i - 1 / C_Na_o); % in mS/(cm^2*mM)
P_K  = g_K  * V_t / E_K  * (1 / C_K_i  - 1 / C_K_o );
P_Cl =-g_Cl * V_t / E_Cl * (1 / C_Cl_i - 1 / C_Cl_o);


%% Part 1: Linear and nonlinear current approximations %%
V_m = -150:150; % in mV

% Calculate current densities for all V_m
J_Na_lin = g_Na * (V_m - E_Na);  % in uA/cm^2
J_K_lin  = g_K  * (V_m - E_K );
J_Cl_lin = g_Cl * (V_m - E_Cl);

figure()
%plot(V_m, 0*V_m, 'k:');
hold on
% recommended style: linear uses dashes, GHK uses solid lines, same color for the same ions
plot(V_m, J_Na_lin, 'r--')
plot(V_m, J_K_lin, 'g--')
plot(V_m, J_Cl_lin, 'b--')
xlabel('V_m (mV)') % V_m in mV as per neuroscience tradition
ylabel('J (\mu A/cm^2)')
title('Problem 1');
%legend('hide')
legend('J_{Na} (linear)', 'J_K (linear)', 'J_{Cl} (linear)', 'Location', 'SouthEast')
hold off


%% Part 2: The resting potential %%
% Derive these equations. You will either need to write the equations write in
% Microsoft Word, Latex, or submit a hard copy of hand-written derivations.

% Computing the resting potentials, in mV
V_r_GHK = V_t*log((P_Na*C_Na_o + P_K*C_K_o + P_Cl*C_Cl_i) / (P_Na*C_Na_i + P_K*C_K_i + P_Cl*C_Cl_o));
V_r_lin = (g_Na*E_Na + g_K*E_K + g_Cl*E_Cl)/(g_Na + g_K + g_Cl);

%% Part 3: GHK membrane dynamics %%

% Linear case

% Na+ is active only during first 500 ms
g_Na_t = @(t) g_Na * (t <= 0.5);
% K+ is active during first 250 ms and last 500 ms
g_K_t  = @(t) g_K  * ((t <= 0.25) + (t >= 0.5));
% Cl- is active during first 250 ms and last 500 ms
g_Cl_t = @(t) g_Cl * ((t <= 0.25) + (t >= 0.5));

% Time derivative of V_m, as a function of t and V_m
dVmdt_lin = @(t,V_m) -1/C_m * (g_Na_t(t)*(V_m-E_Na) + g_K_t(t)*(V_m-E_K) + g_Cl_t(t)*(V_m-E_Cl));
[t_lin, V_m_lin] = ode45(dVmdt_lin, [0 1], V_r_lin);

% Plotting
figure;

subplot(2,1,1) % top 1/2 panel, V_m
plot(t_lin, 0*t_lin, 'k:');
hold on
plot(t_lin, V_m_lin, 'k--')
ylabel('V_m (mV)')
legend('', 'linear', 'Location', 'SouthEast');
title('Problem 3');
hold off

subplot(6,1,4) % bottom 4/6 panel, g_K
plot(t_lin, g_K_t(t_lin), 'k');
ylim([-0.1, g_K+0.2]);
ylabel('g_K (mS/cm^2)');

subplot(6,1,5) % bottom 5/6 panel, g_Na
plot(t_lin, g_Na_t(t_lin), 'k');
ylim([-0.1, g_Na+0.2]);
ylabel('g_{Na} (mS/cm^2)');

subplot(6,1,6) % bottom 6/6 panel, g_Cl
plot(t_lin, g_Cl_t(t_lin), 'k');
ylim([-0.1, g_Cl+0.2]);
ylabel('g_{Cl} (mS/cm^2)');

xlabel('Time (s)')
