From be9c28fca12979294c536edc3e16e3fa26ba3dd8 Mon Sep 17 00:00:00 2001 From: bergolho Date: Fri, 13 Oct 2023 10:51:21 +0100 Subject: [PATCH] Added 'courtemanche-ramirez-1998-RushLarsen'/Changed PMJ delay calculus to use the minLAT of the coupled ventricular cells --- ToDo | 5 +- .../activation_time_and_apd_example.ini | 64 +++ scripts/reader_acm.py | 69 +++ ...urtemanche_ramirez_nattel_1998_Victor_RL.c | 363 +++++++++++++++ ...rtemanche_ramirez_nattel_1998_Victor_RL.cu | 413 ++++++++++++++++++ ...urtemanche_ramirez_nattel_1998_Victor_RL.h | 30 ++ ...ramirez_nattel_1998_Victor_RL_common.inc.c | 313 +++++++++++++ src/monodomain/monodomain_solver.c | 16 +- src/save_mesh_library/save_mesh_purkinje.c | 14 +- 9 files changed, 1279 insertions(+), 8 deletions(-) create mode 100644 example_configs/activation_time_and_apd_example.ini create mode 100644 scripts/reader_acm.py create mode 100644 src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL.c create mode 100644 src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL.cu create mode 100644 src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL.h create mode 100644 src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL_common.inc.c diff --git a/ToDo b/ToDo index 4f300149..35e30af2 100644 --- a/ToDo +++ b/ToDo @@ -32,4 +32,7 @@ When the minimum number of PMJs is not reached the solver will be in an infinite ... . . . . . . . . . . -. . . . . \ No newline at end of file +. . . . . + + +Fix the "save_multiple_cell_state_variables". It cannot find the proper cell center. \ No newline at end of file diff --git a/example_configs/activation_time_and_apd_example.ini b/example_configs/activation_time_and_apd_example.ini new file mode 100644 index 00000000..6036e745 --- /dev/null +++ b/example_configs/activation_time_and_apd_example.ini @@ -0,0 +1,64 @@ +[main] +num_threads=4 +dt_pde=0.01 +simulation_time=5000.0 +abort_on_no_activity=false +use_adaptivity=false + +[update_monodomain] +main_function=update_monodomain_default + +[save_result] +print_rate=10 +mesh_print_rate=100 +mesh_format=ensight +output_dir=./outputs/cable_save_activation_time_and_apd_example +init_function=init_save_with_activation_times +end_function=end_save_with_activation_times +main_function=save_with_activation_times +time_threshold=0.0 +activation_threshold=-50.0 +apd_threshold=-70.0 +save_visible_mask=false +remove_older_simulation=true + +[assembly_matrix] +init_function=set_initial_conditions_fvm +sigma_x=0.0000176 +sigma_y=0.0000176 +sigma_z=0.0000176 +library_file=shared_libs/libdefault_matrix_assembly.so +main_function=homogeneous_sigma_assembly_matrix + +[linear_system_solver] +tolerance=1e-16 +use_preconditioner=no +max_iterations=500 +library_file=shared_libs/libdefault_linear_system_solver.so +use_gpu=no +main_function=conjugate_gradient +init_function=init_conjugate_gradient +end_function=end_conjugate_gradient + +[domain] +name=Cable Mesh with no fibrosis +start_dx=100.0 +start_dy=100.0 +start_dz=100.0 +cable_length=10000.0 +main_function=initialize_grid_with_cable_mesh + +[ode_solver] +adaptive=false +dt=0.01 +use_gpu=false +gpu_id=0 +library_file= shared_libs/libToRORd_dynCl_mixed_endo_mid_epi.so + +[stim_plain] +start = 0.0 +duration = 2.0 +current = -53.0 +x_limit = 500.0 +period=1000.0 +main_function=stim_if_x_less_than diff --git a/scripts/reader_acm.py b/scripts/reader_acm.py new file mode 100644 index 00000000..6d26f898 --- /dev/null +++ b/scripts/reader_acm.py @@ -0,0 +1,69 @@ +import sys + +def read_acm_file (filename, target_center_x, target_center_y, target_center_z): + file = open(filename) + counter_lines = 0 + for line in file: + if (counter_lines > 1): + line = line.strip() + new_line = line.replace(",", " ") + + # Cell geometry section + first_open_bracket_id = new_line.find('[') # find() -> lowest index + first_close_bracket_id = new_line.find(']') + second_open_bracket_id = new_line.rfind('[') # rfind() -> hightest index + second_close_bracket_id = new_line.rfind(']') + cell_geometry_data = new_line[0:first_open_bracket_id] + tokens = cell_geometry_data.split() + center_x, center_y, center_z, dx, dy, dz = float(tokens[0]), float(tokens[1]), float(tokens[2]), float(tokens[3]), float(tokens[4]), float(tokens[5]) + #print("%g %g %g %g %g %g" % (center_x, center_y, center_y, dx, dy, dz)) + + # Activation time section + activation_time_values = [] + activation_time_data = new_line[first_open_bracket_id+1:first_close_bracket_id].split() + for value in activation_time_data: + activation_time_values.append(float(value)) + #print(activation_time_values) + + # APD section + apd_values = [] + apd_data = new_line[second_open_bracket_id+1:second_close_bracket_id].split() + for value in apd_data: + apd_values.append(float(value)) + #print(apd_values) + + # SUCESS: Found the target cell! + if (center_x == target_center_x and center_y == target_center_y and center_z == target_center_z): + return activation_time_values, apd_values + + counter_lines = counter_lines + 1 + file.close() + return activation_time_values, apd_values + +def main(): + + if len(sys.argv) != 5: + print("-------------------------------------------------------------------------") + print("Usage:> python %s " % sys.argv[0]) + print("-------------------------------------------------------------------------") + print(" = Input \".acm\" file with the activation time and APD data") + print(" = Cell center x position") + print(" = Cell center y position") + print(" = Cell center z position") + print("-------------------------------------------------------------------------") + return 1 + + input_file = sys.argv[1] + target_center_x = float(sys.argv[2]) + target_center_y = float(sys.argv[3]) + target_center_z = float(sys.argv[4]) + + activation_time_values, apd_values = read_acm_file(input_file, target_center_x, target_center_y, target_center_z) + print("Activation times:") + print(activation_time_values) + print() + print("APDs:") + print(apd_values) + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL.c b/src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL.c new file mode 100644 index 00000000..1051c17f --- /dev/null +++ b/src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL.c @@ -0,0 +1,363 @@ +#include "courtemanche_ramirez_nattel_1998_Victor_RL.h" +#include + +GET_CELL_MODEL_DATA(init_cell_model_data) { + + if(get_initial_v) + cell_model->initial_v = INITIAL_V; + if(get_neq) + cell_model->number_of_ode_equations = NEQ; +} + +SET_ODE_INITIAL_CONDITIONS_CPU(set_model_initial_conditions_cpu) { + + log_info("Using courtemanche_ramirez_nattel_1998_Victor_RL CPU model\n"); + + uint32_t num_cells = solver->original_num_cells; + solver->sv = (real*)malloc(NEQ*num_cells*sizeof(real)); + + bool adpt = solver->adaptive; + + if(adpt) { + solver->ode_dt = (real*)malloc(num_cells*sizeof(real)); + + OMP(parallel for) + for(int i = 0; i < num_cells; i++) { + solver->ode_dt[i] = solver->min_dt; + } + + solver->ode_previous_dt = (real*)calloc(num_cells, sizeof(real)); + solver->ode_time_new = (real*)calloc(num_cells, sizeof(real)); + log_info("Using Adaptive Euler model to solve the ODEs\n"); + } else { + log_info("Using Rush-Larsen/Euler model to solve the ODEs\n"); + } + + OMP(parallel for) + for(uint32_t i = 0; i < num_cells; i++) { + + real *sv = &solver->sv[i * NEQ]; + if (solver->ode_extra_data == NULL) { + sv[0] = -8.118000e+01f; //V millivolt + sv[1] = 2.908000e-03f; //m dimensionless + sv[2] = 9.649000e-01f; //h dimensionless + sv[3] = 9.775000e-01f; //j dimensionless + sv[4] = 3.043000e-02f; //oa dimensionless + sv[5] = 9.992000e-01f; //oi dimensionless + sv[6] = 4.966000e-03f; //ua dimensionless + sv[7] = 9.986000e-01f; //ui dimensionless + sv[8] = 3.296000e-05f; //xr dimensionless + sv[9] = 1.869000e-02f; //xs dimensionless + sv[10] = 1.367000e-04f; //d dimensionless + sv[11] = 9.996000e-01f; //f dimensionless + sv[12] = 7.755000e-01f; //f_Ca dimensionless + sv[13] = 0.000000e+00f; //u dimensionless + sv[14] = 1.000000e+00f; //v dimensionless + sv[15] = 9.992000e-01f; //w dimensionless + sv[16] = 1.117000e+01f; //Na_i millimolar + sv[17] = 1.390000e+02f; //K_i millimolar + sv[18] = 1.013000e-04f; //Ca_i millimolar + sv[19] = 1.488000e+00f; //Ca_up millimolar + sv[20] = 1.488000e+00f; //Ca_rel millimolar + } + else { + real *extra_parameters = (real *)solver->ode_extra_data; + sv[ 0] = extra_parameters[13]; + sv[ 1] = extra_parameters[14]; + sv[ 2] = extra_parameters[15]; + sv[ 3] = extra_parameters[16]; + sv[ 4] = extra_parameters[17]; + sv[ 5] = extra_parameters[18]; + sv[ 6] = extra_parameters[19]; + sv[ 7] = extra_parameters[20]; + sv[ 8] = extra_parameters[21]; + sv[ 9] = extra_parameters[22]; + sv[10] = extra_parameters[23]; + sv[11] = extra_parameters[24]; + sv[12] = extra_parameters[25]; + sv[13] = extra_parameters[26]; + sv[14] = extra_parameters[27]; + sv[15] = extra_parameters[28]; + sv[16] = extra_parameters[29]; + sv[17] = extra_parameters[30]; + sv[18] = extra_parameters[31]; + sv[19] = extra_parameters[32]; + sv[20] = extra_parameters[33]; + } + } +} + +SOLVE_MODEL_ODES(solve_model_odes_cpu) { + + uint32_t sv_id; + + size_t num_cells_to_solve = ode_solver->num_cells_to_solve; + uint32_t * cells_to_solve = ode_solver->cells_to_solve; + real *sv = ode_solver->sv; + real dt = ode_solver->min_dt; + uint32_t num_steps = ode_solver->num_steps; + + int offset = 34; + bool adpt = ode_solver->adaptive; + real *extra_parameters = NULL; + + if(ode_solver->ode_extra_data) { + extra_parameters = ((real*)ode_solver->ode_extra_data); + } + else { + log_error_and_exit("You need to specify a mask function when using this mixed model!\n"); + } + + #pragma omp parallel for private(sv_id) + for (int i = 0; i < num_cells_to_solve; i++) { + int mapping = extra_parameters[i+offset]; + + if(cells_to_solve) + sv_id = cells_to_solve[i]; + else + sv_id = i; + + if(adpt) { + solve_forward_euler_cpu_adpt(sv + (sv_id * NEQ), stim_currents[i], current_t + dt, sv_id, ode_solver, ode_solver->ode_extra_data, mapping); + } + else { + for (int j = 0; j < num_steps; ++j) { + solve_model_ode_cpu(dt, sv + (sv_id * NEQ), stim_currents[i], ode_solver->ode_extra_data, mapping); + } + } + } +} + +void solve_model_ode_cpu(real dt, real *sv, real stim_current, real *extra_parameters, int mapping) { + + real rY[NEQ], rDY[NEQ]; + + for(int i = 0; i < NEQ; i++) + rY[i] = sv[i]; + + RHS_cpu(rY, rDY, stim_current, dt, extra_parameters, mapping); + + // This model uses the Rush-Larsen method to solve the ODEs + sv[0] = dt*rDY[0] + rY[0]; // Euler + + sv[1] = rDY[1]; // Rush-Larsen + sv[2] = rDY[2]; // Rush-Larsen + sv[3] = rDY[3]; // Rush-Larsen + sv[4] = rDY[4]; // Rush-Larsen + sv[5] = rDY[5]; // Rush-Larsen + sv[6] = rDY[6]; // Rush-Larsen + sv[7] = rDY[7]; // Rush-Larsen + sv[8] = rDY[8]; // Rush-Larsen + sv[9] = rDY[9]; // Rush-Larsen + sv[10] = rDY[10]; // Rush-Larsen + sv[11] = rDY[11]; // Rush-Larsen + sv[12] = rDY[12]; // Rush-Larsen + sv[13] = rDY[13]; // Rush-Larsen + sv[14] = rDY[14]; // Rush-Larsen + sv[15] = rDY[15]; // Rush-Larsen + + sv[16] = dt*rDY[16] + rY[16]; // Euler + sv[17] = dt*rDY[17] + rY[17]; // Euler + sv[18] = dt*rDY[18] + rY[18]; // Euler + sv[19] = dt*rDY[19] + rY[19]; // Euler + sv[20] = dt*rDY[20] + rY[20]; // Euler + +} + +void solve_forward_euler_cpu_adpt(real *sv, real stim_curr, real final_time, int sv_id, struct ode_solver *solver, real *extra_parameters, int mapping) { + + const real _beta_safety_ = 0.8; + int numEDO = NEQ; + + real rDY[numEDO]; + + real _tolerances_[numEDO]; + real _aux_tol = 0.0; + // initializes the variables + solver->ode_previous_dt[sv_id] = solver->ode_dt[sv_id]; + + real edos_old_aux_[numEDO]; + real edos_new_euler_[numEDO]; + real *_k1__ = (real *)malloc(sizeof(real) * numEDO); + real *_k2__ = (real *)malloc(sizeof(real) * numEDO); + real *_k_aux__; + + real *dt = &solver->ode_dt[sv_id]; + real *time_new = &solver->ode_time_new[sv_id]; + real *previous_dt = &solver->ode_previous_dt[sv_id]; + + if(*time_new + *dt > final_time) { + *dt = final_time - *time_new; + } + + RHS_cpu(sv, rDY, stim_curr, *dt, extra_parameters, mapping); + *time_new += *dt; + + for(int i = 0; i < numEDO; i++) { + _k1__[i] = rDY[i]; + } + + const real rel_tol = solver->rel_tol; + const real abs_tol = solver->abs_tol; + + const real __tiny_ = pow(abs_tol, 2.0); + + real min_dt = solver->min_dt; + real max_dt = solver->max_dt; + + while(1) { + + for(int i = 0; i < numEDO; i++) { + // stores the old variables in a vector + edos_old_aux_[i] = sv[i]; + // computes euler method + edos_new_euler_[i] = _k1__[i] * *dt + edos_old_aux_[i]; + // steps ahead to compute the rk2 method + sv[i] = edos_new_euler_[i]; + } + + *time_new += *dt; + RHS_cpu(sv, rDY, stim_curr, *dt, extra_parameters, mapping); + *time_new -= *dt; // step back + + double greatestError = 0.0, auxError = 0.0; + for(int i = 0; i < numEDO; i++) { + // stores the new evaluation + _k2__[i] = rDY[i]; + _aux_tol = fabs(edos_new_euler_[i]) * rel_tol; + _tolerances_[i] = (abs_tol > _aux_tol) ? abs_tol : _aux_tol; + // finds the greatest error between the steps + auxError = fabs(((*dt / 2.0) * (_k1__[i] - _k2__[i])) / _tolerances_[i]); + + greatestError = (auxError > greatestError) ? auxError : greatestError; + } + /// adapt the time step + greatestError += __tiny_; + *previous_dt = *dt; + /// adapt the time step + *dt = _beta_safety_ * (*dt) * sqrt(1.0f / greatestError); + + if(*dt < min_dt) { + *dt = min_dt; + } else if(*dt > max_dt) { + *dt = max_dt; + } + + if(*time_new + *dt > final_time) { + *dt = final_time - *time_new; + } + + // it doesn't accept the solution + if(greatestError >= 1.0f && *dt > min_dt) { + // restore the old values to do it again + for(int i = 0; i < numEDO; i++) { + sv[i] = edos_old_aux_[i]; + } + // throw the results away and compute again + } else { + // it accepts the solutions + if(greatestError >= 1.0) { + printf("Accepting solution with error > %lf \n", greatestError); + } + + _k_aux__ = _k2__; + _k2__ = _k1__; + _k1__ = _k_aux__; + + // it steps the method ahead, with euler solution + for(int i = 0; i < numEDO; i++) { + sv[i] = edos_new_euler_[i]; + } + + if(*time_new + *previous_dt >= final_time) { + if(final_time == *time_new) { + break; + } else if(*time_new < final_time) { + *dt = *previous_dt = final_time - *time_new; + *time_new += *previous_dt; + break; + } + } else { + *time_new += *previous_dt; + } + } + } + + free(_k1__); + free(_k2__); +} + +void RHS_cpu(const real *sv, real *rDY, real stim_current, real dt, real *extra_parameters, int mapping) { + + //State variables + const real V_old_ = sv[0]; + const real m_old_ = sv[1]; + const real h_old_ = sv[2]; + const real j_old_ = sv[3]; + const real oa_old_ = sv[4]; + const real oi_old_ = sv[5]; + const real ua_old_ = sv[6]; + const real ui_old_ = sv[7]; + const real xr_old_ = sv[8]; + const real xs_old_ = sv[9]; + const real d_old_ = sv[10]; + const real f_old_ = sv[11]; + const real f_Ca_old_ = sv[12]; + const real u_old_ = sv[13]; + const real v_old_ = sv[14]; + const real w_old_ = sv[15]; + const real Na_i_old_ = sv[16]; + const real K_i_old_ = sv[17]; + const real Ca_i_old_ = sv[18]; + const real Ca_up_old_ = sv[19]; + const real Ca_rel_old_ = sv[20]; + + #include "courtemanche_ramirez_nattel_1998_Victor_RL_common.inc.c" +} + + +/* +[extra_data] +; Membrane model parameters +;GKur = 1.0 +;GKr = 1.0 +;Gto = 0.4 +;GK1 = 2.0 +;GKs = 1.0 +;GCaL = 0.35 +;GNaK = 1.0 +;GNCX = 1.0 +;GNa = 1.0 +;GK2P = 0.0 +;Gup = 1.0 +;Grel = 1.0 +;Gleak= 1.0 + +; Initial conditions +;Vm = -81.2 +;INa_va = 2.91e-03 +;INa_vi_1 = 9.65e-01 +;INa_vi_2 = 9.78e-01 +;Ito_va = 3.04e-02 +;Ito_vi = 9.99e-01 +;ICaL_va = 1.37e-04 +;ICaL_vi = 9.99e-01 +;ICaL_ci = 7.75e-01 +;IKur_va = 4.96e-03 +;IKur_viS = 9.99e-01 +;IKur_viF = 9.99e-01 +;IKr_va = 3.29e-05 +;IKs_va = 1.87e-02 +;IK2P_va = 0.202895277583553 +;CajSR = 1.49 +;CanSR = 1.49 +;RyRo = 1.06460418123529e-31 +;RyRr = 1.0 +;RyRi = 9.99e-01 +;Nai = 1.12e+01 +;Nao = 140 +;Ki = 1.39e02 +;Ko = 5.4 +;Cai = 1.02e-04 +;Cao = 1.8 +*/ \ No newline at end of file diff --git a/src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL.cu b/src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL.cu new file mode 100644 index 00000000..954d42f6 --- /dev/null +++ b/src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL.cu @@ -0,0 +1,413 @@ +#include "courtemanche_ramirez_nattel_1998_Victor_RL.h" +#include +#include + +__global__ void kernel_set_model_initial_conditions(real *sv, int num_volumes, size_t pitch, bool use_adpt_dt, real min_dt, real *extra_parameters) { + int threadID = blockDim.x * blockIdx.x + threadIdx.x; + + if (threadID < num_volumes) { + + if (extra_parameters == NULL) { + *((real * )((char *) sv + pitch * 0) + threadID) = -8.118000e+01f; //V millivolt + *((real * )((char *) sv + pitch * 1) + threadID) = 2.908000e-03f; //m dimensionless + *((real * )((char *) sv + pitch * 2) + threadID) = 9.649000e-01f; //h dimensionless + *((real * )((char *) sv + pitch * 3) + threadID) = 9.775000e-01f; //j dimensionless + *((real * )((char *) sv + pitch * 4) + threadID) = 3.043000e-02f; //oa dimensionless + *((real * )((char *) sv + pitch * 5) + threadID) = 9.992000e-01f; //oi dimensionless + *((real * )((char *) sv + pitch * 6) + threadID) = 4.966000e-03f; //ua dimensionless + *((real * )((char *) sv + pitch * 7) + threadID) = 9.986000e-01f; //ui dimensionless + *((real * )((char *) sv + pitch * 8) + threadID) = 3.296000e-05f; //xr dimensionless + *((real * )((char *) sv + pitch * 9) + threadID) = 1.869000e-02f; //xs dimensionless + *((real * )((char *) sv + pitch * 10) + threadID) = 1.367000e-04f; //d dimensionless + *((real * )((char *) sv + pitch * 11) + threadID) = 9.996000e-01f; //f dimensionless + *((real * )((char *) sv + pitch * 12) + threadID) = 7.755000e-01f; //f_Ca dimensionless + *((real * )((char *) sv + pitch * 13) + threadID) = 0.0f; //u dimensionless + *((real * )((char *) sv + pitch * 14) + threadID) = 1.000000e+00f; //v dimensionless + *((real * )((char *) sv + pitch * 15) + threadID) = 9.992000e-01f; //w dimensionless + *((real * )((char *) sv + pitch * 16) + threadID) = 1.117000e+01f; //Na_i millimolar + *((real * )((char *) sv + pitch * 17) + threadID) = 1.390000e+02f; //K_i millimolar + *((real * )((char *) sv + pitch * 18) + threadID) = 1.013000e-04f; //Ca_i millimolar + *((real * )((char *) sv + pitch * 19) + threadID) = 1.488000e+00f; //Ca_up millimolar + *((real * )((char *) sv + pitch * 20) + threadID) = 1.488000e+00f; //Ca_rel millimolar + } + else { + *((real * )((char *) sv + pitch * 0) + threadID) = extra_parameters[13]; + *((real * )((char *) sv + pitch * 1) + threadID) = extra_parameters[14]; + *((real * )((char *) sv + pitch * 2) + threadID) = extra_parameters[15]; + *((real * )((char *) sv + pitch * 3) + threadID) = extra_parameters[16]; + *((real * )((char *) sv + pitch * 4) + threadID) = extra_parameters[17]; + *((real * )((char *) sv + pitch * 5) + threadID) = extra_parameters[18]; + *((real * )((char *) sv + pitch * 6) + threadID) = extra_parameters[19]; + *((real * )((char *) sv + pitch * 7) + threadID) = extra_parameters[20]; + *((real * )((char *) sv + pitch * 8) + threadID) = extra_parameters[21]; + *((real * )((char *) sv + pitch * 9) + threadID) = extra_parameters[22]; + *((real * )((char *) sv + pitch * 10) + threadID) = extra_parameters[23]; + *((real * )((char *) sv + pitch * 11) + threadID) = extra_parameters[24]; + *((real * )((char *) sv + pitch * 12) + threadID) = extra_parameters[25]; + *((real * )((char *) sv + pitch * 13) + threadID) = extra_parameters[26]; + *((real * )((char *) sv + pitch * 14) + threadID) = extra_parameters[27]; + *((real * )((char *) sv + pitch * 15) + threadID) = extra_parameters[28]; + *((real * )((char *) sv + pitch * 16) + threadID) = extra_parameters[29]; + *((real * )((char *) sv + pitch * 17) + threadID) = extra_parameters[30]; + *((real * )((char *) sv + pitch * 18) + threadID) = extra_parameters[31]; + *((real * )((char *) sv + pitch * 19) + threadID) = extra_parameters[32]; + *((real * )((char *) sv + pitch * 20) + threadID) = extra_parameters[33]; + } + if(use_adpt_dt) { + *((real *)((char *)sv + pitch * NEQ) + threadID) = min_dt; // dt + *((real *)((char *)sv + pitch * (NEQ + 1)) + threadID) = 0.0; // time_new + *((real *)((char *)sv + pitch * (NEQ + 2)) + threadID) = 0.0; // previous dt + } + } +} + +inline __device__ void RHS_gpu(real *sv, real *rDY, real stim_current, int thread_id, real dt, size_t pitch, bool use_adpt_dt, real *extra_parameters, int mapping) { + + //State variables + real V_old_; //millivolt + real m_old_; //dimensionless + real h_old_; //dimensionless + real j_old_; //dimensionless + real oa_old_; //dimensionless + real oi_old_; //dimensionless + real ua_old_; //dimensionless + real ui_old_; //dimensionless + real xr_old_; //dimensionless + real xs_old_; //dimensionless + real d_old_; //dimensionless + real f_old_; //dimensionless + real f_Ca_old_; //dimensionless + real u_old_; //dimensionless + real v_old_; //dimensionless + real w_old_; //dimensionless + real Na_i_old_; //millimolar + real K_i_old_; //millimolar + real Ca_i_old_; //millimolar + real Ca_up_old_; //millimolar + real Ca_rel_old_; //millimolar + + if(use_adpt_dt) { + V_old_ = sv[0]; + m_old_ = sv[1]; + h_old_ = sv[2]; + j_old_ = sv[3]; + oa_old_ = sv[4]; + oi_old_ = sv[5]; + ua_old_ = sv[6]; + ui_old_ = sv[7]; + xr_old_ = sv[8]; + xs_old_ = sv[9]; + d_old_ = sv[10]; + f_old_ = sv[11]; + f_Ca_old_ = sv[12]; + u_old_ = sv[13]; + v_old_ = sv[14]; + w_old_ = sv[15]; + Na_i_old_ = sv[16]; + K_i_old_ = sv[17]; + Ca_i_old_ = sv[18]; + Ca_up_old_ = sv[19]; + Ca_rel_old_ = sv[20]; + } else { + V_old_ = *((real*)((char*)sv + pitch * 0) + thread_id); + m_old_ = *((real*)((char*)sv + pitch * 1) + thread_id); + h_old_ = *((real*)((char*)sv + pitch * 2) + thread_id); + j_old_ = *((real*)((char*)sv + pitch * 3) + thread_id); + oa_old_ = *((real*)((char*)sv + pitch * 4) + thread_id); + oi_old_ = *((real*)((char*)sv + pitch * 5) + thread_id); + ua_old_ = *((real*)((char*)sv + pitch * 6) + thread_id); + ui_old_ = *((real*)((char*)sv + pitch * 7) + thread_id); + xr_old_ = *((real*)((char*)sv + pitch * 8) + thread_id); + xs_old_ = *((real*)((char*)sv + pitch * 9) + thread_id); + d_old_ = *((real*)((char*)sv + pitch * 10) + thread_id); + f_old_ = *((real*)((char*)sv + pitch * 11) + thread_id); + f_Ca_old_ = *((real*)((char*)sv + pitch * 12) + thread_id); + u_old_ = *((real*)((char*)sv + pitch * 13) + thread_id); + v_old_ = *((real*)((char*)sv + pitch * 14) + thread_id); + w_old_ = *((real*)((char*)sv + pitch * 15) + thread_id); + Na_i_old_ = *((real*)((char*)sv + pitch * 16) + thread_id); + K_i_old_ = *((real*)((char*)sv + pitch * 17) + thread_id); + Ca_i_old_ = *((real*)((char*)sv + pitch * 18) + thread_id); + Ca_up_old_ = *((real*)((char*)sv + pitch * 19) + thread_id); + Ca_rel_old_ = *((real*)((char*)sv + pitch * 20) + thread_id); + } + + #include "courtemanche_ramirez_nattel_1998_Victor_RL_common.inc.c" + +} + +extern "C" SET_ODE_INITIAL_CONDITIONS_GPU(set_model_initial_conditions_gpu) { + + size_t pitch_h; + + uint8_t use_adpt_dt = (uint8_t)solver->adaptive; + + log_info("Using GPU model implemented in %s\n", __FILE__); + + uint32_t num_volumes = solver->original_num_cells; + + if(use_adpt_dt) { + log_info("Using Adaptive Euler model to solve the ODEs\n"); + } else { + log_info("Using Rush-Larsen/Euler model to solve the ODEs\n"); + } + + // execution configuration + const int GRID = (num_volumes + BLOCK_SIZE - 1) / BLOCK_SIZE; + + size_t size = num_volumes * sizeof(real); + + if(use_adpt_dt) + check_cuda_error(cudaMallocPitch((void **)&(solver->sv), &pitch_h, size, (size_t)NEQ + 3)); + else + check_cuda_error(cudaMallocPitch((void **)&(solver->sv), &pitch_h, size, (size_t)NEQ)); + + real *extra_parameters = NULL; + real *extra_parameters_device = NULL; + + if(solver->ode_extra_data) { + extra_parameters = (real *)solver->ode_extra_data; + check_cuda_error(cudaMalloc((void **)&extra_parameters_device, solver->extra_data_size)); + check_cuda_error(cudaMemcpy(extra_parameters_device, extra_parameters, solver->extra_data_size, cudaMemcpyHostToDevice)); + } + + kernel_set_model_initial_conditions<<>>(solver->sv, num_volumes, pitch_h, use_adpt_dt, solver->min_dt, extra_parameters_device); + + check_cuda_error(cudaPeekAtLastError()); + cudaDeviceSynchronize(); + + return pitch_h; +} + +extern "C" SOLVE_MODEL_ODES(solve_model_odes_gpu) { + + size_t num_cells_to_solve = ode_solver->num_cells_to_solve; + uint32_t * cells_to_solve = ode_solver->cells_to_solve; + real *sv = ode_solver->sv; + real dt = ode_solver->min_dt; + uint32_t num_steps = ode_solver->num_steps; + + // execution configuration + const int GRID = ((int)num_cells_to_solve + BLOCK_SIZE - 1) / BLOCK_SIZE; + + size_t stim_currents_size = sizeof(real) * num_cells_to_solve; + size_t cells_to_solve_size = sizeof(uint32_t) * num_cells_to_solve; + + real *stims_currents_device; + check_cuda_error(cudaMalloc((void **)&stims_currents_device, stim_currents_size)); + check_cuda_error(cudaMemcpy(stims_currents_device, stim_currents, stim_currents_size, cudaMemcpyHostToDevice)); + + // the array cells to solve is passed when we are using and adaptive mesh + uint32_t *cells_to_solve_device = NULL; + if(cells_to_solve != NULL) { + check_cuda_error(cudaMalloc((void **)&cells_to_solve_device, cells_to_solve_size)); + check_cuda_error(cudaMemcpy(cells_to_solve_device, cells_to_solve, cells_to_solve_size, cudaMemcpyHostToDevice)); + } + + // Get the extra_data array + real *extra_data = NULL; + real *extra_data_device = NULL; + if(ode_solver->ode_extra_data) { + extra_data = (real*)ode_solver->ode_extra_data; + } + else { + log_error_and_exit("You need to specify a mask function when using a mixed model!\n"); + } + check_cuda_error(cudaMalloc((void **)&extra_data_device, ode_solver->extra_data_size)); + check_cuda_error(cudaMemcpy(extra_data_device, extra_data, ode_solver->extra_data_size, cudaMemcpyHostToDevice)); + + solve_gpu<<>>(current_t, dt, sv, stims_currents_device, cells_to_solve_device, num_cells_to_solve, + num_steps, ode_solver->pitch, ode_solver->adaptive, ode_solver->abs_tol, + ode_solver->rel_tol, ode_solver->max_dt, extra_data_device); + + check_cuda_error(cudaPeekAtLastError()); + + check_cuda_error(cudaFree(stims_currents_device)); + if(cells_to_solve_device) + check_cuda_error(cudaFree(cells_to_solve_device)); +} + + +inline __device__ void solve_forward_euler_gpu_adpt(real *sv, real stim_curr, real final_time, int thread_id, size_t pitch, real abstol, real reltol, real min_dt, real max_dt, real *extra_parameters, int mapping) { + + #define DT *((real *)((char *)sv + pitch * (NEQ)) + thread_id) + #define TIME_NEW *((real *)((char *)sv + pitch * (NEQ+1)) + thread_id) + #define PREVIOUS_DT *((real *)((char *)sv + pitch * (NEQ+2)) + thread_id) + + real rDY[NEQ]; + + real _tolerances_[NEQ]; + real _aux_tol = 0.0; + real dt = DT; + real time_new = TIME_NEW; + real previous_dt = PREVIOUS_DT; + + real edos_old_aux_[NEQ]; + real edos_new_euler_[NEQ]; + real _k1__[NEQ]; + real _k2__[NEQ]; + real _k_aux__[NEQ]; + real sv_local[NEQ]; + + const real _beta_safety_ = 0.8; + + const real __tiny_ = pow(abstol, 2.0); + + if(time_new + dt > final_time) { + dt = final_time - time_new; + } + + for(int i = 0; i < NEQ; i++) { + sv_local[i] = *((real *)((char *)sv + pitch * i) + thread_id); + } + + RHS_gpu(sv_local, rDY, stim_curr, thread_id, dt, pitch, true, extra_parameters, mapping); + time_new += dt; + + for(int i = 0; i < NEQ; i++) { + _k1__[i] = rDY[i]; + } + + while(1) { + + for(int i = 0; i < NEQ; i++) { + // stores the old variables in a vector + edos_old_aux_[i] = sv_local[i]; + // //computes euler method + edos_new_euler_[i] = _k1__[i] * dt + edos_old_aux_[i]; + // steps ahead to compute the rk2 method + sv_local[i] = edos_new_euler_[i]; + } + + time_new += dt; + + RHS_gpu(sv_local, rDY, stim_curr, thread_id, dt, pitch, true, extra_parameters, mapping); + time_new -= dt; // step back + + real greatestError = 0.0, auxError = 0.0; + + for(int i = 0; i < NEQ; i++) { + + // stores the new evaluation + _k2__[i] = rDY[i]; + _aux_tol = fabs(edos_new_euler_[i]) * reltol; + _tolerances_[i] = (abstol > _aux_tol) ? abstol : _aux_tol; + + // finds the greatest error between the steps + auxError = fabs(((dt / 2.0) * (_k1__[i] - _k2__[i])) / _tolerances_[i]); + + greatestError = (auxError > greatestError) ? auxError : greatestError; + } + + /// adapt the time step + greatestError += __tiny_; + previous_dt = dt; + + /// adapt the time step + dt = _beta_safety_ * dt * sqrt(1.0f / greatestError); + + if(dt < min_dt) { + dt = min_dt; + } + else if(dt > max_dt) { + dt = max_dt; + } + + if(time_new + dt > final_time) { + dt = final_time - time_new; + } + + // it doesn't accept the solution or accept and risk a NaN + if(greatestError >= 1.0f && dt > min_dt) { + // restore the old values to do it again + for(int i = 0; i < NEQ; i++) { + sv_local[i] = edos_old_aux_[i]; + } + + } else { + for(int i = 0; i < NEQ; i++) { + _k_aux__[i] = _k2__[i]; + _k2__[i] = _k1__[i]; + _k1__[i] = _k_aux__[i]; + } + + for(int i = 0; i < NEQ; i++) { + sv_local[i] = edos_new_euler_[i]; + } + + if(time_new + previous_dt >= final_time) { + if(final_time == time_new) { + break; + } else if(time_new < final_time) { + dt = previous_dt = final_time - time_new; + time_new += previous_dt; + break; + } + } else { + time_new += previous_dt; + } + } + } + + for(int i = 0; i < NEQ; i++) { + *((real *)((char *)sv + pitch * i) + thread_id) = sv_local[i]; + } + + DT = dt; + TIME_NEW = time_new; + PREVIOUS_DT = previous_dt; +} + +// Solving the model for each cell in the tissue matrix ni x nj +__global__ void solve_gpu(real cur_time, real dt, real *sv, real *stim_currents, uint32_t *cells_to_solve, + uint32_t num_cells_to_solve, int num_steps, size_t pitch, bool use_adpt, + real abstol, real reltol, real max_dt, real *extra_parameters) { + int threadID = blockDim.x * blockIdx.x + threadIdx.x; + int sv_id; + + int offset = 34; + int mapping = extra_parameters[threadID+offset]; + + // Each thread solves one cell model + if(threadID < num_cells_to_solve) { + if(cells_to_solve) + sv_id = cells_to_solve[threadID]; + else + sv_id = threadID; + + if(!use_adpt) { + real rDY[NEQ]; + + for(int n = 0; n < num_steps; ++n) { + + RHS_gpu(sv, rDY, stim_currents[threadID], sv_id, dt, pitch, false, extra_parameters, mapping); + + // Rush-Larsen/Euler + *((real *)((char *)sv + pitch * 0) + sv_id) = dt * rDY[0] + *((real *)((char *)sv + pitch * 0) + sv_id); + *((real *)((char *)sv + pitch * 1) + sv_id) = rDY[1]; + *((real *)((char *)sv + pitch * 2) + sv_id) = rDY[2]; + *((real *)((char *)sv + pitch * 3) + sv_id) = rDY[3]; + *((real *)((char *)sv + pitch * 4) + sv_id) = rDY[4]; + *((real *)((char *)sv + pitch * 5) + sv_id) = rDY[5]; + *((real *)((char *)sv + pitch * 6) + sv_id) = rDY[6]; + *((real *)((char *)sv + pitch * 7) + sv_id) = rDY[7]; + *((real *)((char *)sv + pitch * 8) + sv_id) = rDY[8]; + *((real *)((char *)sv + pitch * 9) + sv_id) = rDY[9]; + *((real *)((char *)sv + pitch * 10) + sv_id) = rDY[10]; + *((real *)((char *)sv + pitch * 11) + sv_id) = rDY[11]; + *((real *)((char *)sv + pitch * 12) + sv_id) = rDY[12]; + *((real *)((char *)sv + pitch * 13) + sv_id) = rDY[13]; + *((real *)((char *)sv + pitch * 14) + sv_id) = rDY[14]; + *((real *)((char *)sv + pitch * 15) + sv_id) = rDY[15]; + *((real *)((char *)sv + pitch * 16) + sv_id) = dt * rDY[16] + *((real *)((char *)sv + pitch * 16) + sv_id); + *((real *)((char *)sv + pitch * 17) + sv_id) = dt * rDY[17] + *((real *)((char *)sv + pitch * 17) + sv_id); + *((real *)((char *)sv + pitch * 18) + sv_id) = dt * rDY[18] + *((real *)((char *)sv + pitch * 18) + sv_id); + *((real *)((char *)sv + pitch * 19) + sv_id) = dt * rDY[19] + *((real *)((char *)sv + pitch * 19) + sv_id); + *((real *)((char *)sv + pitch * 20) + sv_id) = dt * rDY[20] + *((real *)((char *)sv + pitch * 20) + sv_id); + } + } else { + solve_forward_euler_gpu_adpt(sv, stim_currents[threadID], cur_time + max_dt, sv_id, pitch, abstol, reltol, dt, max_dt, extra_parameters, mapping); + } + } +} + diff --git a/src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL.h b/src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL.h new file mode 100644 index 00000000..da2d38d1 --- /dev/null +++ b/src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL.h @@ -0,0 +1,30 @@ +#ifndef MONOALG3D_MODEL_COURTEMANCHE_RAMIREZ_NATTEL_1998_VICTOR_RL_H +#define MONOALG3D_MODEL_COURTEMANCHE_RAMIREZ_NATTEL_1998_VICTOR_RL_H + +#include "../model_common.h" + +#define NEQ 21 +#define INITIAL_V (-81.180000f) + +#ifdef __CUDACC__ + +#include "../../gpu_utils/gpu_utils.h" + +inline __device__ void RHS_gpu(real *sv, real *rDY_, real stim_current, int threadID_, real dt, size_t pitch, bool use_adpt_dt, real *extra_parameters, int mapping); + +__global__ void kernel_set_model_initial_conditions(real *sv, int num_volumes, size_t pitch, bool use_adpt, real min_dt, real *extra_parameters); + +__global__ void solve_gpu(real cur_time, real dt, real *sv, real* stim_currents, + uint32_t *cells_to_solve, uint32_t num_cells_to_solve, + int num_steps, size_t pitch, bool use_adpt, + real abstol, real reltol, real max_dt, real *extra_parameters); + +#endif + +void RHS_cpu(const real *sv, real *rDY_, real stim_current, real dt, real *extra_parameters, int mapping); +inline void solve_forward_euler_cpu_adpt(real *sv, real stim_curr, real final_time, int thread_id, struct ode_solver *solver, real *extra_parameters, int mapping); + +void solve_model_ode_cpu(real dt, real *sv, real stim_current, real *extra_parameters, int mapping); + +#endif //MONOALG3D_MODEL_COURTEMANCHE_RAMIREZ_NATTEL_1998_VICTOR_RL_H + diff --git a/src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL_common.inc.c b/src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL_common.inc.c new file mode 100644 index 00000000..cca29d09 --- /dev/null +++ b/src/models_library/courtemanche_ramirez_nattel_1998/courtemanche_ramirez_nattel_1998_Victor_RL_common.inc.c @@ -0,0 +1,313 @@ + +#define IFNUMBER_1(name)if((V_old_==(-4.713000000000000e+01))) { (name) = 3.200000000000000e+00; } else{ (name) = ((3.200000000000000e-01*(V_old_+4.713000000000000e+01))/(1.000000000000000e+00-exp(((-1.000000000000000e-01)*(V_old_+4.713000000000000e+01))))); } +#define IFNUMBER_2(name)if((V_old_<(-4.000000000000000e+01))) { (name) = (1.350000000000000e-01*exp(((V_old_+8.000000000000000e+01)/(-6.800000000000000e+00)))); } else{ (name) = 0.000000000000000e+00; } +#define IFNUMBER_3(name)if((V_old_<(-4.000000000000000e+01))) { (name) = ((3.560000000000000e+00*exp((7.900000000000000e-02*V_old_)))+(3.100000000000000e+05*exp((3.500000000000000e-01*V_old_)))); } else{ (name) = (1.000000000000000e+00/(1.300000000000000e-01*(1.000000000000000e+00+exp(((V_old_+1.066000000000000e+01)/(-1.110000000000000e+01)))))); } +#define IFNUMBER_4(name)if((V_old_<(-4.000000000000000e+01))) { (name) = (((((-1.271400000000000e+05)*exp((2.444000000000000e-01*V_old_)))-(3.474000000000000e-05*exp(((-4.391000000000000e-02)*V_old_))))*(V_old_+3.778000000000000e+01))/(1.000000000000000e+00+exp((3.110000000000000e-01*(V_old_+7.923000000000000e+01))))); } else{ (name) = 0.000000000000000e+00; } +#define IFNUMBER_5(name)if((V_old_<(-4.000000000000000e+01))) { (name) = ((1.212000000000000e-01*exp(((-1.052000000000000e-02)*V_old_)))/(1.000000000000000e+00+exp(((-1.378000000000000e-01)*(V_old_+4.014000000000000e+01))))); } else{ (name) = ((3.000000000000000e-01*exp(((-2.535000000000000e-07)*V_old_)))/(1.000000000000000e+00+exp(((-1.000000000000000e-01)*(V_old_+3.200000000000000e+01))))); } +#define IFNUMBER_6(name)if((fabs((V_old_+1.410000000000000e+01))<1.000000000000000e-10)) { (name) = 1.500000000000000e-03; } else{ (name) = ((3.000000000000000e-04*(V_old_+1.410000000000000e+01))/(1.000000000000000e+00-exp(((V_old_+1.410000000000000e+01)/(-5.000000000000000e+00))))); } +#define IFNUMBER_7(name)if((fabs((V_old_-3.332800000000000e+00))<1.000000000000000e-10)) { (name) = 3.783611800000000e-04; } else{ (name) = ((7.389800000000000e-05*(V_old_-3.332800000000000e+00))/(exp(((V_old_-3.332800000000000e+00)/5.123700000000000e+00))-1.000000000000000e+00)); } +#define IFNUMBER_8(name)if((fabs((V_old_-1.990000000000000e+01))<1.000000000000000e-10)) { (name) = 6.800000000000000e-04; } else{ (name) = ((4.000000000000000e-05*(V_old_-1.990000000000000e+01))/(1.000000000000000e+00-exp(((V_old_-1.990000000000000e+01)/(-1.700000000000000e+01))))); } +#define IFNUMBER_9(name)if((fabs((V_old_-1.990000000000000e+01))<1.000000000000000e-10)) { (name) = 3.150000000000000e-04; } else{ (name) = ((3.500000000000000e-05*(V_old_-1.990000000000000e+01))/(exp(((V_old_-1.990000000000000e+01)/9.000000000000000e+00))-1.000000000000000e+00)); } +#define IFNUMBER_10(name)if((fabs((V_old_+1.000000000000000e+01))<1.000000000000000e-10)) { (name) = (4.579000000000000e+00/(1.000000000000000e+00+exp(((V_old_+1.000000000000000e+01)/(-6.240000000000000e+00))))); } else{ (name) = ((1.000000000000000e+00-exp(((V_old_+1.000000000000000e+01)/(-6.240000000000000e+00))))/(3.500000000000000e-02*(V_old_+1.000000000000000e+01)*(1.000000000000000e+00+exp(((V_old_+1.000000000000000e+01)/(-6.240000000000000e+00)))))); } +#define IFNUMBER_11(name)if((fabs((V_old_-7.900000000000000e+00))<1.000000000000000e-10)) { (name) = ((6.000000000000000e+00*2.000000000000000e-01)/1.300000000000000e+00); } else{ (name) = ((6.000000000000000e+00*(1.000000000000000e+00-exp(((-(V_old_-7.900000000000000e+00))/5.000000000000000e+00))))/((1.000000000000000e+00+(3.000000000000000e-01*exp(((-(V_old_-7.900000000000000e+00))/5.000000000000000e+00))))*1.000000000000000e+00*(V_old_-7.900000000000000e+00))); } //Parameters + +// VGM changes: +real ICaL_Block = 0.0; +real IK1_Block = 0.0; +real IKr_Block = 0.0; +real IKs_Block = 0.0; +real IKur_Block = 0.0; +real INa_Block = 0.0; +real INaK_Block = 0.0; +real INCX_Block = 0.0; +real Ito_Block = 0.0; +real IK2P_Block = 0.0; +real GKur = extra_parameters[0] * (1 - IKur_Block); +real GKr = extra_parameters[1] * (1 - IKr_Block); +real Gto = extra_parameters[2] * (1 - Ito_Block); +real GK1 = extra_parameters[3] * (1 - IK1_Block); +real GKs = extra_parameters[4] * (1 - IKs_Block); +real GCaL = extra_parameters[5] * (1 - ICaL_Block); +real GNaK = extra_parameters[6] * (1 - INaK_Block); +real GNCX = extra_parameters[7] * (1 - INCX_Block); +real GNa = extra_parameters[8] * (1 - INa_Block); +real GK2P = 0; //extra_parameters[9] * (1 - IK2P_Block); +real Gup = extra_parameters[10]; +real Grel = extra_parameters[11]; +real Gleak= extra_parameters[12]; +real GCaP = 1.0; +real GNab = 1.0; +real GCab = 1.0; +// End changes + +const real Cm = 1.000000000000000e+02f; +const real g_Na = 7.800000000000000e+00f; +const real R = 8.314299999999999e+00f; +const real T = 3.100000000000000e+02f; +const real F = 9.648670000000000e+01f; +const real Na_o = 1.400000000000000e+02f; +const real K_o = 5.400000000000000e+00f; +const real g_K1 = 9.000000000000000e-02f; +const real g_to = 1.652000000000000e-01f; +const real K_Q10 = 3.000000000000000e+00f; +const real g_Kr = 2.941176500000000e-02f; +const real g_Ks = 1.294117600000000e-01f; +const real g_Ca_L = 1.237500000000000e-01f; +const real i_NaK_max = 5.993387400000000e-01f; +const real Km_Na_i = 1.000000000000000e+01f; +const real Km_K_o = 1.500000000000000e+00f; +const real Ca_o = 1.800000000000000e+00f; +const real g_B_Na = 6.744375000000000e-04f; +const real g_B_Ca = 1.131000000000000e-03f; +const real g_B_K = 0.000000000000000e+00f; +const real I_NaCa_max = 1.600000000000000e+03f; +const real gamma = 3.500000000000000e-01f; +const real K_mNa = 8.750000000000000e+01f; +const real K_mCa = 1.380000000000000e+00f; +const real K_sat = 1.000000000000000e-01f; +const real i_CaP_max = 2.750000000000000e-01f; +const real K_rel = 3.000000000000000e+01f; +const real tau_tr = 1.800000000000000e+02f; +const real I_up_max = 5.000000000000000e-03f; +const real K_up = 9.200000000000000e-04f; +const real Ca_up_max = 1.500000000000000e+01f; +const real CMDN_max = 5.000000000000000e-02f; +const real Km_CMDN = 2.380000000000000e-03f; +const real TRPN_max = 7.000000000000001e-02f; +const real Km_TRPN = 5.000000000000000e-04f; +const real CSQN_max = 1.000000000000000e+01f; +const real Km_CSQN = 8.000000000000000e-01f; +const real V_cell = 2.010000000000000e+04f; + +// real calc_E_Na = (((R*T)/F)*log((Na_o/Na_i_old_))); //3 +// real calc_alpha_m = 0.0f; +// IFNUMBER_1(calc_alpha_m); //4 +// real calc_beta_m = (8.000000000000000e-02*exp(((-V_old_)/1.100000000000000e+01))); //5 +// real calc_alpha_h = 0.0f; +// IFNUMBER_2(calc_alpha_h); //9 +// real calc_beta_h = 0.0f; +// IFNUMBER_3(calc_beta_h); //10 +// real calc_alpha_j = 0.0f; +// IFNUMBER_4(calc_alpha_j); //14 +// real calc_beta_j = 0.0f; +// IFNUMBER_5(calc_beta_j); //15 +// real calc_E_K = (((R*T)/F)*log((K_o/K_i_old_))); //19 +// real calc_alpha_oa = (6.500000000000000e-01*pow((exp(((V_old_-(-1.000000000000000e+01))/(-8.500000000000000e+00)))+exp((((V_old_-(-1.000000000000000e+01))-4.000000000000000e+01)/(-5.900000000000000e+01)))),(-1.000000000000000e+00))); //22 +// real calc_beta_oa = (6.500000000000000e-01*pow((2.500000000000000e+00+exp((((V_old_-(-1.000000000000000e+01))+7.200000000000000e+01)/1.700000000000000e+01))),(-1.000000000000000e+00))); //23 +// real calc_oa_infinity = pow((1.000000000000000e+00+exp((((V_old_-(-1.000000000000000e+01))+1.047000000000000e+01)/(-1.754000000000000e+01)))),(-1.000000000000000e+00)); //25 +// real calc_alpha_oi = pow((1.853000000000000e+01+(1.000000000000000e+00*exp((((V_old_-(-1.000000000000000e+01))+1.037000000000000e+02)/1.095000000000000e+01)))),(-1.000000000000000e+00)); //27 +// real calc_beta_oi = pow((3.556000000000000e+01+(1.000000000000000e+00*exp((((V_old_-(-1.000000000000000e+01))-8.740000000000000e+00)/(-7.440000000000000e+00))))),(-1.000000000000000e+00)); //28 +// real calc_oi_infinity = pow((1.000000000000000e+00+exp((((V_old_-(-1.000000000000000e+01))+3.310000000000000e+01)/5.300000000000000e+00))),(-1.000000000000000e+00)); //30 +// real calc_g_Kur = (5.000000000000000e-03+(5.000000000000000e-02/(1.000000000000000e+00+exp(((V_old_-1.500000000000000e+01)/(-1.300000000000000e+01)))))); //32 +// real calc_alpha_ua = (6.500000000000000e-01*pow((exp(((V_old_-(-1.000000000000000e+01))/(-8.500000000000000e+00)))+exp((((V_old_-(-1.000000000000000e+01))-4.000000000000000e+01)/(-5.900000000000000e+01)))),(-1.000000000000000e+00))); //34 +// real calc_beta_ua = (6.500000000000000e-01*pow((2.500000000000000e+00+exp((((V_old_-(-1.000000000000000e+01))+7.200000000000000e+01)/1.700000000000000e+01))),(-1.000000000000000e+00))); //35 +// real calc_ua_infinity = pow((1.000000000000000e+00+exp((((V_old_-(-1.000000000000000e+01))+2.030000000000000e+01)/(-9.600000000000000e+00)))),(-1.000000000000000e+00)); //37 +// real calc_alpha_ui = pow((2.100000000000000e+01+(1.000000000000000e+00*exp((((V_old_-(-1.000000000000000e+01))-1.950000000000000e+02)/(-2.800000000000000e+01))))),(-1.000000000000000e+00)); //39 +// real calc_beta_ui = (1.000000000000000e+00/exp((((V_old_-(-1.000000000000000e+01))-1.680000000000000e+02)/(-1.600000000000000e+01)))); //40 +// real calc_ui_infinity = pow((1.000000000000000e+00+exp((((V_old_-(-1.000000000000000e+01))-1.094500000000000e+02)/2.748000000000000e+01))),(-1.000000000000000e+00)); //42 +// real calc_alpha_xr = 0.0f; +// IFNUMBER_6(calc_alpha_xr); //45 +// real calc_beta_xr = 0.0f; +// IFNUMBER_7(calc_beta_xr); //46 +// real calc_xr_infinity = pow((1.000000000000000e+00+exp(((V_old_+1.410000000000000e+01)/(-6.500000000000000e+00)))),(-1.000000000000000e+00)); //48 +// real calc_alpha_xs = 0.0f; +// IFNUMBER_8(calc_alpha_xs); //51 +// real calc_beta_xs = 0.0f; +// IFNUMBER_9(calc_beta_xs); //52 +// real calc_xs_infinity = pow((1.000000000000000e+00+exp(((V_old_-1.990000000000000e+01)/(-1.270000000000000e+01)))),(-5.000000000000000e-01)); //54 +// real calc_i_Ca_L = (Cm*g_Ca_L*d_old_*f_old_*f_Ca_old_*(V_old_-6.500000000000000e+01)); //56 +// real calc_d_infinity = pow((1.000000000000000e+00+exp(((V_old_+1.000000000000000e+01)/(-8.000000000000000e+00)))),(-1.000000000000000e+00)); //57 +// real calc_tau_d = 0.0f; +// IFNUMBER_10(calc_tau_d); //58 +// real calc_f_infinity = (exp(((-(V_old_+2.800000000000000e+01))/6.900000000000000e+00))/(1.000000000000000e+00+exp(((-(V_old_+2.800000000000000e+01))/6.900000000000000e+00)))); //60 +// real calc_tau_f = (9.000000000000000e+00*pow(((1.970000000000000e-02*exp(((-pow(3.370000000000000e-02,2.000000000000000e+00))*pow((V_old_+1.000000000000000e+01),2.000000000000000e+00))))+2.000000000000000e-02),(-1.000000000000000e+00))); //61 +// real calc_f_Ca_infinity = pow((1.000000000000000e+00+(Ca_i_old_/3.500000000000000e-04)),(-1.000000000000000e+00)); //63 +// real calc_tau_f_Ca = 2.000000000000000e+00; //64 +// real calc_sigma = ((1.000000000000000e+00/7.000000000000000e+00)*(exp((Na_o/6.730000000000000e+01))-1.000000000000000e+00)); //66 +// real calc_E_Ca = (((R*T)/(2.000000000000000e+00*F))*log((Ca_o/Ca_i_old_))); //69 +// real calc_i_NaCa = ((Cm*I_NaCa_max*((exp(((gamma*F*V_old_)/(R*T)))*pow(Na_i_old_,3.000000000000000e+00)*Ca_o)-(exp((((gamma-1.000000000000000e+00)*F*V_old_)/(R*T)))*pow(Na_o,3.000000000000000e+00)*Ca_i_old_)))/((pow(K_mNa,3.000000000000000e+00)+pow(Na_o,3.000000000000000e+00))*(K_mCa+Ca_o)*(1.000000000000000e+00+(K_sat*exp((((gamma-1.000000000000000e+00)*V_old_*F)/(R*T))))))); //73 +// real calc_i_CaP = ((Cm*i_CaP_max*Ca_i_old_)/(5.000000000000000e-04+Ca_i_old_)); //74 +// real calc_i_rel = (K_rel*pow(u_old_,2.000000000000000e+00)*v_old_*w_old_*(Ca_rel_old_-Ca_i_old_)); //76 +// real calc_tau_u = 8.000000000000000e+00; //77 +// real calc_tau_w = 0.0f; +// IFNUMBER_11(calc_tau_w); //83 +// real calc_w_infinity = (1.000000000000000e+00-pow((1.000000000000000e+00+exp(((-(V_old_-4.000000000000000e+01))/1.700000000000000e+01))),(-1.000000000000000e+00))); //84 +// real calc_i_tr = ((Ca_up_old_-Ca_rel_old_)/tau_tr); //86 +// real calc_i_up = (I_up_max/(1.000000000000000e+00+(K_up/Ca_i_old_))); //87 +// real calc_i_up_leak = ((I_up_max*Ca_up_old_)/Ca_up_max); //88 +// //real calc_Ca_CMDN = ((CMDN_max*Ca_i_old_)/(Ca_i_old_+Km_CMDN)); //89 +// //real calc_Ca_TRPN = ((TRPN_max*Ca_i_old_)/(Ca_i_old_+Km_TRPN)); //90 +// //real calc_Ca_CSQN = ((CSQN_max*Ca_rel_old_)/(Ca_rel_old_+Km_CSQN)); //91 +// real calc_V_i = (V_cell*6.800000000000000e-01); //92 +// real calc_V_rel = (4.800000000000000e-03*V_cell); //93 +// real calc_V_up = (5.520000000000000e-02*V_cell); //94 +// real calc_B2 = (1.000000000000000e+00+((TRPN_max*Km_TRPN)/pow((Ca_i_old_+Km_TRPN),2.000000000000000e+00))+((CMDN_max*Km_CMDN)/pow((Ca_i_old_+Km_CMDN),2.000000000000000e+00))); //99 +// real calc_m_inf = (calc_alpha_m/(calc_alpha_m+calc_beta_m)); //6 +// real calc_tau_m = (1.000000000000000e+00/(calc_alpha_m+calc_beta_m)); //7 +// real calc_h_inf = (calc_alpha_h/(calc_alpha_h+calc_beta_h)); //11 +// real calc_tau_h = (1.000000000000000e+00/(calc_alpha_h+calc_beta_h)); //12 +// real calc_j_inf = (calc_alpha_j/(calc_alpha_j+calc_beta_j)); //16 +// real calc_tau_j = (1.000000000000000e+00/(calc_alpha_j+calc_beta_j)); //17 +// real calc_i_K1 = ((Cm*g_K1*(V_old_-calc_E_K))/(1.000000000000000e+00+exp((7.000000000000001e-02*(V_old_+8.000000000000000e+01))))); //20 +// real calc_i_to = (Cm*g_to*pow(oa_old_,3.000000000000000e+00)*oi_old_*(V_old_-calc_E_K)); //21 +// real calc_tau_oa = (pow((calc_alpha_oa+calc_beta_oa),(-1.000000000000000e+00))/K_Q10); //24 +// real calc_tau_oi = (pow((calc_alpha_oi+calc_beta_oi),(-1.000000000000000e+00))/K_Q10); //29 +// real calc_i_Kur = (Cm*calc_g_Kur*pow(ua_old_,3.000000000000000e+00)*ui_old_*(V_old_-calc_E_K)); //33 +// real calc_tau_ua = (pow((calc_alpha_ua+calc_beta_ua),(-1.000000000000000e+00))/K_Q10); //36 +// real calc_tau_ui = (pow((calc_alpha_ui+calc_beta_ui),(-1.000000000000000e+00))/K_Q10); //41 +// real calc_i_Kr = ((Cm*g_Kr*xr_old_*(V_old_-calc_E_K))/(1.000000000000000e+00+exp(((V_old_+1.500000000000000e+01)/2.240000000000000e+01)))); //44 +// real calc_tau_xr = pow((calc_alpha_xr+calc_beta_xr),(-1.000000000000000e+00)); //47 +// real calc_i_Ks = (Cm*g_Ks*pow(xs_old_,2.000000000000000e+00)*(V_old_-calc_E_K)); //50 +// real calc_tau_xs = (5.000000000000000e-01*pow((calc_alpha_xs+calc_beta_xs),(-1.000000000000000e+00))); //53 +// real calc_f_NaK = pow((1.000000000000000e+00+(1.245000000000000e-01*exp((((-1.000000000000000e-01)*F*V_old_)/(R*T))))+(3.650000000000000e-02*calc_sigma*exp((((-F)*V_old_)/(R*T))))),(-1.000000000000000e+00)); //67 +// real calc_i_B_Na = (Cm*g_B_Na*(V_old_-calc_E_Na)); //70 +// real calc_i_B_Ca = (Cm*g_B_Ca*(V_old_-calc_E_Ca)); //71 +// real calc_i_B_K = (Cm*g_B_K*(V_old_-calc_E_K)); //72 +// real calc_i_Na = (Cm*g_Na*pow(m_old_,3.000000000000000e+00)*h_old_*j_old_*(V_old_-calc_E_Na)); //2 +// real calc_i_NaK = ((((Cm*i_NaK_max*calc_f_NaK*1.000000000000000e+00)/(1.000000000000000e+00+pow((Km_Na_i/Na_i_old_),1.500000000000000e+00)))*K_o)/(K_o+Km_K_o)); //68 +// real calc_Fn = (1.000000000000000e+03*((1.000000000000000e-15*calc_V_rel*calc_i_rel)-((1.000000000000000e-15/(2.000000000000000e+00*F))*((5.000000000000000e-01*calc_i_Ca_L)-(2.000000000000000e-01*calc_i_NaCa))))); //75 +// real calc_B1 = ((((2.000000000000000e+00*calc_i_NaCa)-(calc_i_CaP+calc_i_Ca_L+calc_i_B_Ca))/(2.000000000000000e+00*calc_V_i*F))+(((calc_V_up*(calc_i_up_leak-calc_i_up))+(calc_i_rel*calc_V_rel))/calc_V_i)); //98 +// real calc_u_infinity = pow((1.000000000000000e+00+exp(((-(calc_Fn-3.417500000000000e-13))/1.367000000000000e-15))),(-1.000000000000000e+00)); //78 +// real calc_tau_v = (1.910000000000000e+00+(2.090000000000000e+00*pow((1.000000000000000e+00+exp(((-(calc_Fn-3.417500000000000e-13))/1.367000000000000e-15))),(-1.000000000000000e+00)))); //80 +// real calc_v_infinity = (1.000000000000000e+00-pow((1.000000000000000e+00+exp(((-(calc_Fn-6.834999999999999e-14))/1.367000000000000e-15))),(-1.000000000000000e+00))); //81 + +real calc_E_Na = (((R*T)/F)*log((Na_o/Na_i_old_))); //3 +real calc_alpha_m = 0.0f; +IFNUMBER_1(calc_alpha_m); //4 +real calc_beta_m = (8.000000000000000e-02*exp(((-V_old_)/1.100000000000000e+01))); //5 +real calc_alpha_h = 0.0f; +IFNUMBER_2(calc_alpha_h); //9 +real calc_beta_h = 0.0f; +IFNUMBER_3(calc_beta_h); //10 +real calc_alpha_j = 0.0f; +IFNUMBER_4(calc_alpha_j); //14 +real calc_beta_j = 0.0f; +IFNUMBER_5(calc_beta_j); //15 +real calc_E_K = (((R*T)/F)*log((K_o/K_i_old_))); //19 +real calc_alpha_oa = (6.500000000000000e-01*pow((exp(((V_old_-(-1.000000000000000e+01))/(-8.500000000000000e+00)))+exp((((V_old_-(-1.000000000000000e+01))-4.000000000000000e+01)/(-5.900000000000000e+01)))),(-1.000000000000000e+00))); //22 +real calc_beta_oa = (6.500000000000000e-01*pow((2.500000000000000e+00+exp((((V_old_-(-1.000000000000000e+01))+7.200000000000000e+01)/1.700000000000000e+01))),(-1.000000000000000e+00))); //23 +real calc_oa_infinity = pow((1.000000000000000e+00+exp((((V_old_-(-1.000000000000000e+01))+1.047000000000000e+01)/(-1.754000000000000e+01)))),(-1.000000000000000e+00)); //25 +real calc_alpha_oi = pow((1.853000000000000e+01+(1.000000000000000e+00*exp((((V_old_-(-1.000000000000000e+01))+1.037000000000000e+02)/1.095000000000000e+01)))),(-1.000000000000000e+00)); //27 +real calc_beta_oi = pow((3.556000000000000e+01+(1.000000000000000e+00*exp((((V_old_-(-1.000000000000000e+01))-8.740000000000000e+00)/(-7.440000000000000e+00))))),(-1.000000000000000e+00)); //28 +real calc_oi_infinity = pow((1.000000000000000e+00+exp((((V_old_-(-1.000000000000000e+01))+3.310000000000000e+01)/5.300000000000000e+00))),(-1.000000000000000e+00)); //30 +real calc_g_Kur = (5.000000000000000e-03+(5.000000000000000e-02/(1.000000000000000e+00+exp(((V_old_-1.500000000000000e+01)/(-1.300000000000000e+01)))))); //32 +real calc_alpha_ua = (6.500000000000000e-01*pow((exp(((V_old_-(-1.000000000000000e+01))/(-8.500000000000000e+00)))+exp((((V_old_-(-1.000000000000000e+01))-4.000000000000000e+01)/(-5.900000000000000e+01)))),(-1.000000000000000e+00))); //34 +real calc_beta_ua = (6.500000000000000e-01*pow((2.500000000000000e+00+exp((((V_old_-(-1.000000000000000e+01))+7.200000000000000e+01)/1.700000000000000e+01))),(-1.000000000000000e+00))); //35 +real calc_ua_infinity = pow((1.000000000000000e+00+exp((((V_old_-(-1.000000000000000e+01))+2.030000000000000e+01)/(-9.600000000000000e+00)))),(-1.000000000000000e+00)); //37 +real calc_alpha_ui = pow((2.100000000000000e+01+(1.000000000000000e+00*exp((((V_old_-(-1.000000000000000e+01))-1.950000000000000e+02)/(-2.800000000000000e+01))))),(-1.000000000000000e+00)); //39 +real calc_beta_ui = (1.000000000000000e+00/exp((((V_old_-(-1.000000000000000e+01))-1.680000000000000e+02)/(-1.600000000000000e+01)))); //40 +real calc_ui_infinity = pow((1.000000000000000e+00+exp((((V_old_-(-1.000000000000000e+01))-1.094500000000000e+02)/2.748000000000000e+01))),(-1.000000000000000e+00)); //42 +real calc_alpha_xr = 0.0f; +IFNUMBER_6(calc_alpha_xr); //45 +real calc_beta_xr = 0.0f; +IFNUMBER_7(calc_beta_xr); //46 +real calc_xr_infinity = pow((1.000000000000000e+00+exp(((V_old_+1.410000000000000e+01)/(-6.500000000000000e+00)))),(-1.000000000000000e+00)); //48 +real calc_alpha_xs = 0.0f; +IFNUMBER_8(calc_alpha_xs); //51 +real calc_beta_xs = 0.0f; +IFNUMBER_9(calc_beta_xs); //52 +real calc_xs_infinity = pow((1.000000000000000e+00+exp(((V_old_-1.990000000000000e+01)/(-1.270000000000000e+01)))),(-5.000000000000000e-01)); //54 +real calc_i_Ca_L = (Cm*GCaL*g_Ca_L*d_old_*f_old_*f_Ca_old_*(V_old_-6.500000000000000e+01)); //56 +real calc_d_infinity = pow((1.000000000000000e+00+exp(((V_old_+1.000000000000000e+01)/(-8.000000000000000e+00)))),(-1.000000000000000e+00)); //57 +real calc_tau_d = 0.0f; +IFNUMBER_10(calc_tau_d); //58 +real calc_f_infinity = (exp(((-(V_old_+2.800000000000000e+01))/6.900000000000000e+00))/(1.000000000000000e+00+exp(((-(V_old_+2.800000000000000e+01))/6.900000000000000e+00)))); //60 +real calc_tau_f = (9.000000000000000e+00*pow(((1.970000000000000e-02*exp(((-pow(3.370000000000000e-02,2.000000000000000e+00))*pow((V_old_+1.000000000000000e+01),2.000000000000000e+00))))+2.000000000000000e-02),(-1.000000000000000e+00))); //61 +real calc_f_Ca_infinity = pow((1.000000000000000e+00+(Ca_i_old_/3.500000000000000e-04)),(-1.000000000000000e+00)); //63 +real calc_tau_f_Ca = 2.000000000000000e+00; //64 +real calc_sigma = ((1.000000000000000e+00/7.000000000000000e+00)*(exp((Na_o/6.730000000000000e+01))-1.000000000000000e+00)); //66 +real calc_E_Ca = (((R*T)/(2.000000000000000e+00*F))*log((Ca_o/Ca_i_old_))); //69 +real calc_i_NaCa = ((Cm*GNCX*I_NaCa_max*((exp(((gamma*F*V_old_)/(R*T)))*pow(Na_i_old_,3.000000000000000e+00)*Ca_o)-(exp((((gamma-1.000000000000000e+00)*F*V_old_)/(R*T)))*pow(Na_o,3.000000000000000e+00)*Ca_i_old_)))/((pow(K_mNa,3.000000000000000e+00)+pow(Na_o,3.000000000000000e+00))*(K_mCa+Ca_o)*(1.000000000000000e+00+(K_sat*exp((((gamma-1.000000000000000e+00)*V_old_*F)/(R*T))))))); //73 +real calc_i_CaP = ((Cm*GCaP*i_CaP_max*Ca_i_old_)/(5.000000000000000e-04+Ca_i_old_)); //74 +real calc_i_rel = (Grel*K_rel*pow(u_old_,2.000000000000000e+00)*v_old_*w_old_*(Ca_rel_old_-Ca_i_old_)); //76 +real calc_tau_u = 8.000000000000000e+00; //77 +real calc_tau_w = 0.0f; +IFNUMBER_11(calc_tau_w); //83 +real calc_w_infinity = (1.000000000000000e+00-pow((1.000000000000000e+00+exp(((-(V_old_-4.000000000000000e+01))/1.700000000000000e+01))),(-1.000000000000000e+00))); //84 +real calc_i_tr = ((Ca_up_old_-Ca_rel_old_)/tau_tr); //86 +real calc_i_up = (Gup*I_up_max/(1.000000000000000e+00+(K_up/Ca_i_old_))); //87 +real calc_i_up_leak = (Gleak*(I_up_max*Ca_up_old_)/Ca_up_max); //88 +//real calc_Ca_CMDN = ((CMDN_max*Ca_i_old_)/(Ca_i_old_+Km_CMDN)); //89 +//real calc_Ca_TRPN = ((TRPN_max*Ca_i_old_)/(Ca_i_old_+Km_TRPN)); //90 +//real calc_Ca_CSQN = ((CSQN_max*Ca_rel_old_)/(Ca_rel_old_+Km_CSQN)); //91 +real calc_V_i = (V_cell*6.800000000000000e-01); //92 +real calc_V_rel = (4.800000000000000e-03*V_cell); //93 +real calc_V_up = (5.520000000000000e-02*V_cell); //94 +real calc_B2 = (1.000000000000000e+00+((TRPN_max*Km_TRPN)/pow((Ca_i_old_+Km_TRPN),2.000000000000000e+00))+((CMDN_max*Km_CMDN)/pow((Ca_i_old_+Km_CMDN),2.000000000000000e+00))); //99 +real calc_m_inf = (calc_alpha_m/(calc_alpha_m+calc_beta_m)); //6 +real calc_tau_m = (1.000000000000000e+00/(calc_alpha_m+calc_beta_m)); //7 +real calc_h_inf = (calc_alpha_h/(calc_alpha_h+calc_beta_h)); //11 +real calc_tau_h = (1.000000000000000e+00/(calc_alpha_h+calc_beta_h)); //12 +real calc_j_inf = (calc_alpha_j/(calc_alpha_j+calc_beta_j)); //16 +real calc_tau_j = (1.000000000000000e+00/(calc_alpha_j+calc_beta_j)); //17 +real calc_i_K1 = ((Cm*GK1*g_K1*(V_old_-calc_E_K))/(1.000000000000000e+00+exp((7.000000000000001e-02*(V_old_+8.000000000000000e+01))))); //20 +real calc_i_to = (Cm*Gto*g_to*pow(oa_old_,3.000000000000000e+00)*oi_old_*(V_old_-calc_E_K)); //21 +real calc_tau_oa = (pow((calc_alpha_oa+calc_beta_oa),(-1.000000000000000e+00))/K_Q10); //24 +real calc_tau_oi = (pow((calc_alpha_oi+calc_beta_oi),(-1.000000000000000e+00))/K_Q10); //29 +real calc_i_Kur = (Cm*GKur*calc_g_Kur*pow(ua_old_,3.000000000000000e+00)*ui_old_*(V_old_-calc_E_K)); //33 +real calc_tau_ua = (pow((calc_alpha_ua+calc_beta_ua),(-1.000000000000000e+00))/K_Q10); //36 +real calc_tau_ui = (pow((calc_alpha_ui+calc_beta_ui),(-1.000000000000000e+00))/K_Q10); //41 +real calc_i_Kr = ((Cm*GKr*g_Kr*xr_old_*(V_old_-calc_E_K))/(1.000000000000000e+00+exp(((V_old_+1.500000000000000e+01)/2.240000000000000e+01)))); //44 +real calc_tau_xr = pow((calc_alpha_xr+calc_beta_xr),(-1.000000000000000e+00)); //47 +real calc_i_Ks = (Cm*GKs*g_Ks*pow(xs_old_,2.000000000000000e+00)*(V_old_-calc_E_K)); //50 +real calc_tau_xs = (5.000000000000000e-01*pow((calc_alpha_xs+calc_beta_xs),(-1.000000000000000e+00))); //53 +real calc_f_NaK = pow((1.000000000000000e+00+(1.245000000000000e-01*exp((((-1.000000000000000e-01)*F*V_old_)/(R*T))))+(3.650000000000000e-02*calc_sigma*exp((((-F)*V_old_)/(R*T))))),(-1.000000000000000e+00)); //67 +real calc_i_B_Na = (Cm*GNab*g_B_Na*(V_old_-calc_E_Na)); //70 +real calc_i_B_Ca = (Cm*GCab*g_B_Ca*(V_old_-calc_E_Ca)); //71 +real calc_i_B_K = (Cm*g_B_K*(V_old_-calc_E_K)); //72 +real calc_i_Na = (Cm*GNa*g_Na*pow(m_old_,3.000000000000000e+00)*h_old_*j_old_*(V_old_-calc_E_Na)); //2 +real calc_i_NaK = ((((Cm*GNaK*i_NaK_max*calc_f_NaK*1.000000000000000e+00)/(1.000000000000000e+00+pow((Km_Na_i/Na_i_old_),1.500000000000000e+00)))*K_o)/(K_o+Km_K_o)); //68 +real calc_Fn = (1.000000000000000e+03*((1.000000000000000e-15*calc_V_rel*calc_i_rel)-((1.000000000000000e-15/(2.000000000000000e+00*F))*((5.000000000000000e-01*calc_i_Ca_L)-(2.000000000000000e-01*calc_i_NaCa))))); //75 +real calc_B1 = ((((2.000000000000000e+00*calc_i_NaCa)-(calc_i_CaP+calc_i_Ca_L+calc_i_B_Ca))/(2.000000000000000e+00*calc_V_i*F))+(((calc_V_up*(calc_i_up_leak-calc_i_up))+(calc_i_rel*calc_V_rel))/calc_V_i)); //98 +real calc_u_infinity = pow((1.000000000000000e+00+exp(((-(calc_Fn-3.417500000000000e-13))/1.367000000000000e-15))),(-1.000000000000000e+00)); //78 +real calc_tau_v = (1.910000000000000e+00+(2.090000000000000e+00*pow((1.000000000000000e+00+exp(((-(calc_Fn-3.417500000000000e-13))/1.367000000000000e-15))),(-1.000000000000000e+00)))); //80 +real calc_v_infinity = (1.000000000000000e+00-pow((1.000000000000000e+00+exp(((-(calc_Fn-6.834999999999999e-14))/1.367000000000000e-15))),(-1.000000000000000e+00))); //81 + + +// Full Euler +//rDY[0] = ((-(calc_i_Na+calc_i_K1+calc_i_to+calc_i_Kur+calc_i_Kr+calc_i_Ks+calc_i_B_Na+calc_i_B_Ca+calc_i_NaK+calc_i_CaP+calc_i_NaCa+calc_i_Ca_L+stim_current))/Cm); +//rDY[1] = ((calc_m_inf-m_old_)/calc_tau_m); +//rDY[2] = ((calc_h_inf-h_old_)/calc_tau_h); +//rDY[3] = ((calc_j_inf-j_old_)/calc_tau_j); +//rDY[4] = ((calc_oa_infinity-oa_old_)/calc_tau_oa); +//rDY[5] = ((calc_oi_infinity-oi_old_)/calc_tau_oi); +//rDY[6] = ((calc_ua_infinity-ua_old_)/calc_tau_ua); +//rDY[7] = ((calc_ui_infinity-ui_old_)/calc_tau_ui); +//rDY[8] = ((calc_xr_infinity-xr_old_)/calc_tau_xr); +//rDY[9] = ((calc_xs_infinity-xs_old_)/calc_tau_xs); +//rDY[10] = ((calc_d_infinity-d_old_)/calc_tau_d); +//rDY[11] = ((calc_f_infinity-f_old_)/calc_tau_f); +//rDY[12] = ((calc_f_Ca_infinity-f_Ca_old_)/calc_tau_f_Ca); +//rDY[13] = ((calc_u_infinity-u_old_)/calc_tau_u); +//rDY[14] = ((calc_v_infinity-v_old_)/calc_tau_v); +//rDY[15] = ((calc_w_infinity-w_old_)/calc_tau_w); +//rDY[16] = ((((-3.000000000000000e+00)*calc_i_NaK)-((3.000000000000000e+00*calc_i_NaCa)+calc_i_B_Na+calc_i_Na))/(calc_V_i*F)); +//rDY[17] = (((2.000000000000000e+00*calc_i_NaK)-(calc_i_K1+calc_i_to+calc_i_Kur+calc_i_Kr+calc_i_Ks+calc_i_B_K))/(calc_V_i*F)); +//rDY[18] = (calc_B1/calc_B2); +//rDY[19] = (calc_i_up-(calc_i_up_leak+((calc_i_tr*calc_V_rel)/calc_V_up))); +//rDY[20] = ((calc_i_tr-calc_i_rel)*pow((1.000000000000000e+00+((CSQN_max*Km_CSQN)/pow((Ca_rel_old_+Km_CSQN),2.000000000000000e+00))),(-1.000000000000000e+00))); + +// Euler + Rush-Larsen +rDY[0] = ((-(calc_i_Na+calc_i_K1+calc_i_to+calc_i_Kur+calc_i_Kr+calc_i_Ks+calc_i_B_Na+calc_i_B_Ca+calc_i_NaK+calc_i_CaP+calc_i_NaCa+calc_i_Ca_L+stim_current))/Cm); + +rDY[1] = calc_m_inf + (m_old_-calc_m_inf)*exp(-dt/calc_tau_m); +rDY[2] = calc_h_inf + (h_old_-calc_h_inf)*exp(-dt/calc_tau_h); +rDY[3] = calc_j_inf + (j_old_-calc_j_inf)*exp(-dt/calc_tau_j); +rDY[4] = calc_oa_infinity + (oa_old_-calc_oa_infinity)*exp(-dt/calc_tau_oa); +rDY[5] = calc_oi_infinity + (oi_old_-calc_oi_infinity)*exp(-dt/calc_tau_oi); +rDY[6] = calc_ua_infinity + (ua_old_-calc_ua_infinity)*exp(-dt/calc_tau_ua); +rDY[7] = calc_ui_infinity + (ui_old_-calc_ui_infinity)*exp(-dt/calc_tau_ui); +rDY[8] = calc_xr_infinity + (xr_old_-calc_xr_infinity)*exp(-dt/calc_tau_xr); +rDY[9] = calc_xs_infinity + (xs_old_-calc_xs_infinity)*exp(-dt/calc_tau_xs); +rDY[10] = calc_d_infinity + (d_old_-calc_d_infinity)*exp(-dt/calc_tau_d); +rDY[11] = calc_f_infinity + (f_old_-calc_f_infinity)*exp(-dt/calc_tau_f); +rDY[12] = calc_f_Ca_infinity + (f_Ca_old_-calc_f_Ca_infinity)*exp(-dt/calc_tau_f_Ca); +rDY[13] = calc_u_infinity + (u_old_-calc_u_infinity)*exp(-dt/calc_tau_u); +rDY[14] = calc_v_infinity + (v_old_-calc_v_infinity)*exp(-dt/calc_tau_v); +rDY[15] = calc_w_infinity + (w_old_-calc_w_infinity)*exp(-dt/calc_tau_w); + +rDY[16] = ((((-3.000000000000000e+00)*calc_i_NaK)-((3.000000000000000e+00*calc_i_NaCa)+calc_i_B_Na+calc_i_Na))/(calc_V_i*F)); +rDY[17] = (((2.000000000000000e+00*calc_i_NaK)-(calc_i_K1+calc_i_to+calc_i_Kur+calc_i_Kr+calc_i_Ks+calc_i_B_K))/(calc_V_i*F)); +rDY[18] = (calc_B1/calc_B2); +rDY[19] = (calc_i_up-(calc_i_up_leak+((calc_i_tr*calc_V_rel)/calc_V_up))); +rDY[20] = ((calc_i_tr-calc_i_rel)*pow((1.000000000000000e+00+((CSQN_max*Km_CSQN)/pow((Ca_rel_old_+Km_CSQN),2.000000000000000e+00))),(-1.000000000000000e+00))); diff --git a/src/monodomain/monodomain_solver.c b/src/monodomain/monodomain_solver.c index b900ebee..282286fe 100644 --- a/src/monodomain/monodomain_solver.c +++ b/src/monodomain/monodomain_solver.c @@ -1537,6 +1537,7 @@ void write_pmj_delay (struct grid *the_grid, struct config *config, struct termi //fprintf(output_file,"====================== PULSE %u ======================\n", k+1); for(uint32_t i = 0; i < num_terminals; i++) { + uint32_t term_id = i; bool is_terminal_active = the_terminals[i].active; // [PURKINJE] Get the informaion from the Purkinje cell @@ -1566,7 +1567,9 @@ void write_pmj_delay (struct grid *the_grid, struct config *config, struct termi // Calculate the mean LAT of the tissue cells surrounding the Purkinje cell real_cpu mean_tissue_lat = 0.0; + real_cpu min_tissue_lat = __DBL_MAX__; uint32_t cur_pulse = k; + uint32_t min_tissue_lat_id = 0; for(uint32_t j = 0; j < number_tissue_cells; j++) { cell_coordinates.x = tissue_cells[j]->center.x; @@ -1588,8 +1591,12 @@ void write_pmj_delay (struct grid *the_grid, struct config *config, struct termi cur_pulse = 0; return; } - //fprintf(output_file,"\tTissue cell %u --> LAT = %g\n", j, activation_times_array_tissue[cur_pulse]); mean_tissue_lat += activation_times_array_tissue[cur_pulse]; + if (activation_times_array_tissue[cur_pulse] < min_tissue_lat) { + min_tissue_lat = activation_times_array_tissue[cur_pulse]; + min_tissue_lat_id = tissue_cells[j]->sv_position; + } + } if (is_terminal_active) { @@ -1597,8 +1604,11 @@ void write_pmj_delay (struct grid *the_grid, struct config *config, struct termi real_cpu pmj_delay = (mean_tissue_lat - purkinje_lat); - // pulse_id, terminal_id, purkinje_lat, mean_tissue_lat, pmj_delay, is_active - fprintf(output_file,"%d,%u,%g,%g,%g,%d\n", k, i, purkinje_lat, mean_tissue_lat, pmj_delay, (int)is_terminal_active); + // version 1 = pulse_id, terminal_id, purkinje_lat, mean_tissue_lat, pmj_delay, is_active + //fprintf(output_file,"%d,%u,%g,%g,%g,%d\n", k, term_id, purkinje_lat, mean_tissue_lat, pmj_delay, (int)is_terminal_active); + + // version 2 = pulse_id, purkinje_terminal_id, min_tissue_coupled_cell_id, purkinje_lat, min_tissue_lat, pmj_delay, is_active + fprintf(output_file,"%d,%u,%u,%g,%g,%g,%d\n", k, term_id, min_tissue_lat_id, purkinje_lat, min_tissue_lat, pmj_delay, (int)is_terminal_active); //log_info("[purkinje_coupling] Terminal %u (%g,%g,%g) [Pulse %d] -- Purkinje LAT = %g ms -- Tissue mean LAT = %g ms -- PMJ delay = %g ms [Active = %d]\n", i, // purkinje_cells[purkinje_index]->center.x, purkinje_cells[purkinje_index]->center.y, purkinje_cells[purkinje_index]->center.z, k, diff --git a/src/save_mesh_library/save_mesh_purkinje.c b/src/save_mesh_library/save_mesh_purkinje.c index c067131a..07089023 100644 --- a/src/save_mesh_library/save_mesh_purkinje.c +++ b/src/save_mesh_library/save_mesh_purkinje.c @@ -624,14 +624,20 @@ SAVE_MESH(save_multiple_cell_state_variables) { for (uint32_t i = 0; i < params->num_tissue_cells; i++) { if(params->tissue_cell_sv_positions[i] == -1) { if(!the_grid->adaptive) { - FOR_EACH_CELL(the_grid) { - if(cell->center.x == params->tissue_cell_centers[i*3] && cell->center.y == params->tissue_cell_centers[i*3+1] && cell->center.z == params->tissue_cell_centers[i*3+2]) { - params->tissue_cell_sv_positions[i] = cell->sv_position; - break; + FOR_EACH_CELL(the_grid) { + if (cell->active) { + if(cell->center.x == params->tissue_cell_centers[i*3] && cell->center.y == params->tissue_cell_centers[i*3+1] && cell->center.z == params->tissue_cell_centers[i*3+2]) { + params->tissue_cell_sv_positions[i] = cell->sv_position; + break; + } } } } } + printf("Found cell %u with coordinates: (%lf %lf %lf)\n", params->tissue_cell_sv_positions[i], \ + params->tissue_cell_centers[i*3], \ + params->tissue_cell_centers[i*3+1], \ + params->tissue_cell_centers[i*3+2]); } if(ode_solver->gpu) {