/*
 *
 *  Koivumaki Fibroblast Model
 *  Modified by: Jorge Sanchez Arciniegas
 *
 * Koivumäki, J., Clark, R. B., Belke, D., Kondo, C., Fedak, P., Maleckar, M. M., et al. (2014). 
 * Na+ current expression in human atrial myofibroblasts: Identity and functional roles. 
 * Front. Physiol. 5 JUL, 1–14. doi:10.3389/fphys.2014.00275
 *
 *
 */

#include<stdio.h>
#include<math.h>
#include<vector>

using namespace std;

void KOIFb(vector<double>& Y, double* currents, double dt, double i_Stim, int rmp_slt, int ina_slt) {

	//Select model RMP
	switch (rmp_slt) {
	    case 1: // -35 mV
	    	Kvf_mod = 1.75; K1f_mod = 3; NaKf_mod = 2.75; Nabf_mod = 8.1275; rkv_mod = 15; skv_mod = 18;
	        break;
	    case 2: // -45 mV
	    	Kvf_mod = 0.75; K1f_mod = 0.282; NaKf_mod = 1; Nabf_mod = 1.75; rkv_mod = 20.8; skv_mod = 23;
	        break;
	    case 3: // -65 mV
	    	Kvf_mod = 0.5; K1f_mod = 0.07; NaKf_mod = 0.35; Nabf_mod = 0.44; rkv_mod = 20; skv_mod = 23;
	        break;
	 }

	// Select model variant (different Nav1.5 properties)
	switch (ina_slt) {
	    case 1://'INa_Fb_Poulet'
	        P_Na_Fb = 0.0000225; Na_m_Vhalf_Fb = 29.9; Na_m_hill_Fb = -3.1; Na_h_Vhalf_Fb = 64.7; Na_h_hill_Fb = 6.9;
	        break;
	    case 2://'INa_Fb_Koivumaki'
	        P_Na_Fb = 0.0000125; Na_m_Vhalf_Fb = 16; Na_m_hill_Fb = -8.21; Na_h_Vhalf_Fb = 65; Na_h_hill_Fb = 8.4;
	        break;
	    case 3://'INa_Fb_Chatelier'
	        P_Na_Fb = 0.0000225; Na_m_Vhalf_Fb = 41.7; Na_m_hill_Fb = -7.6; Na_h_Vhalf_Fb = 83.9; Na_h_hill_Fb = 7.8;
	        break;
	    case 4://'no_INa'
	        P_Na_Fb = 0; Na_m_Vhalf_Fb = 41.7; Na_m_hill_Fb = -7.6; Na_h_Vhalf_Fb = 83.9; Na_h_hill_Fb = 7.8;
	        break;
	 }


	// Numerical parameters

	const double T = 306.15;   // kelvin (in membrane)
	const double F = 96487.0;   // coulomb_per_mole (in membrane)
	const double R = 8314.0;   // millijoule_per_mole_kelvin (in membrane)

	// Temperature dependencies
	const double q10exp = (T - 306.15)/10;
	const double q10PNa = 1.45; // Milburn et al. 1995
	const double q10IK1 = 1.4;
	const double q10NKA = 1.63; // Grandi et al. 2011
	const double q10tauNa = 3;

	// Fibroblast parameters

	const double Cm_Fb = 6.3;//pF

	const double Ko_Fb = 5.3581;
	const double Nao_Fb = 130.0110;

	const double Vol_i_Fb = 0.00137;

	const double B_Fb = -200.0;
	const double V_rev_Fb = -150.0;
	const double K_mK_Fb = 1.0;
	const double K_mNa_Fb = 11.0;

	double g_Kv_Fb = 0.37e-3*Cm_Fb*Kvf_mod;
	double g_K1_Fb = 0.25e-3*Cm_Fb*K1f_mod;
	double g_NKA_Fb = 5.39e-3*Cm_Fb*NaKf_mod;
	double g_Nab_Fb = 0.0217e-3*Cm_Fb*Nabf_mod;

	// Analytical equations for the fibroblast

	// Nernst Potentials
	double E_K_Fb = R*T/F * log(Ko_Fb/Y[1]);
	double E_Na_Fb = R*T/F * log(Nao_Fb/Y[2]);

	// INa
	double INaminf_Fb = 1/(1 + exp((Y[0] + Na_m_Vhalf_Fb - (4.3 + 5.0)/20*(T - 294.15))/Na_m_hill_Fb));
	double INahinf_Fb = 1/(1 + exp((Y[0] + Na_h_Vhalf_Fb - (4.7 + 11)/20*(T - 294.15))/Na_h_hill_Fb));
	double INamtau_Fb = 1.0/pow(q10tauNa,q10exp) * (0.000042*exp( -1.0*pow(((Y[0]+25.57)/28.8),2) ) + 0.000024);
	double INah1tau_Fb = 1.0/pow(q10tauNa,q10exp) * (0.03/(1+exp((Y[0]+35.1)/3.2)) + 0.0003);
	double INah2tau_Fb = 1.0/pow(q10tauNa,q10exp) * (0.12/(1+exp((Y[0]+35.1)/3.2)) + 0.003);

	// Time and voltage dependent K+ current, IKv
	double Kv_r_tau_Fb = 0.0203 + 0.1380*exp(-((Y[0] + 20)/25.9)*((Y[0] + 20)/25.9));
	double Kv_s_tau_Fb = 1.574 + 5.268*exp(-((Y[0] + 23)/22.7)*((Y[0] + 23)/22.7));
	double Kv_r_inf_Fb = 1/(1 + exp(-(Y[0] + rkv_mod)/11));//20 - 20
	double Kv_s_inf_Fb = 1/(1 + exp((Y[0] + skv_mod)/7));//23 - 20

	// INa
	double INa_Fb = pow(q10PNa,q10exp) * P_Na_Fb * pow(Y[3],3) * (0.9*Y[4] + 0.1*Y[5]) 
					* Nao_Fb * Y[0] * pow(F,2)/(R*T) * (exp((Y[0]-E_Na_Fb)*F/R/T) - 1) / (exp(Y[0]*F/R/T) - 1);

	// Time and voltage dependent K+ current, IKv
	double IKv_Fb = g_Kv_Fb*Y[6]*Y[7]*(Y[0] - E_K_Fb);

	// Inward-rectifier K+ current, IK1
	double IK1_Fb = pow(q10IK1,q10exp) * g_K1_Fb * pow(Y[1],0.4457) * (Y[0] - E_K_Fb) 
					/ (1 + exp(1.5*(Y[0] - E_K_Fb + 3.6)*F/R/T));

	// Na+-K+ ATPase current, INKA
	double INKA_Fb = pow(q10NKA,q10exp) * g_NKA_Fb*(Ko_Fb/(Ko_Fb + K_mK_Fb)) * ((pow(Y[2], 1.5))/((pow(Y[2], 1.5)) 
					+ (pow(K_mNa_Fb,1.5))))*(Y[0] - V_rev_Fb)/(Y[0] - B_Fb);

	// Na+ background current
	double INab_Fb = g_Nab_Fb*(Y[0] - E_Na_Fb);

	// Intracellular Ion Concentrations
	Y[1] += dt*(-(IK1_Fb + IKv_Fb - 2*INKA_Fb )/(Vol_i_Fb*F)); 
	Y[2] += dt*(-(INa_Fb + INab_Fb + 3*INKA_Fb )/(Vol_i_Fb*F)); 

	// INa
	Y[3] = INaminf_Fb - (INaminf_Fb - Y[3])*exp(-dt/INamtau_Fb);
	Y[4] = INahinf_Fb - (INahinf_Fb - Y[4])*exp(-dt/INah1tau_Fb);
	Y[5] = INahinf_Fb - (INahinf_Fb - Y[5])*exp(-dt/INah2tau_Fb);

	// Time and voltage dependent K+ current, IKv
	Y[6] = Kv_r_inf_Fb - (Kv_r_inf_Fb - Y[6])*exp(-dt/Kv_r_tau_Fb);
	Y[7] = Kv_s_inf_Fb - (Kv_s_inf_Fb - Y[7])*exp(-dt/Kv_s_tau_Fb);

	// Membrane voltage
	Y[0] += dt*(-(INa_Fb + IKv_Fb + IK1_Fb + INab_Fb + INKA_Fb + i_Stim)/Cm_Fb); // currents are in (pA)
	*currents = INa_Fb + IKv_Fb + IK1_Fb + INab_Fb + INKA_Fb + i_Stim;

}