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/example_configs/cable_mesh_with_ToRORd_Land_endo_mid_epi.ini b/example_configs/cable_mesh_with_ToRORd_Land_endo_mid_epi.ini new file mode 100644 index 00000000..5671a78e --- /dev/null +++ b/example_configs/cable_mesh_with_ToRORd_Land_endo_mid_epi.ini @@ -0,0 +1,79 @@ +# =========================================================================== +# Author: Lucas Berg +# Description: Simple simulation to test the ToRORd_Land model in a cable. +# When no [extra_data] information is provided all cells are +# considered to be control ENDO. +# ============================================================================ +[main] +num_threads=6 +dt_pde=0.01 +simulation_time=1000.0 +abort_on_no_activity=false +use_adaptivity=false + +[update_monodomain] +main_function=update_monodomain_default + +[save_result] +print_rate=100 +output_dir=./outputs/cable_ToRORd_Land_endo_baseline_cpu_en +add_timestamp=false +binary=true +main_function=save_as_ensight + +;[save_result] +;print_rate=1 +;output_dir=./outputs/cable_ToRORd_Land_endo_baseline_cpu +;file_name=./outputs/cable_ToRORd_Land_endo_baseline_cpu/full_trace_first_cell_endo.txt +;output_dir=./outputs/cable_ToRORd_Land_mid_baseline_cpu +;file_name=./outputs/cable_ToRORd_Land_mid_baseline_cpu/full_trace_first_cell_mid.txt +;output_dir=./outputs/cable_ToRORd_Land_epi_baseline_cpu +;file_name=./outputs/cable_ToRORd_Land_epi_baseline_cpu/full_trace_first_cell_epi.txt +;main_function=save_one_cell_state_variables +;init_function=init_save_one_cell_state_variables +;end_function=end_save_one_cell_state_variables +;cell_center_x=50 +;cell_center_y=50 +;cell_center_z=50 +;remove_older_simulation=true + +[assembly_matrix] +init_function=set_initial_conditions_fvm +sigma_x=0.000176 +sigma_y=0.000176 +sigma_z=0.000176 +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=yes +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=20000.0 +main_function=initialize_grid_with_cable_mesh + +[ode_solver] +adaptive=false +dt=0.01 +use_gpu=no +gpu_id=0 +library_file= shared_libs/libToRORd_Land_mixed_endo_mid_epi.so + +[stim_plain] +start = 0.0 +duration = 1.0 +;period = 800.0 +current = -53.0 +x_limit = 500.0 +main_function = stim_if_x_less_than diff --git a/example_configs/cable_mesh_with_ToRORd_fkatp_endo_mid_epi.ini b/example_configs/cable_mesh_with_ToRORd_fkatp_endo_mid_epi.ini index cdc9713f..11100129 100644 --- a/example_configs/cable_mesh_with_ToRORd_fkatp_endo_mid_epi.ini +++ b/example_configs/cable_mesh_with_ToRORd_fkatp_endo_mid_epi.ini @@ -1,26 +1,27 @@ - # ==================================================================== # Author: Lucas Berg -# Description: +# Description: Simple simulation to test the ToRORd_fkatp model in a cable. +# When no [extra_data] information is provided all cells are +# considered to be control ENDO. # ==================================================================== [main] -num_threads=4 -dt_pde=0.02 -simulation_time=600.0 +num_threads=6 +dt_pde=0.01 +simulation_time=500.0 abort_on_no_activity=false use_adaptivity=false [update_monodomain] main_function=update_monodomain_default +; Save using ensight format and store the state variables as CELL_DATA [save_result] -print_rate=50 -output_dir=./outputs/cable_ToRORd_endo_mid_epi_cpu -main_function=save_as_vtu -init_function=init_save_as_vtk_or_vtu -end_function=end_save_as_vtk_or_vtu -save_pvd=true -file_prefix=V +print_rate=100 +output_dir=./outputs/cable_ToRORd_fkatp_endo_baseline_cpu_en +add_timestamp=false +binary=true +save_ode_state_variables=true +main_function=save_as_ensight #[save_result] #print_rate=1 @@ -44,7 +45,7 @@ main_function=homogeneous_sigma_assembly_matrix [linear_system_solver] tolerance=1e-16 use_preconditioner=no -max_iterations=200 +max_iterations=500 library_file=shared_libs/libdefault_linear_system_solver.so use_gpu=no main_function=conjugate_gradient @@ -60,19 +61,19 @@ cable_length=20000.0 main_function=initialize_grid_with_cable_mesh [ode_solver] -adaptive=true -dt=0.0001 +adaptive=false +dt=0.01 use_gpu=no gpu_id=0 library_file= shared_libs/libToRORd_fkatp_mixed_endo_mid_epi.so [stim_plain] start = 0.0 -duration = 2.0 -period = 1000.0 +duration = 1.0 +;period = 1000.0 current = -53.0 x_limit = 500.0 main_function=stim_if_x_less_than [extra_data] -main_function=set_extra_data_mixed_model_epi_mid_endo +main_function=set_extra_data_mixed_torord_fkatp_epi_mid_endo diff --git a/example_configs/plain_mesh_ToRORd_Land_endo_mid_epi.ini b/example_configs/plain_mesh_ToRORd_Land_endo_mid_epi.ini new file mode 100644 index 00000000..0ae22063 --- /dev/null +++ b/example_configs/plain_mesh_ToRORd_Land_endo_mid_epi.ini @@ -0,0 +1,66 @@ +[main] +num_threads=6 +dt_pde=0.01 +simulation_time=500.0 +abort_on_no_activity=false +use_adaptivity=false + +[update_monodomain] +main_function=update_monodomain_default + +[save_result] +print_rate=100 +output_dir=./outputs/plain_ToRORd_Land_endo_mid_epi_gpu_en +add_timestamp=false +binary=true +save_ode_state_variables=false +main_function=save_as_ensight + +[assembly_matrix] +init_function=set_initial_conditions_fvm +sigma_x=0.000176 +sigma_y=0.000176 +sigma_z=0.000176 +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=yes +main_function=conjugate_gradient +init_function=init_conjugate_gradient +end_function=end_conjugate_gradient + +[domain] +name=Plain Mesh +num_layers=1 +start_dx=100.0 +start_dy=100.0 +start_dz=100.0 +side_length=20000 +main_function=initialize_grid_with_square_mesh + +[ode_solver] +adaptive=false +;reltol=1e-10 +;abstol=1e-10 +dt=0.01 +use_gpu=yes +gpu_id=0 +library_file= shared_libs/libToRORd_Land_mixed_endo_mid_epi.so + +[stim_corner] +start = 0.0 +duration = 1.0 +current = -53.0 +min_x = 0.0 +max_x = 1000.0 +min_y = 0.0 +max_y = 1000.0 +main_function=stim_x_y_limits + +[extra_data] +main_function=set_extra_data_mixed_torord_Land_epi_mid_endo diff --git a/example_configs/plain_mesh_ToRORd_fkatp_endo_mid_epi.ini b/example_configs/plain_mesh_ToRORd_fkatp_endo_mid_epi.ini index 0c7fe066..989e3edf 100644 --- a/example_configs/plain_mesh_ToRORd_fkatp_endo_mid_epi.ini +++ b/example_configs/plain_mesh_ToRORd_fkatp_endo_mid_epi.ini @@ -1,25 +1,16 @@ [main] num_threads=6 -dt_pde=0.02 -simulation_time=600.0 +dt_pde=0.01 +simulation_time=500.0 abort_on_no_activity=false use_adaptivity=false [update_monodomain] main_function=update_monodomain_default -#[save_result] -#print_rate=50 -#output_dir=./outputs/plain_ToRORd_endo_mid_epi -#main_function=save_as_vtu -#init_function=init_save_as_vtk_or_vtu -#end_function=end_save_as_vtk_or_vtu -#save_pvd=true -#file_prefix=V - [save_result] -print_rate=50 -output_dir=./outputs/plain_ToRORd_endo_mid_epi_ensight +print_rate=100 +output_dir=./outputs/plain_ToRORd_fkatp_endo_mid_epi_gpu_en add_timestamp=false binary=true save_ode_state_variables=false @@ -36,7 +27,7 @@ main_function=homogeneous_sigma_assembly_matrix [linear_system_solver] tolerance=1e-16 use_preconditioner=no -max_iterations=200 +max_iterations=500 library_file=shared_libs/libdefault_linear_system_solver.so use_gpu=yes main_function=conjugate_gradient @@ -56,9 +47,6 @@ main_function=initialize_grid_with_square_mesh dt=0.01 use_gpu=yes gpu_id=0 -;adaptive=false -;abstol=1e-12 -;reltol=1e-5 library_file= shared_libs/libToRORd_fkatp_mixed_endo_mid_epi.so [stim_corner] diff --git a/example_configs/trace_state_vector_cable_torord.ini b/example_configs/trace_state_vector_cable_torord.ini deleted file mode 100644 index 037c654b..00000000 --- a/example_configs/trace_state_vector_cable_torord.ini +++ /dev/null @@ -1,71 +0,0 @@ -# ==================================================================== -# Author: Lucas Berg -# Description: This example demonstrate how to save the state-vector -# trace of any given cellular model. -# ==================================================================== - -[main] -num_threads=4 -dt_pde=0.02 -simulation_time=40000 -abort_on_no_activity=false -use_adaptivity=false -#quiet=true - -[update_monodomain] -main_function=update_monodomain_default - -[save_result] -print_rate=100 -output_dir=./outputs/trace_sv_torord -file_name=./outputs/trace_sv_torord/trace_state_vector_torord.dat -cell_center_x=5050 -cell_center_y=50 -cell_center_z=50 -file_prefix=V -library_file=shared_libs/libdefault_save_mesh_purkinje.so -main_function=save_one_cell_state_variables -init_function=init_save_one_cell_state_variables -end_function=end_save_one_cell_state_variables -remove_older_simulation=true - -[assembly_matrix] -sigma_x=0.00005336 -sigma_y=0.00005336 -sigma_z=0.00005336 -library_file=shared_libs/libdefault_matrix_assembly.so -main_function=homogeneous_sigma_assembly_matrix -init_function=set_initial_conditions_fvm - -[linear_system_solver] -tolerance=1e-16 -use_preconditioner=no -use_gpu=yes -max_iterations=200 -library_file=shared_libs/libdefault_linear_system_solver.so -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=true -dt=0.0001 -use_gpu=yes -gpu_id=0 -library_file=shared_libs/libToRORd_fkatp_endo_2019.so - -[stim_plain] -start = 0.0 -duration = 1.0 -current = -53.0 -period = 1000.0 -x_limit = 500.0 -main_function = stim_if_x_less_than diff --git a/example_configs/trace_state_vector_cell_ToRORd_Land.ini b/example_configs/trace_state_vector_cell_ToRORd_Land.ini new file mode 100644 index 00000000..4530ce24 --- /dev/null +++ b/example_configs/trace_state_vector_cell_ToRORd_Land.ini @@ -0,0 +1,148 @@ +# ==================================================================== +# Author: Lucas Berg +# Description: Simple simulation to test the ToRORd_Land model. +# - Trace the state-vector of a single cell. +# - 20 pulses at BCL=800ms +# - Model solved using fixed timestep Euler/Rush-Larsen +# - Apply different configurations to the ToRORd_Land model using +# current modifiers in the [extra_data] section +# ==================================================================== +[main] +num_threads=1 +dt_pde=0.01 +simulation_time=16000 +abort_on_no_activity=false +use_adaptivity=false + +[update_monodomain] +main_function=update_monodomain_default + +;[save_result] +;print_rate=100 +;output_dir=./outputs/cable_ToRORd_Land_endo_baseline_cpu_en +;add_timestamp=false +;binary=true +;main_function=save_as_ensight + +[save_result] +print_rate=1 +output_dir=./outputs/cell_ToRORd_Land_endo_baseline_fixed_dt_0,01_cpu +file_name=./outputs/cell_ToRORd_Land_endo_baseline_fixed_dt_0,01_cpu/full_trace_first_cell_endo.txt +;output_dir=./outputs/cell_ToRORd_Land_mid_baseline_fixed_dt_0,01_cpu +;file_name=./outputs/cell_ToRORd_Land_mid_baseline_fixed_dt_0,01_cpu/full_trace_first_cell_mid.txt +;output_dir=./outputs/cell_ToRORd_Land_epi_baseline_fixed_dt_0,01_cpu +;file_name=./outputs/cell_ToRORd_Land_epi_baseline_fixed_dt_0,01_cpu/full_trace_first_cell_epi.txt +;output_dir=./outputs/cell_ToRORd_Land_endo_bz1_fixed_dt_0,01_cpu +;file_name=./outputs/cell_ToRORd_Land_endo_bz1_fixed_dt_0,01_cpu/full_trace_first_cell_endo.txt +;output_dir=./outputs/cell_ToRORd_Land_endo_bz2_fixed_dt_0,01_cpu +;file_name=./outputs/cell_ToRORd_Land_endo_bz2_fixed_dt_0,01_cpu/full_trace_first_cell_endo.txt +;output_dir=./outputs/cell_ToRORd_Land_endo_bz3_fixed_dt_0,01_cpu +;file_name=./outputs/cell_ToRORd_Land_endo_bz3_fixed_dt_0,01_cpu/full_trace_first_cell_endo.txt +;output_dir=./outputs/cell_ToRORd_Land_mid_bz1_fixed_dt_0,01_cpu +;file_name=./outputs/cell_ToRORd_Land_mid_bz1_fixed_dt_0,01_cpu/full_trace_first_cell_mid.txt +;output_dir=./outputs/cell_ToRORd_Land_mid_bz2_fixed_dt_0,01_cpu +;file_name=./outputs/cell_ToRORd_Land_mid_bz2_fixed_dt_0,01_cpu/full_trace_first_cell_mid.txt +;output_dir=./outputs/cell_ToRORd_Land_mid_bz3_fixed_dt_0,01_cpu +;file_name=./outputs/cell_ToRORd_Land_mid_bz3_fixed_dt_0,01_cpu/full_trace_first_cell_mid.txt +;output_dir=./outputs/cell_ToRORd_Land_epi_bz1_fixed_dt_0,01_cpu +;file_name=./outputs/cell_ToRORd_Land_epi_bz1_fixed_dt_0,01_cpu/full_trace_first_cell_epi.txt +;output_dir=./outputs/cell_ToRORd_Land_epi_bz2_fixed_dt_0,01_cpu +;file_name=./outputs/cell_ToRORd_Land_epi_bz2_fixed_dt_0,01_cpu/full_trace_first_cell_epi.txt +;output_dir=./outputs/cell_ToRORd_Land_epi_bz3_fixed_dt_0,01_cpu +;file_name=./outputs/cell_ToRORd_Land_epi_bz3_fixed_dt_0,01_cpu/full_trace_first_cell_epi.txt +main_function=save_one_cell_state_variables +init_function=init_save_one_cell_state_variables +end_function=end_save_one_cell_state_variables +cell_center_x=50 +cell_center_y=50 +cell_center_z=50 +remove_older_simulation=true + +[assembly_matrix] +init_function=set_initial_conditions_fvm +sigma_x=0.000176 +sigma_y=0.000176 +sigma_z=0.000176 +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=100.0 +main_function=initialize_grid_with_cable_mesh + +[ode_solver] +adaptive=false +;reltol=1e-10 +;abstol=1e-10 +dt=0.01 +use_gpu=no +gpu_id=0 +library_file= shared_libs/libToRORd_Land_mixed_endo_mid_epi.so + +; Solving using Euler/Rush-Larsen +[ode_solver] +adaptive=false +dt=0.01 +use_gpu=no +gpu_id=0 +library_file= shared_libs/libToRORd_Land_mixed_endo_mid_epi.so + +; Solving using Adaptive Euler +;[ode_solver] +;adaptive=true +;reltol=1e-10 +;abstol=1e-10 +;use_gpu=no +;gpu_id=0 +;library_file= shared_libs/libToRORd_Land_mixed_endo_mid_epi.so + +[stim_plain] +start = 0.0 +duration = 1.0 +period = 800.0 +current = -53.0 +x_limit = 500.0 +main_function=stim_if_x_less_than + +[extra_data] +; ---------------------------- +; bz1 +;Ito_Multiplier = 0.1 +;IKs_Multiplier = 0.2 +;IK1_Multiplier = 0.3 +;IKr_Multiplier = 0.7 +;INa_Multiplier = 0.4 +;ICaL_Multiplier = 0.64 +;IKCa_Multiplier = 1.0 +; ---------------------------- +; bz2 +;IKs_Multiplier = 0.2 +;IKr_Multiplier = 0.3 +;INa_Multiplier = 0.38 +;ICaL_Multiplier = 0.31 +;IKCa_Multiplier = 1.0 +; ---------------------------- +; bz3 +;Ito_Multiplier = 0.0 +;IK1_Multiplier = 0.6 +;INa_Multiplier = 0.4 +;ICaL_Multiplier = 0.64 +;aCaMK_Multiplier = 1.5 +;taurelp_Multiplier = 6.0 +;ICab_Multiplier = 1.33 +; ---------------------------- +main_function=set_extra_data_mixed_torord_Land_same_celltype \ No newline at end of file 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/extra_data_library/extra_data.c b/src/extra_data_library/extra_data.c index ef99322b..1dc64b8e 100644 --- a/src/extra_data_library/extra_data.c +++ b/src/extra_data_library/extra_data.c @@ -10,6 +10,67 @@ #include "../domains_library/mesh_info_data.h" #include "helper_functions.h" +struct uv { + real_cpu u; + real_cpu v; +}; + +struct original_volumes_hash_entry { + struct point_3d key; + struct uv value; +}; + +struct cells_hash_entry { + struct point_3d key; + struct cell_node *value; +}; + +struct point_3d find_next_point(struct original_volumes_hash_entry *original_volumes, int l, struct point_3d p, struct point_3d *n, int np) { + + double dist; + double min_dist = 100000.0; + int min_idx = -1; + + int continue_loop = false; + struct point_3d aux; + + for(int i = 0; i < l; i++) { + + aux = original_volumes[i].key; + double ax = aux.x; + double ay = aux.y; + + for(int j = 0; j < np; j++) { + if(ax == n[j].x && ay == n[j].y) { + continue_loop = true; + break; + } + } + + if(continue_loop) { + continue_loop = false; + continue; + } + + double a = ax - p.x; + double b = ay - p.y; + + dist = a*a + b*b; + + if(dist < min_dist) { + min_dist = dist; + min_idx = i; + } + } + + if(min_idx > -1) { + return original_volumes[min_idx].key; + } + else { + return POINT3D(0,0,0); + } + +} SET_EXTRA_DATA(set_extra_data_for_fibrosis_sphere) { @@ -125,211 +186,11 @@ SET_EXTRA_DATA(set_extra_data_for_benchmark) { return (void*)initial_conditions; } -SET_EXTRA_DATA (set_mixed_model_if_x_less_than) -{ - uint32_t num_active_cells = the_grid->num_active_cells; - - *extra_data_size = sizeof(uint32_t)*(num_active_cells); - - uint32_t *mapping = (uint32_t*)malloc(*extra_data_size); - - struct cell_node ** ac = the_grid->active_cells; - - real x_limit = 0.0; - GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real, x_limit, config, "x_limit"); - - int i; - bool inside; - - OMP(parallel for) - for (i = 0; i < num_active_cells; i++) - { - real center_x = ac[i]->center.x; - - inside = (center_x <= x_limit); - - if (inside) - mapping[i] = 0; - else - mapping[i] = 1; - } - - return (void*)mapping; -} - -SET_EXTRA_DATA (set_mixed_model_purkinje_and_tissue) -{ - uint32_t num_active_tissue_cells = the_grid->num_active_cells; - uint32_t num_active_purkinje_cells = the_grid->purkinje->num_active_purkinje_cells; - uint32_t num_active_cells = num_active_tissue_cells + num_active_purkinje_cells; - - *extra_data_size = sizeof(uint32_t)*(num_active_cells + 1); - - uint32_t *mapping = (uint32_t*)malloc(*extra_data_size); - - int i; - - // Purkinje section - OMP(parallel for) - for (i = 0; i < num_active_purkinje_cells; i++) { - mapping[i] = 0; - } - - // Tissue section - OMP(parallel for) - for (i = num_active_purkinje_cells; i < num_active_cells; i++) { - mapping[i] = 1; - } - - return (void*)mapping; -} - -// Initial condition - 'libToRORd_fkatp_mixed_endo_mid_epi.so' + transmurality (cable and cuboid) -SET_EXTRA_DATA(set_extra_data_mixed_torord_fkatp_epi_mid_endo) { - - uint32_t num_active_cells = the_grid->num_active_cells; - real side_length = the_grid->mesh_side_length.x; - struct cell_node ** ac = the_grid->active_cells; - - // The percentages were taken from the ToRORd paper (Transmural experiment) - real side_length_endo = side_length*0.45; - real side_length_mid = side_length_endo + side_length*0.25; - real side_length_epi = side_length_mid + side_length*0.3; - - struct extra_data_for_torord *extra_data = NULL; - - extra_data = set_common_torord_data(config, num_active_cells); - - OMP(parallel for) - for (int i = 0; i < num_active_cells; i++) { - - real center_x = ac[i]->center.x; - - // ENDO - if (center_x < side_length_endo) - extra_data->transmurality[i] = 0.0; - // MID - else if (center_x >= side_length_endo && center_x < side_length_mid) - extra_data->transmurality[i] = 1.0; - // EPI - else - extra_data->transmurality[i] = 2.0; - } - - SET_EXTRA_DATA_SIZE(sizeof(struct extra_data_for_torord)); - - return (void*)extra_data; - -} - -// Initial condition - 'libToRORd_dynCl_mixed_endo_mid_epi.so' + transmurality (cable and cuboid) -SET_EXTRA_DATA(set_extra_data_mixed_torord_dynCl_epi_mid_endo) { - - uint32_t num_active_cells = the_grid->num_active_cells; - real side_length = the_grid->mesh_side_length.x; - struct cell_node ** ac = the_grid->active_cells; - - // The percentages were taken from the ToRORd paper (Transmural experiment) - real side_length_endo = side_length*0.45; - real side_length_mid = side_length_endo + side_length*0.25; - real side_length_epi = side_length_mid + side_length*0.3; - - struct extra_data_for_torord *extra_data = NULL; - - extra_data = set_common_torord_dyncl_data(config, num_active_cells); - - OMP(parallel for) - for (int i = 0; i < num_active_cells; i++) { - - real center_x = ac[i]->center.x; - - // ENDO - if (center_x < side_length_endo) - extra_data->transmurality[i] = 0.0; - // MID - else if (center_x >= side_length_endo && center_x < side_length_mid) - extra_data->transmurality[i] = 1.0; - // EPI - else - extra_data->transmurality[i] = 2.0; - } - - SET_EXTRA_DATA_SIZE(sizeof(struct extra_data_for_torord)); - - return (void*)extra_data; - -} - -struct uv { - real_cpu u; - real_cpu v; -}; - -struct original_volumes_hash_entry { - struct point_3d key; - struct uv value; -}; - -struct cells_hash_entry { - struct point_3d key; - struct cell_node *value; -}; - - -struct point_3d find_next_point(struct original_volumes_hash_entry *original_volumes, int l, struct point_3d p, struct point_3d *n, int np) { - - double dist; - double min_dist = 100000.0; - int min_idx = -1; - - int continue_loop = false; - struct point_3d aux; - - for(int i = 0; i < l; i++) { - - aux = original_volumes[i].key; - double ax = aux.x; - double ay = aux.y; - - for(int j = 0; j < np; j++) { - if(ax == n[j].x && ay == n[j].y) { - continue_loop = true; - break; - } - } - - if(continue_loop) { - continue_loop = false; - continue; - } - - double a = ax - p.x; - double b = ay - p.y; - - dist = a*a + b*b; - - if(dist < min_dist) { - min_dist = dist; - min_idx = i; - } - } - - if(min_idx > -1) { - return original_volumes[min_idx].key; - } - else { - return POINT3D(0,0,0); - } - -} - SET_EXTRA_DATA(set_extra_data_for_spiral_fhn) { - size_t num_sv_entries = 2; uint32_t num_active_cells = the_grid->num_active_cells; - *extra_data_size = num_active_cells * num_sv_entries; real *sv_cpu; @@ -509,6 +370,202 @@ SET_EXTRA_DATA(set_extra_data_for_spiral_fhn) { return sv_cpu; } +SET_EXTRA_DATA (set_mixed_model_if_x_less_than) +{ + uint32_t num_active_cells = the_grid->num_active_cells; + + *extra_data_size = sizeof(uint32_t)*(num_active_cells); + + uint32_t *mapping = (uint32_t*)malloc(*extra_data_size); + + struct cell_node ** ac = the_grid->active_cells; + + real x_limit = 0.0; + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real, x_limit, config, "x_limit"); + + int i; + bool inside; + + OMP(parallel for) + for (i = 0; i < num_active_cells; i++) + { + real center_x = ac[i]->center.x; + + inside = (center_x <= x_limit); + + if (inside) + mapping[i] = 0; + else + mapping[i] = 1; + } + + return (void*)mapping; +} + +SET_EXTRA_DATA (set_mixed_model_purkinje_and_tissue) +{ + uint32_t num_active_tissue_cells = the_grid->num_active_cells; + uint32_t num_active_purkinje_cells = the_grid->purkinje->num_active_purkinje_cells; + uint32_t num_active_cells = num_active_tissue_cells + num_active_purkinje_cells; + + *extra_data_size = sizeof(uint32_t)*(num_active_cells + 1); + + uint32_t *mapping = (uint32_t*)malloc(*extra_data_size); + + int i; + + // Purkinje section + OMP(parallel for) + for (i = 0; i < num_active_purkinje_cells; i++) { + mapping[i] = 0; + } + + // Tissue section + OMP(parallel for) + for (i = num_active_purkinje_cells; i < num_active_cells; i++) { + mapping[i] = 1; + } + + return (void*)mapping; +} + +// Initial condition - 'libToRORd_fkatp_mixed_endo_mid_epi.so' + transmurality + current modifiers (plain and cuboid) +SET_EXTRA_DATA(set_extra_data_mixed_torord_fkatp_epi_mid_endo) { + + uint32_t num_active_cells = the_grid->num_active_cells; + real side_length = the_grid->mesh_side_length.x; + struct cell_node ** ac = the_grid->active_cells; + + // The percentages were taken from the ToRORd paper (Transmural experiment) + real side_length_endo = side_length*0.45; + real side_length_mid = side_length_endo + side_length*0.25; + real side_length_epi = side_length_mid + side_length*0.3; + + struct extra_data_for_torord *extra_data = NULL; + extra_data = set_common_torord_data(config, num_active_cells); + + OMP(parallel for) + for (int i = 0; i < num_active_cells; i++) { + + real center_x = ac[i]->center.x; + + // ENDO + if (center_x < side_length_endo) + extra_data->transmurality[i] = 0.0; + // MID + else if (center_x >= side_length_endo && center_x < side_length_mid) + extra_data->transmurality[i] = 1.0; + // EPI + else + extra_data->transmurality[i] = 2.0; + } + + SET_EXTRA_DATA_SIZE(sizeof(struct extra_data_for_torord)); + + return (void*)extra_data; + +} + +// Initial condition - 'libToRORd_dynCl_mixed_endo_mid_epi.so' + transmurality + current modifiers (plain and cuboid) +SET_EXTRA_DATA(set_extra_data_mixed_torord_dynCl_epi_mid_endo) { + + uint32_t num_active_cells = the_grid->num_active_cells; + real side_length = the_grid->mesh_side_length.x; + struct cell_node ** ac = the_grid->active_cells; + + // The percentages were taken from the ToRORd paper (Transmural experiment) + real side_length_endo = side_length*0.45; + real side_length_mid = side_length_endo + side_length*0.25; + real side_length_epi = side_length_mid + side_length*0.3; + + struct extra_data_for_torord *extra_data = NULL; + extra_data = set_common_torord_dyncl_data(config, num_active_cells); + + OMP(parallel for) + for (int i = 0; i < num_active_cells; i++) { + + real center_x = ac[i]->center.x; + + // ENDO + if (center_x < side_length_endo) + extra_data->transmurality[i] = 0.0; + // MID + else if (center_x >= side_length_endo && center_x < side_length_mid) + extra_data->transmurality[i] = 1.0; + // EPI + else + extra_data->transmurality[i] = 2.0; + } + + SET_EXTRA_DATA_SIZE(sizeof(struct extra_data_for_torord)); + + return (void*)extra_data; +} + +// Initial condition - 'libToRORd_Land_mixed_endo_mid_epi.so' + transmurality + current modifiers (plain and cuboid) +SET_EXTRA_DATA(set_extra_data_mixed_torord_Land_epi_mid_endo) { + + uint32_t num_active_cells = the_grid->num_active_cells; + real side_length = the_grid->mesh_side_length.x; + struct cell_node ** ac = the_grid->active_cells; + + // The percentages were taken from the ToRORd paper (Transmural experiment) + real side_length_endo = side_length*0.45; + real side_length_mid = side_length_endo + side_length*0.25; + real side_length_epi = side_length_mid + side_length*0.3; + + struct extra_data_for_torord_land *extra_data = NULL; + extra_data = set_common_torord_Land_data(config, num_active_cells); + + OMP(parallel for) + for (int i = 0; i < num_active_cells; i++) { + + real center_x = ac[i]->center.x; + + // ENDO + if (center_x < side_length_endo) + extra_data->transmurality[i] = 0.0; + // MID + else if (center_x >= side_length_endo && center_x < side_length_mid) + extra_data->transmurality[i] = 1.0; + // EPI + else + extra_data->transmurality[i] = 2.0; + } + + SET_EXTRA_DATA_SIZE(sizeof(struct extra_data_for_torord_land)); + + return (void*)extra_data; + +} + +// Initial condition - 'libToRORd_Land_mixed_endo_mid_epi.so' + current modifiers (cell, plain and cuboid) +SET_EXTRA_DATA(set_extra_data_mixed_torord_Land_same_celltype) { + + uint32_t num_active_cells = the_grid->num_active_cells; + struct cell_node ** ac = the_grid->active_cells; + + struct extra_data_for_torord_land *extra_data = NULL; + extra_data = set_common_torord_Land_data(config, num_active_cells); + + // All cells will be the same type + OMP(parallel for) + for (int i = 0; i < num_active_cells; i++) { + // ENDO + extra_data->transmurality[i] = 0.0; + // MID + //extra_data->transmurality[i] = 1.0; + // EPI + //extra_data->transmurality[i] = 2.0; + } + + SET_EXTRA_DATA_SIZE(sizeof(struct extra_data_for_torord_land)); + + return (void*)extra_data; + +} + +// Current modifiers - 'libtrovato_2020.so' SET_EXTRA_DATA(set_extra_data_trovato) { if (!the_grid->purkinje) { diff --git a/src/extra_data_library/helper_functions.c b/src/extra_data_library/helper_functions.c index 6110a4a9..8d709992 100644 --- a/src/extra_data_library/helper_functions.c +++ b/src/extra_data_library/helper_functions.c @@ -482,6 +482,234 @@ struct extra_data_for_torord * set_common_torord_dyncl_data (struct config *conf return extra_data; } +struct extra_data_for_torord_land * set_common_torord_Land_data (struct config *config, uint32_t num_cells) { + const uint32_t neq = 49; + struct extra_data_for_torord_land *extra_data = MALLOC_ONE_TYPE(struct extra_data_for_torord_land); + + real INa_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, INa_Multiplier, config, "INa_Multiplier"); + real INaL_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, INaL_Multiplier, config, "INaL_Multiplier"); + real INaCa_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, INaCa_Multiplier, config, "INaCa_Multiplier"); + real INaK_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, INaK_Multiplier, config, "INaK_Multiplier"); + real INab_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, INab_Multiplier, config, "INab_Multiplier"); + real Ito_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, Ito_Multiplier, config, "Ito_Multiplier"); + real IKr_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, IKr_Multiplier, config, "IKr_Multiplier"); + real IKs_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, IKs_Multiplier, config, "IKs_Multiplier"); + real IK1_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, IK1_Multiplier, config, "IK1_Multiplier"); + real IKb_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, IKb_Multiplier, config, "IKb_Multiplier"); + real IKCa_Multiplier = 0.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, IKCa_Multiplier, config, "IKCa_Multiplier"); + real ICaL_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, ICaL_Multiplier, config, "ICaL_Multiplier"); + real ICab_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, ICab_Multiplier, config, "ICab_Multiplier"); + real IpCa_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, IpCa_Multiplier, config, "IpCa_Multiplier"); + real ICaCl_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, ICaCl_Multiplier, config, "ICaCl_Multiplier"); + real IClb_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, IClb_Multiplier, config, "IClb_Multiplier"); + real Jrel_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, Jrel_Multiplier, config, "Jrel_Multiplier"); + real Jup_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, Jup_Multiplier, config, "Jup_Multiplier"); + real aCaMK_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, aCaMK_Multiplier, config, "aCaMK_Multiplier"); + real taurelp_Multiplier = 1.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, taurelp_Multiplier, config, "taurelp_Multiplier"); + + extra_data->INa_Multiplier = INa_Multiplier; + extra_data->INaL_Multiplier = INaL_Multiplier; + extra_data->INaCa_Multiplier = INaCa_Multiplier; + extra_data->INaK_Multiplier = INaK_Multiplier; + extra_data->INab_Multiplier = INab_Multiplier; + extra_data->Ito_Multiplier = Ito_Multiplier; + extra_data->IKr_Multiplier = IKr_Multiplier; + extra_data->IKs_Multiplier = IKs_Multiplier; + extra_data->IK1_Multiplier = IK1_Multiplier; + extra_data->IKb_Multiplier = IKb_Multiplier; + extra_data->IKCa_Multiplier = IKCa_Multiplier; + extra_data->ICaL_Multiplier = ICaL_Multiplier; + extra_data->ICab_Multiplier = ICab_Multiplier; + extra_data->IpCa_Multiplier = IpCa_Multiplier; + extra_data->ICaCl_Multiplier = ICaCl_Multiplier; + extra_data->IClb_Multiplier = IClb_Multiplier; + extra_data->Jrel_Multiplier = Jrel_Multiplier; + extra_data->Jup_Multiplier = Jup_Multiplier; + extra_data->aCaMK_Multiplier = aCaMK_Multiplier; + extra_data->taurelp_Multiplier = taurelp_Multiplier; + + extra_data->initial_ss_endo = MALLOC_ARRAY_OF_TYPE(real, neq); + extra_data->initial_ss_epi = MALLOC_ARRAY_OF_TYPE(real, neq); + extra_data->initial_ss_mid = MALLOC_ARRAY_OF_TYPE(real, neq); + + // Default initial conditions for ENDO cell (from original Matlab script) + extra_data->initial_ss_endo[0] = -8.863699e+01; + extra_data->initial_ss_endo[1] = 1.189734e+01; + extra_data->initial_ss_endo[2] = 1.189766e+01; + extra_data->initial_ss_endo[3] = 1.412345e+02; + extra_data->initial_ss_endo[4] = 1.412344e+02; + extra_data->initial_ss_endo[5] = 7.267473e-05; + extra_data->initial_ss_endo[6] = 6.337870e-05; + extra_data->initial_ss_endo[7] = 1.532653e+00; + extra_data->initial_ss_endo[8] = 1.533946e+00; + extra_data->initial_ss_endo[9] = 8.280078e-04; + extra_data->initial_ss_endo[10] = 6.665272e-01; + extra_data->initial_ss_endo[11] = 8.260208e-01; + extra_data->initial_ss_endo[12] = 8.260560e-01; + extra_data->initial_ss_endo[13] = 8.258509e-01; + extra_data->initial_ss_endo[14] = 1.668686e-04; + extra_data->initial_ss_endo[15] = 5.228306e-01; + extra_data->initial_ss_endo[16] = 2.859696e-01; + extra_data->initial_ss_endo[17] = 9.591370e-04; + extra_data->initial_ss_endo[18] = 9.996012e-01; + extra_data->initial_ss_endo[19] = 5.934016e-01; + extra_data->initial_ss_endo[20] = 4.886961e-04; + extra_data->initial_ss_endo[21] = 9.996011e-01; + extra_data->initial_ss_endo[22] = 6.546687e-01; + extra_data->initial_ss_endo[23] = 9.500075e-32; + extra_data->initial_ss_endo[24] = 1.000000e+00; + extra_data->initial_ss_endo[25] = 9.392580e-01; + extra_data->initial_ss_endo[26] = 1.000000e+00; + extra_data->initial_ss_endo[27] = 9.998984e-01; + extra_data->initial_ss_endo[28] = 9.999783e-01; + extra_data->initial_ss_endo[29] = 4.448162e-04; + extra_data->initial_ss_endo[30] = 7.550725e-04; + extra_data->initial_ss_endo[31] = 1.000000e+00; + extra_data->initial_ss_endo[32] = 1.000000e+00; + extra_data->initial_ss_endo[33] = 2.424047e-01; + extra_data->initial_ss_endo[34] = 1.795377e-04; + extra_data->initial_ss_endo[35] = -6.883086e-25; + extra_data->initial_ss_endo[36] = 1.117498e-02; + extra_data->initial_ss_endo[37] = 9.980366e-01; + extra_data->initial_ss_endo[38] = 8.588018e-04; + extra_data->initial_ss_endo[39] = 7.097447e-04; + extra_data->initial_ss_endo[40] = 3.812617e-04; + extra_data->initial_ss_endo[41] = 1.357116e-05; + extra_data->initial_ss_endo[42] = 2.302525e-23; + extra_data->initial_ss_endo[43] = 1.561941e-04; + extra_data->initial_ss_endo[44] = 2.351289e-04; + extra_data->initial_ss_endo[45] = 8.077631e-03; + extra_data->initial_ss_endo[46] = 9.993734e-01; + extra_data->initial_ss_endo[47] = 0.000000e+00; + extra_data->initial_ss_endo[48] = 0.000000e+00; + + // Default initial conditions for EPI cell (from original Matlab script) + extra_data->initial_ss_epi[0] = -8.904628e+01; + extra_data->initial_ss_epi[1] = 1.272190e+01; + extra_data->initial_ss_epi[2] = 1.272220e+01; + extra_data->initial_ss_epi[3] = 1.422490e+02; + extra_data->initial_ss_epi[4] = 1.422489e+02; + extra_data->initial_ss_epi[5] = 6.541058e-05; + extra_data->initial_ss_epi[6] = 5.684431e-05; + extra_data->initial_ss_epi[7] = 1.809117e+00; + extra_data->initial_ss_epi[8] = 1.809702e+00; + extra_data->initial_ss_epi[9] = 7.581821e-04; + extra_data->initial_ss_epi[10] = 6.798398e-01; + extra_data->initial_ss_epi[11] = 8.341502e-01; + extra_data->initial_ss_epi[12] = 8.341883e-01; + extra_data->initial_ss_epi[13] = 8.340817e-01; + extra_data->initial_ss_epi[14] = 1.543877e-04; + extra_data->initial_ss_epi[15] = 5.382951e-01; + extra_data->initial_ss_epi[16] = 3.027694e-01; + extra_data->initial_ss_epi[17] = 9.330351e-04; + extra_data->initial_ss_epi[18] = 9.996287e-01; + extra_data->initial_ss_epi[19] = 9.996262e-01; + extra_data->initial_ss_epi[20] = 4.753907e-04; + extra_data->initial_ss_epi[21] = 9.996287e-01; + extra_data->initial_ss_epi[22] = 9.996285e-01; + extra_data->initial_ss_epi[23] = 1.742134e-37; + extra_data->initial_ss_epi[24] = 1.000000e+00; + extra_data->initial_ss_epi[25] = 9.479522e-01; + extra_data->initial_ss_epi[26] = 1.000000e+00; + extra_data->initial_ss_epi[27] = 9.999327e-01; + extra_data->initial_ss_epi[28] = 9.999829e-01; + extra_data->initial_ss_epi[29] = 2.915447e-04; + extra_data->initial_ss_epi[30] = 5.026045e-04; + extra_data->initial_ss_epi[31] = 1.000000e+00; + extra_data->initial_ss_epi[32] = 1.000000e+00; + extra_data->initial_ss_epi[33] = 2.288155e-01; + extra_data->initial_ss_epi[34] = 1.714978e-04; + extra_data->initial_ss_epi[35] = -1.131190e-26; + extra_data->initial_ss_epi[36] = 1.295052e-02; + extra_data->initial_ss_epi[37] = 9.981944e-01; + extra_data->initial_ss_epi[38] = 8.342321e-04; + extra_data->initial_ss_epi[39] = 6.838658e-04; + extra_data->initial_ss_epi[40] = 2.778785e-04; + extra_data->initial_ss_epi[41] = 9.667759e-06; + extra_data->initial_ss_epi[42] = 8.169304e-24; + extra_data->initial_ss_epi[43] = 1.259996e-04; + extra_data->initial_ss_epi[44] = 1.899522e-04; + extra_data->initial_ss_epi[45] = 6.551494e-03; + extra_data->initial_ss_epi[46] = 9.994940e-01; + extra_data->initial_ss_epi[47] = 0.000000e+00; + extra_data->initial_ss_epi[48] = 0.000000e+00; + + // Default initial conditions for MID cell (from original Matlab script) + extra_data->initial_ss_mid[0] = -8.953800e+01; + extra_data->initial_ss_mid[1] = 1.492920e+01; + extra_data->initial_ss_mid[2] = 1.492967e+01; + extra_data->initial_ss_mid[3] = 1.448447e+02; + extra_data->initial_ss_mid[4] = 1.448447e+02; + extra_data->initial_ss_mid[5] = 7.502288e-05; + extra_data->initial_ss_mid[6] = 6.107636e-05; + extra_data->initial_ss_mid[7] = 1.790435e+00; + extra_data->initial_ss_mid[8] = 1.794842e+00; + extra_data->initial_ss_mid[9] = 6.819365e-04; + extra_data->initial_ss_mid[10] = 6.953807e-01; + extra_data->initial_ss_mid[11] = 8.434888e-01; + extra_data->initial_ss_mid[12] = 8.435208e-01; + extra_data->initial_ss_mid[13] = 8.432262e-01; + extra_data->initial_ss_mid[14] = 1.406211e-04; + extra_data->initial_ss_mid[15] = 5.453149e-01; + extra_data->initial_ss_mid[16] = 2.924967e-01; + extra_data->initial_ss_mid[17] = 9.026127e-04; + extra_data->initial_ss_mid[18] = 9.996593e-01; + extra_data->initial_ss_mid[19] = 5.631197e-01; + extra_data->initial_ss_mid[20] = 4.598833e-04; + extra_data->initial_ss_mid[21] = 9.996593e-01; + extra_data->initial_ss_mid[22] = 6.236964e-01; + extra_data->initial_ss_mid[23] = -1.314189e-33; + extra_data->initial_ss_mid[24] = 1.000000e+00; + extra_data->initial_ss_mid[25] = 9.204086e-01; + extra_data->initial_ss_mid[26] = 1.000000e+00; + extra_data->initial_ss_mid[27] = 9.997620e-01; + extra_data->initial_ss_mid[28] = 9.999625e-01; + extra_data->initial_ss_mid[29] = 3.853595e-04; + extra_data->initial_ss_mid[30] = 8.535292e-04; + extra_data->initial_ss_mid[31] = 1.000000e+00; + extra_data->initial_ss_mid[32] = 1.000000e+00; + extra_data->initial_ss_mid[33] = 2.664151e-01; + extra_data->initial_ss_mid[34] = 1.623107e-04; + extra_data->initial_ss_mid[35] = 1.209762e-24; + extra_data->initial_ss_mid[36] = 1.782437e-02; + extra_data->initial_ss_mid[37] = 9.979720e-01; + extra_data->initial_ss_mid[38] = 8.053991e-04; + extra_data->initial_ss_mid[39] = 6.781800e-04; + extra_data->initial_ss_mid[40] = 5.265363e-04; + extra_data->initial_ss_mid[41] = 1.789565e-05; + extra_data->initial_ss_mid[42] = 7.059162e-23; + extra_data->initial_ss_mid[43] = 1.670654e-04; + extra_data->initial_ss_mid[44] = 2.506794e-04; + extra_data->initial_ss_mid[45] = 8.602625e-03; + extra_data->initial_ss_mid[46] = 9.993314e-01; + extra_data->initial_ss_mid[47] = 0.000000e+00; + extra_data->initial_ss_mid[48] = 0.000000e+00; + + extra_data->transmurality = MALLOC_ARRAY_OF_TYPE(real, num_cells); + + return extra_data; +} + struct extra_data_for_trovato * set_common_trovato_data (struct config *config, uint32_t num_cells) { struct extra_data_for_trovato *extra_data = MALLOC_ONE_TYPE(struct extra_data_for_trovato); diff --git a/src/extra_data_library/helper_functions.h b/src/extra_data_library/helper_functions.h index 8685624a..b15fdadb 100644 --- a/src/extra_data_library/helper_functions.h +++ b/src/extra_data_library/helper_functions.h @@ -46,6 +46,33 @@ struct extra_data_for_torord { real *transmurality; }; +struct extra_data_for_torord_land { + real INa_Multiplier; + real INaL_Multiplier; + real INaCa_Multiplier; + real INaK_Multiplier; + real INab_Multiplier; + real Ito_Multiplier; + real IKr_Multiplier; + real IKs_Multiplier; + real IK1_Multiplier; + real IKb_Multiplier; + real IKCa_Multiplier; + real ICaL_Multiplier; + real ICab_Multiplier; + real IpCa_Multiplier; + real ICaCl_Multiplier; + real IClb_Multiplier; + real Jrel_Multiplier; + real Jup_Multiplier; + real aCaMK_Multiplier; + real taurelp_Multiplier; + real *initial_ss_endo; + real *initial_ss_epi; + real *initial_ss_mid; + real *transmurality; +}; + struct extra_data_for_trovato { real GNa_Multiplier; real GNaL_Multiplier; @@ -82,6 +109,7 @@ struct extra_data_for_trovato { struct extra_data_for_fibrosis * set_common_schemia_data(struct config *config, uint32_t num_cells); struct extra_data_for_torord * set_common_torord_data (struct config *config, uint32_t num_cells); struct extra_data_for_torord * set_common_torord_dyncl_data (struct config *config, uint32_t num_cells); +struct extra_data_for_torord_land * set_common_torord_Land_data (struct config *config, uint32_t num_cells); struct extra_data_for_trovato * set_common_trovato_data (struct config *config, uint32_t num_cells); #endif // MONOALG3D_C_EXTRA_DATA_HELPER_FUNCTIONS_H diff --git a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.c b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.c new file mode 100644 index 00000000..bf32adbd --- /dev/null +++ b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.c @@ -0,0 +1,682 @@ +#include "ToRORd_Land_mixed_endo_mid_epi.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 ToRORd_Land_2020 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 timestep model to solve the ODEs\n"); + } else { + log_info("Using Fixed timestep to solve the ODEs\n"); + } + + real *initial_endo = NULL; + real *initial_epi = NULL; + real *initial_mid = NULL; + real *transmurality = NULL; + if(solver->ode_extra_data) { + struct extra_data_for_torord_land *extra_data = (struct extra_data_for_torord_land*)solver->ode_extra_data; + initial_endo = extra_data->initial_ss_endo; + initial_epi = extra_data->initial_ss_epi; + initial_mid = extra_data->initial_ss_mid; + transmurality = extra_data->transmurality; + + OMP(parallel for) + for(uint32_t i = 0; i < num_cells; i++){ + + real *sv = &solver->sv[i * NEQ]; + + for (int j = 0; j < NEQ; j++) { + if (transmurality[i] == ENDO) + sv[j] = initial_endo[j]; + else if (transmurality[i] == EPI) + sv[j] = initial_epi[j]; + else + sv[j] = initial_mid[j]; + } + } + } + else { + log_info("[INFO] You should supply a mask function to tag the cells when using this mixed model!\n"); + log_info("[INFO] Considering all cells ENDO!\n"); + + OMP(parallel for) + for(uint32_t i = 0; i < num_cells; i++){ + + real *sv = &solver->sv[i * NEQ]; + + // Default initial conditions for ENDO cell (from original Matlab script) + sv[0] = -8.863699e+01; + sv[1] = 1.189734e+01; + sv[2] = 1.189766e+01; + sv[3] = 1.412345e+02; + sv[4] = 1.412344e+02; + sv[5] = 7.267473e-05; + sv[6] = 6.337870e-05; + sv[7] = 1.532653e+00; + sv[8] = 1.533946e+00; + sv[9] = 8.280078e-04; + sv[10] = 6.665272e-01; + sv[11] = 8.260208e-01; + sv[12] = 8.260560e-01; + sv[13] = 8.258509e-01; + sv[14] = 1.668686e-04; + sv[15] = 5.228306e-01; + sv[16] = 2.859696e-01; + sv[17] = 9.591370e-04; + sv[18] = 9.996012e-01; + sv[19] = 5.934016e-01; + sv[20] = 4.886961e-04; + sv[21] = 9.996011e-01; + sv[22] = 6.546687e-01; + sv[23] = 9.500075e-32; + sv[24] = 1.000000e+00; + sv[25] = 9.392580e-01; + sv[26] = 1.000000e+00; + sv[27] = 9.998984e-01; + sv[28] = 9.999783e-01; + sv[29] = 4.448162e-04; + sv[30] = 7.550725e-04; + sv[31] = 1.000000e+00; + sv[32] = 1.000000e+00; + sv[33] = 2.424047e-01; + sv[34] = 1.795377e-04; + sv[35] = -6.883086e-25; + sv[36] = 1.117498e-02; + sv[37] = 9.980366e-01; + sv[38] = 8.588018e-04; + sv[39] = 7.097447e-04; + sv[40] = 3.812617e-04; + sv[41] = 1.357116e-05; + sv[42] = 2.302525e-23; + sv[43] = 1.561941e-04; + sv[44] = 2.351289e-04; + sv[45] = 8.077631e-03; + sv[46] = 9.993734e-01; + sv[47] = 0.000000e+00; + sv[48] = 0.000000e+00; + + // Default initial conditions for MID cell (from original Matlab script) + //sv[0] = -8.953800e+01; + //sv[1] = 1.492920e+01; + //sv[2] = 1.492967e+01; + //sv[3] = 1.448447e+02; + //sv[4] = 1.448447e+02; + //sv[5] = 7.502288e-05; + //sv[6] = 6.107636e-05; + //sv[7] = 1.790435e+00; + //sv[8] = 1.794842e+00; + //sv[9] = 6.819365e-04; + //sv[10] = 6.953807e-01; + //sv[11] = 8.434888e-01; + //sv[12] = 8.435208e-01; + //sv[13] = 8.432262e-01; + //sv[14] = 1.406211e-04; + //sv[15] = 5.453149e-01; + //sv[16] = 2.924967e-01; + //sv[17] = 9.026127e-04; + //sv[18] = 9.996593e-01; + //sv[19] = 5.631197e-01; + //sv[20] = 4.598833e-04; + //sv[21] = 9.996593e-01; + //sv[22] = 6.236964e-01; + //sv[23] = -1.314189e-33; + //sv[24] = 1.000000e+00; + //sv[25] = 9.204086e-01; + //sv[26] = 1.000000e+00; + //sv[27] = 9.997620e-01; + //sv[28] = 9.999625e-01; + //sv[29] = 3.853595e-04; + //sv[30] = 8.535292e-04; + //sv[31] = 1.000000e+00; + //sv[32] = 1.000000e+00; + //sv[33] = 2.664151e-01; + //sv[34] = 1.623107e-04; + //sv[35] = 1.209762e-24; + //sv[36] = 1.782437e-02; + //sv[37] = 9.979720e-01; + //sv[38] = 8.053991e-04; + //sv[39] = 6.781800e-04; + //sv[40] = 5.265363e-04; + //sv[41] = 1.789565e-05; + //sv[42] = 7.059162e-23; + //sv[43] = 1.670654e-04; + //sv[44] = 2.506794e-04; + //sv[45] = 8.602625e-03; + //sv[46] = 9.993314e-01; + //sv[47] = 0.000000e+00; + //sv[48] = 0.000000e+00; + + // Default initial conditions for EPI cell (from original Matlab script) + //sv[0] = -8.904628e+01; + //sv[1] = 1.272190e+01; + //sv[2] = 1.272220e+01; + //sv[3] = 1.422490e+02; + //sv[4] = 1.422489e+02; + //sv[5] = 6.541058e-05; + //sv[6] = 5.684431e-05; + //sv[7] = 1.809117e+00; + //sv[8] = 1.809702e+00; + //sv[9] = 7.581821e-04; + //sv[10] = 6.798398e-01; + //sv[11] = 8.341502e-01; + //sv[12] = 8.341883e-01; + //sv[13] = 8.340817e-01; + //sv[14] = 1.543877e-04; + //sv[15] = 5.382951e-01; + //sv[16] = 3.027694e-01; + //sv[17] = 9.330351e-04; + //sv[18] = 9.996287e-01; + //sv[19] = 9.996262e-01; + //sv[20] = 4.753907e-04; + //sv[21] = 9.996287e-01; + //sv[22] = 9.996285e-01; + //sv[23] = 1.742134e-37; + //sv[24] = 1.000000e+00; + //sv[25] = 9.479522e-01; + //sv[26] = 1.000000e+00; + //sv[27] = 9.999327e-01; + //sv[28] = 9.999829e-01; + //sv[29] = 2.915447e-04; + //sv[30] = 5.026045e-04; + //sv[31] = 1.000000e+00; + //sv[32] = 1.000000e+00; + //sv[33] = 2.288155e-01; + //sv[34] = 1.714978e-04; + //sv[35] = -1.131190e-26; + //sv[36] = 1.295052e-02; + //sv[37] = 9.981944e-01; + //sv[38] = 8.342321e-04; + //sv[39] = 6.838658e-04; + //sv[40] = 2.778785e-04; + //sv[41] = 9.667759e-06; + //sv[42] = 8.169304e-24; + //sv[43] = 1.259996e-04; + //sv[44] = 1.899522e-04; + //sv[45] = 6.551494e-03; + //sv[46] = 9.994940e-01; + //sv[47] = 0.000000e+00; + //sv[48] = 0.000000e+00; + } + } +} + +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; + bool adpt = ode_solver->adaptive; + + // Get the extra parameters + int num_extra_parameters = 20; + real extra_par[num_extra_parameters]; + real *transmurality = NULL; + if (ode_solver->ode_extra_data) { + struct extra_data_for_torord_land *extra_data = (struct extra_data_for_torord_land*)ode_solver->ode_extra_data; + extra_par[0] = extra_data->INa_Multiplier; + extra_par[1] = extra_data->INaL_Multiplier; + extra_par[2] = extra_data->INaCa_Multiplier; + extra_par[3] = extra_data->INaK_Multiplier; + extra_par[4] = extra_data->INab_Multiplier; + extra_par[5] = extra_data->Ito_Multiplier; + extra_par[6] = extra_data->IKr_Multiplier; + extra_par[7] = extra_data->IKs_Multiplier; + extra_par[8] = extra_data->IK1_Multiplier; + extra_par[9] = extra_data->IKb_Multiplier; + extra_par[10] = extra_data->IKCa_Multiplier; + extra_par[11] = extra_data->ICaL_Multiplier; + extra_par[12] = extra_data->ICab_Multiplier; + extra_par[13] = extra_data->IpCa_Multiplier; + extra_par[14] = extra_data->ICaCl_Multiplier; + extra_par[15] = extra_data->IClb_Multiplier; + extra_par[16] = extra_data->Jrel_Multiplier; + extra_par[17] = extra_data->Jup_Multiplier; + extra_par[18] = extra_data->aCaMK_Multiplier; + extra_par[19] = extra_data->taurelp_Multiplier; + transmurality = extra_data->transmurality; + } + else { + // Default: initialize all current modifiers + for (uint32_t i = 0; i < num_extra_parameters; i++) { + if (i == 9) + extra_par[i] = 0.0; + else + extra_par[i] = 1.0; + } + } + + // Solve the ODEs + OMP(parallel for private(sv_id)) + for (u_int32_t i = 0; i < num_cells_to_solve; i++) { + + if(cells_to_solve) + sv_id = cells_to_solve[i]; + else + sv_id = i; + + if(adpt) { + if (ode_solver->ode_extra_data) { + solve_forward_euler_cpu_adpt(sv + (sv_id * NEQ), stim_currents[i], transmurality[i], current_t + dt, sv_id, ode_solver, extra_par); + } + else { + solve_forward_euler_cpu_adpt(sv + (sv_id * NEQ), stim_currents[i], 0.0, current_t + dt, sv_id, ode_solver, extra_par); + } + } + else { + for (int j = 0; j < num_steps; ++j) { + if (ode_solver->ode_extra_data) { + solve_model_ode_cpu(dt, sv + (sv_id * NEQ), stim_currents[i], transmurality[i], extra_par); + } + else { + solve_model_ode_cpu(dt, sv + (sv_id * NEQ), stim_currents[i], 0.0, extra_par); + } + } + } + } +} + +void solve_model_ode_cpu(real dt, real *sv, real stim_current, real transmurality, real const *extra_params) { + + const real TOLERANCE = 1e-8; + real rY[NEQ], rDY[NEQ]; + + for(int i = 0; i < NEQ; i++) + rY[i] = sv[i]; + + // Compute 'a', 'b' coefficients alongside 'rhs' + real a[NEQ], b[NEQ]; + RHS_RL_cpu(a, b, sv, rDY, stim_current, dt, transmurality, extra_params); + + // Solve variables based on its type: + // Non-linear = Euler + // Hodkin-Huxley = Rush-Larsen || Euler (if 'a' coefficient is too small) + SOLVE_EQUATION_EULER_CPU(0); // v + SOLVE_EQUATION_EULER_CPU(1); // nai + SOLVE_EQUATION_EULER_CPU(2); // nass + SOLVE_EQUATION_EULER_CPU(3); // ki + SOLVE_EQUATION_EULER_CPU(4); // kss + SOLVE_EQUATION_EULER_CPU(5); // cai + SOLVE_EQUATION_EULER_CPU(6); // cass + SOLVE_EQUATION_EULER_CPU(7); // cansr + SOLVE_EQUATION_EULER_CPU(8); // cajsr + SOLVE_EQUATION_RUSH_LARSEN_CPU(9); // m + SOLVE_EQUATION_RUSH_LARSEN_CPU(10); // hp + SOLVE_EQUATION_RUSH_LARSEN_CPU(11); // h + SOLVE_EQUATION_RUSH_LARSEN_CPU(12); // j + SOLVE_EQUATION_RUSH_LARSEN_CPU(13); // jp + SOLVE_EQUATION_RUSH_LARSEN_CPU(14); // mL + SOLVE_EQUATION_RUSH_LARSEN_CPU(15); // hL + SOLVE_EQUATION_RUSH_LARSEN_CPU(16); // hLp + SOLVE_EQUATION_RUSH_LARSEN_CPU(17); // a + SOLVE_EQUATION_RUSH_LARSEN_CPU(18); // iF + SOLVE_EQUATION_RUSH_LARSEN_CPU(19); // iS + SOLVE_EQUATION_RUSH_LARSEN_CPU(20); // ap + SOLVE_EQUATION_RUSH_LARSEN_CPU(21); // iFp + SOLVE_EQUATION_RUSH_LARSEN_CPU(22); // iSp + SOLVE_EQUATION_RUSH_LARSEN_CPU(23); // d + SOLVE_EQUATION_RUSH_LARSEN_CPU(24); // ff + SOLVE_EQUATION_RUSH_LARSEN_CPU(25); // fs + SOLVE_EQUATION_RUSH_LARSEN_CPU(26); // fcaf + SOLVE_EQUATION_RUSH_LARSEN_CPU(27); // fcas + SOLVE_EQUATION_RUSH_LARSEN_CPU(28); // jca + SOLVE_EQUATION_EULER_CPU(29); // nca + SOLVE_EQUATION_EULER_CPU(30); // nca_i + SOLVE_EQUATION_RUSH_LARSEN_CPU(31); // ffp + SOLVE_EQUATION_RUSH_LARSEN_CPU(32); // fcafp + SOLVE_EQUATION_RUSH_LARSEN_CPU(33); // xs1 + SOLVE_EQUATION_RUSH_LARSEN_CPU(34); // xs2 + SOLVE_EQUATION_RUSH_LARSEN_CPU(35); // Jrel_np + SOLVE_EQUATION_EULER_CPU(36); // CaMKt + SOLVE_EQUATION_EULER_CPU(37); // ikr_c0 + SOLVE_EQUATION_EULER_CPU(38); // ikr_c1 + SOLVE_EQUATION_EULER_CPU(39); // ikr_c2 + SOLVE_EQUATION_EULER_CPU(40); // ikr_o + SOLVE_EQUATION_EULER_CPU(41); // ikr_i + SOLVE_EQUATION_RUSH_LARSEN_CPU(42); // Jrel_p + // --------------------------------------------------- + // Land-Niederer + SOLVE_EQUATION_EULER_CPU(43); // XS + SOLVE_EQUATION_EULER_CPU(44); // XW + SOLVE_EQUATION_EULER_CPU(45); // Ca_TRPN + SOLVE_EQUATION_EULER_CPU(46); // TmBlocked + SOLVE_EQUATION_EULER_CPU(47); // ZETAS + SOLVE_EQUATION_EULER_CPU(48); // ZETAW +} + +void solve_forward_euler_cpu_adpt(real *sv, real stim_curr, real transmurality, real final_time, int sv_id, struct ode_solver *solver, real const *extra_params) { + + 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]; + + // Keep 'dt' inside the adaptive interval + if(*time_new + *dt > final_time) { + *dt = final_time - *time_new; + } + + RHS_cpu(sv, rDY, stim_curr, *dt, transmurality, extra_params); + *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, transmurality, extra_params); + *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 transmurality, real const *extra_params) { + + // Current modifiers + real INa_Multiplier = extra_params[0]; + real INaL_Multiplier = extra_params[1]; + real INaCa_Multiplier = extra_params[2]; + real INaK_Multiplier = extra_params[3]; + real INab_Multiplier = extra_params[4]; + real Ito_Multiplier = extra_params[5]; + real IKr_Multiplier = extra_params[6]; + real IKs_Multiplier = extra_params[7]; + real IK1_Multiplier = extra_params[8]; + real IKb_Multiplier = extra_params[9]; + real IKCa_Multiplier = extra_params[10]; + real ICaL_Multiplier = extra_params[11]; + real ICab_Multiplier = extra_params[12]; + real IpCa_Multiplier = extra_params[13]; + real ICaCl_Multiplier = extra_params[14]; + real IClb_Multiplier = extra_params[15]; + real Jrel_Multiplier = extra_params[16]; + real Jup_Multiplier = extra_params[17]; + real aCaMK_Multiplier = extra_params[18]; + real taurelp_Multiplier = extra_params[19]; + + // Get the celltype for the current cell + real celltype = transmurality; + + // Get the stimulus current from the current cell + real calc_I_stim = stim_current; + + // State variables (same order as the original Matlab script) + real v = sv[0]; + real nai = sv[1]; + real nass = sv[2]; + real ki = sv[3]; + real kss = sv[4]; + real cai = sv[5]; + real cass = sv[6]; + real cansr = sv[7]; + real cajsr = sv[8]; + real m = sv[9]; + real hp = sv[10]; + real h = sv[11]; + real j = sv[12]; + real jp = sv[13]; + real mL = sv[14]; + real hL = sv[15]; + real hLp = sv[16]; + real a = sv[17]; + real iF = sv[18]; + real iS = sv[19]; + real ap = sv[20]; + real iFp = sv[21]; + real iSp = sv[22]; + + // ical + real d = sv[23]; + real ff = sv[24]; + real fs = sv[25]; + real fcaf = sv[26]; + real fcas = sv[27]; + real jca = sv[28]; + real nca = sv[29]; + real nca_i = sv[30]; + real ffp = sv[31]; + real fcafp = sv[32]; + + real xs1 = sv[33]; + real xs2 = sv[34]; + real Jrel_np = sv[35]; + real CaMKt = sv[36]; + + // new MM ICaL states + real ikr_c0 = sv[37]; + real ikr_c1 = sv[38]; + real ikr_c2 = sv[39]; + real ikr_o = sv[40]; + real ikr_i = sv[41]; + real Jrel_p = sv[42]; + + const real cli = 24; // Intracellular Cl [mM] + const real clo = 150; // Extracellular Cl [mM] +// ----------------------------------------------------- + // Land-Niederer model + real XS = fmaxf(0,sv[43]); + real XW = fmaxf(0,sv[44]); + real Ca_TRPN = fmaxf(0,sv[45]); + real TmBlocked = sv[46]; + real ZETAS = sv[47]; + real ZETAW = sv[48]; + + #include "ToRORd_Land_mixed_endo_mid_epi.common.c" +} + +void RHS_RL_cpu(real *a_, real *b_, const real *sv, real *rDY_, real stim_current, real dt, \ + real transmurality, real const *extra_params) { + + // Current modifiers + real INa_Multiplier = extra_params[0]; + real INaL_Multiplier = extra_params[1]; + real INaCa_Multiplier = extra_params[2]; + real INaK_Multiplier = extra_params[3]; + real INab_Multiplier = extra_params[4]; + real Ito_Multiplier = extra_params[5]; + real IKr_Multiplier = extra_params[6]; + real IKs_Multiplier = extra_params[7]; + real IK1_Multiplier = extra_params[8]; + real IKb_Multiplier = extra_params[9]; + real IKCa_Multiplier = extra_params[10]; + real ICaL_Multiplier = extra_params[11]; + real ICab_Multiplier = extra_params[12]; + real IpCa_Multiplier = extra_params[13]; + real ICaCl_Multiplier = extra_params[14]; + real IClb_Multiplier = extra_params[15]; + real Jrel_Multiplier = extra_params[16]; + real Jup_Multiplier = extra_params[17]; + real aCaMK_Multiplier = extra_params[18]; + real taurelp_Multiplier = extra_params[19]; + + // Get the celltype for the current cell + real celltype = transmurality; + + // Get the stimulus current from the current cell + real calc_I_stim = stim_current; + + // State variables (same order as the original Matlab script) + real v = sv[0]; + real nai = sv[1]; + real nass = sv[2]; + real ki = sv[3]; + real kss = sv[4]; + real cai = sv[5]; + real cass = sv[6]; + real cansr = sv[7]; + real cajsr = sv[8]; + real m = sv[9]; + real hp = sv[10]; + real h = sv[11]; + real j = sv[12]; + real jp = sv[13]; + real mL = sv[14]; + real hL = sv[15]; + real hLp = sv[16]; + real a = sv[17]; + real iF = sv[18]; + real iS = sv[19]; + real ap = sv[20]; + real iFp = sv[21]; + real iSp = sv[22]; + + // ical + real d = sv[23]; + real ff = sv[24]; + real fs = sv[25]; + real fcaf = sv[26]; + real fcas = sv[27]; + real jca = sv[28]; + real nca = sv[29]; + real nca_i = sv[30]; + real ffp = sv[31]; + real fcafp = sv[32]; + + real xs1 = sv[33]; + real xs2 = sv[34]; + real Jrel_np = sv[35]; + real CaMKt = sv[36]; + + // new MM ICaL states + real ikr_c0 = sv[37]; + real ikr_c1 = sv[38]; + real ikr_c2 = sv[39]; + real ikr_o = sv[40]; + real ikr_i = sv[41]; + real Jrel_p = sv[42]; + + const real cli = 24; // Intracellular Cl [mM] + const real clo = 150; // Extracellular Cl [mM] +// ----------------------------------------------------- + real XS = fmaxf(0,sv[43]); + real XW = fmaxf(0,sv[44]); + real Ca_TRPN = fmaxf(0,sv[45]); + real TmBlocked = sv[46]; + real ZETAS = sv[47]; + real ZETAW = sv[48]; + + #include "ToRORd_Land_mixed_endo_mid_epi_RL.common.c" +} diff --git a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.common.c b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.common.c new file mode 100644 index 00000000..065a8c08 --- /dev/null +++ b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.common.c @@ -0,0 +1,830 @@ +// Avoid NaN errors +//if (v >= 100.0) v = 100.0; +//if (v <= -100.0) v = -100.0; + +//real aux, aux2, aux3; + +// Changeable parameters +real nao = 140.0; +real cao = 1.8; +real ko = 5.0; +// Localization of ICaL and NCX: the fraction in junctional subspace +real ICaL_fractionSS = 0.8; +real INaCa_fractionSS = 0.35; +// Additional parameters +real ca50 = 0.805; +real trpnmax = 0.07; + +// INPUT CODE: +int mode = 0; // 0 = "intact", 1 = "skinned" +real lambda = 1.0; +real lambda_rate = 0.0; + +// EC parameters +real perm50 = 0.35; +real TRPN_n = 2; +real koff = 0.1; +real dr = 0.25; +real wfrac = 0.5; +real TOT_A = 25; +real ktm_unblock = 0.021; +real beta_1 = -2.4; +real beta_0 = 2.3; +real gamma = 0.0085; +real gamma_wu = 0.615; +real phi = 2.23; + +real nperm; +real Tref; +real nu; +real mu; +if (mode == 1) { + nperm = 2.2; + //ca50=2.5; + Tref = 40.5; + nu = 1; + mu = 1; +} +else { + nperm = 2.036; + //ca50=0.805; + //ca50 = 0.5; + //ca50 = 0.7; + Tref = 120; + nu = 7; + mu = 3; +} + +real k_ws = 0.004; +real k_uw = 0.026; + +real lambda_min = 0.87; +real lambda_max = 1.2; + +//real dydt[6] = {0,0,0,0,0,0}; + +k_ws = k_ws*mu; +k_uw = k_uw*nu; + +real cdw = phi*k_uw*(1-dr)*(1-wfrac)/((1-dr)*wfrac); +real cds = phi*k_ws*(1-dr)*wfrac/dr; +real k_wu = k_uw*(1/wfrac-1)-k_ws; +real k_su = k_ws*(1/dr-1)*wfrac; +real A = (0.25*TOT_A)/((1-dr)*wfrac+dr)*(dr/0.25); + +//real lambda0 = (lambda_max < lambda) ? lambda_max : lambda; +real lambda0 = fminf(lambda_max, lambda); +//real lambda_aux = (lambda_min < lambda0) ? lambda_min : lambda0; +//aux = 1+beta_0*(lambda0+lambda_aux-(1+lambda_min)); +//real Lfac = (aux > 0) ? aux : 0; +real Lfac = fmaxf(0,1+beta_0*(lambda0+fminf(lambda_min,lambda0)-(1+lambda_min))); + +real XU = (1-TmBlocked)-XW-XS; // unattached available xb = all - tm blocked - already prepowerstroke - already post-poststroke - no overlap +real xb_ws = k_ws*XW; +real xb_uw = k_uw*XU; +real xb_wu = k_wu*XW; +real xb_su = k_su*XS; + +//aux = ((ZETAS>0)*ZETAS > (ZETAS<-1)*(-ZETAS-1)) ? (ZETAS>0)*ZETAS : (ZETAS<-1)*(-ZETAS-1); +real gamma_rate=gamma*fmaxf((ZETAS>0)*ZETAS,(ZETAS<-1)*(-ZETAS-1)); +real xb_su_gamma=gamma_rate*XS; +real gamma_rate_w=gamma_wu*abs(ZETAW); // weak xbs don't like being strained +real xb_wu_gamma=gamma_rate_w*XW; + +//dydt[0] = xb_ws-xb_su-xb_su_gamma; +//dydt[1] = xb_uw-xb_wu-xb_ws-xb_wu_gamma; +real dXS = xb_ws-xb_su-xb_su_gamma; +real dXW = xb_uw-xb_wu-xb_ws-xb_wu_gamma; + +//aux = (lambda-1 < 0.2) ? (lambda-1) : 0.2; +//ca50 = ca50+beta_1*aux; +ca50=ca50+beta_1*fminf(0.2,lambda-1); +//dydt[2] = koff*(pow(((cai*1000)/ca50),TRPN_n)*(1-Ca_TRPN)-Ca_TRPN); // untouched +real dCa_TRPN = koff*(pow(((cai*1000)/ca50),TRPN_n)*(1-Ca_TRPN)-Ca_TRPN); + +real XSSS = dr*0.5; +real XWSS = (1-dr)*wfrac*0.5; +real ktm_block = ktm_unblock*(pow(perm50,nperm))*0.5/(0.5-XSSS-XWSS); + +//aux = pow(Ca_TRPN,-(nperm/2)); +//aux2 = pow(Ca_TRPN,(nperm/2)); +//aux3 = (100 < aux) ? 100 : aux; +//dydt[3] = ktm_block*aux3*XU-ktm_unblock*aux2*TmBlocked; +//real dTmBlocked = ktm_block*aux3*XU-ktm_unblock*aux2*TmBlocked; +real dTmBlocked = ktm_block*fminf(100, pow(Ca_TRPN,-(nperm/2)))*XU - ktm_unblock*( pow(Ca_TRPN,(nperm/2)))*TmBlocked; + +// velocity dependence -- assumes distortion resets on W->S +//dydt[4] = A*lambda_rate-cds*ZETAS; // - gamma_rate * ZETAS; +//dydt[5] = A*lambda_rate-cdw*ZETAW; // - gamma_rate_w * ZETAW; +real dZETAS = A*lambda_rate-cds*ZETAS; // - gamma_rate * ZETAS; +real dZETAW = A*lambda_rate-cdw*ZETAW; // - gamma_rate_w * ZETAW; + +// Active Force +real Ta = Lfac*(Tref/dr)*((ZETAS+1)*XS+(ZETAW)*XW); + +// physical constants +real R=8314.0; +real T=310.0; +real F=96485.0; + +// cell geometry +real L=0.01; +real rad=0.0011; +real vcell=1000*3.14*rad*rad*L; +real Ageo=2*3.14*rad*rad+2*3.14*rad*L; +real Acap=2*Ageo; +real vmyo=0.68*vcell; +real vnsr=0.0552*vcell; +real vjsr=0.0048*vcell; +real vss=0.02*vcell; + +real fkatp = 0.0; +real gkatp = 4.3195; + +// CaMK constants +real KmCaMK=0.15; +real aCaMK=0.05*aCaMK_Multiplier; +real bCaMK=0.00068; +real CaMKo=0.05; +real KmCaM=0.0015; +// update CaMK +real CaMKb=CaMKo*(1.0-CaMKt)/(1.0+KmCaM/cass); +real CaMKa=CaMKb+CaMKt; +real dCaMKt=aCaMK*CaMKb*(CaMKb+CaMKt)-bCaMK*CaMKt; // Euler + +// reversal potentials +real ENa=(R*T/F)*log(nao/nai); +real EK=(R*T/F)*log(ko/ki); +real PKNa=0.01833; +real EKs=(R*T/F)*log((ko+PKNa*nao)/(ki+PKNa*nai)); + +// convenient shorthand calculations +real vffrt=v*F*F/(R*T); +real vfrt=v*F/(R*T); +real frt = F/(R*T); + +real K_o_n = 5.0; +real A_atp = 2.0; +real K_atp = 0.25; +real akik = pow((ko / K_o_n), 0.24); +real bkik = (1.0 / (1.0 + pow((A_atp / K_atp), 2.0))); + +real fINap=(1.0/(1.0+KmCaMK/CaMKa)); +real fINaLp=(1.0/(1.0+KmCaMK/CaMKa)); +real fItop=(1.0/(1.0+KmCaMK/CaMKa)); +real fICaLp=(1.0/(1.0+KmCaMK/CaMKa)); + +// INa formulations +// The Grandi implementation updated with INa phosphorylation. +// m gate +real mss = 1 / (pow(1 + exp( -(56.86 + v) / 9.03 ),2)); +real taum = 0.1292 * exp(-pow((v+45.79)/15.54,2)) + 0.06487 * exp(-pow((v-4.823)/51.12,2)); +real dm = (mss - m) / taum; // Rush-Larsen + +// h gate +real ah = (v >= -40) ? (0) : (0.057 * exp( -(v + 80) / 6.8 )); +real bh = (v >= -40) ? (0.77 / (0.13*(1 + exp( -(v + 10.66) / 11.1 )))) : ((2.7 * exp( 0.079 * v) + 3.1*pow(10,5) * exp(0.3485 * v))); +real tauh = 1 / (ah + bh); +real hss = 1 / (pow(1 + exp( (v + 71.55)/7.43 ),2)); +real dh = (hss - h) / tauh; // Rush-Larsen +// j gate +real aj = (v >= -40) ? (0) : (((-2.5428 * pow(10,4)*exp(0.2444*v) - 6.948*pow(10,-6) * exp(-0.04391*v)) * (v + 37.78)) / (1 + exp( 0.311 * (v + 79.23) ))); +real bj = (v >= -40) ? ((0.6 * exp( 0.057 * v)) / (1 + exp( -0.1 * (v + 32) ))) : ((0.02424 * exp( -0.01052 * v )) / (1 + exp( -0.1378 * (v + 40.14) ))); +real tauj = 1 / (aj + bj); +real jss = 1 / pow((1 + exp( (v + 71.55)/7.43 )),2); +real dj = (jss - j) / tauj; // Rush-Larsen + +// h phosphorylated +real hssp = 1 / pow((1 + exp( (v + 71.55 + 6)/7.43 )),2); +real dhp = (hssp - hp) / tauh; // Rush-Larsen +// j phosphorylated +real taujp = 1.46 * tauj; +real djp = (jss - jp) / taujp; // Rush-Larsen +real GNa = 11.7802; +real INa = INa_Multiplier*GNa*(v-ENa)*pow(m,3.0)*((1.0-fINap)*h*j+fINap*hp*jp); + +// INaL +// calculate INaL +real mLss=1.0/(1.0+exp((-(v+42.85))/5.264)); +real tm = 0.1292 * exp(-pow(((v+45.79)/15.54),2)) + 0.06487 * exp(-pow(((v-4.823)/51.12),2)); +real tmL=tm; +real dmL=(mLss-mL)/tmL; // Rush-Larsen +real hLss=1.0/(1.0+exp((v+87.61)/7.488)); +real thL=200.0; +real dhL=(hLss-hL)/thL; // Rush-Larsen +real hLssp=1.0/(1.0+exp((v+93.81)/7.488)); +real thLp=3.0*thL; +real dhLp=(hLssp-hLp)/thLp; // Rush-Larsen +real GNaL=0.0279*INaL_Multiplier; +if (celltype==EPI) GNaL=GNaL*0.6; +real INaL = GNaL*(v-ENa)*mL*((1.0-fINaLp)*hL+fINaLp*hLp); + +// ITo +// calculate Ito +real ass=1.0/(1.0+exp((-(v-14.34))/14.82)); +real ta=1.0515/(1.0/(1.2089*(1.0+exp(-(v-18.4099)/29.3814)))+3.5/(1.0+exp((v+100.0)/29.3814))); +real da=(ass-a)/ta; // Rush-Larsen +real iss=1.0/(1.0+exp((v+43.94)/5.711)); +real delta_epi = (celltype == EPI) ? 1.0-(0.95/(1.0+exp((v+70.0)/5.0))) : 1.0; +real tiF=4.562+1/(0.3933*exp((-(v+100.0))/100.0)+0.08004*exp((v+50.0)/16.59)); +real tiS=23.62+1/(0.001416*exp((-(v+96.52))/59.05)+1.780e-8*exp((v+114.1)/8.079)); +tiF=tiF*delta_epi; +tiS=tiS*delta_epi; +real AiF=1.0/(1.0+exp((v-213.6)/151.2)); +real AiS=1.0-AiF; +real diF=(iss-iF)/tiF; // Rush-Larsen +real diS=(iss-iS)/tiS; // Rush-Larsen +real i=AiF*iF+AiS*iS; +real assp=1.0/(1.0+exp((-(v-24.34))/14.82)); +real dap=(assp-ap)/ta; // Rush-Larsen +real dti_develop=1.354+1.0e-4/(exp((v-167.4)/15.89)+exp(-(v-12.23)/0.2154)); +real dti_recover=1.0-0.5/(1.0+exp((v+70.0)/20.0)); +real tiFp=dti_develop*dti_recover*tiF; +real tiSp=dti_develop*dti_recover*tiS; +real diFp=(iss-iFp)/tiFp; // Rush-Larsen +real diSp=(iss-iSp)/tiSp; // Rush-Larsen +real ip=AiF*iFp+AiS*iSp; +real Gto=0.16*Ito_Multiplier; +Gto = (celltype == EPI || celltype == MID) ? Gto*2.0 : Gto; +real Ito = Gto*(v-EK)*((1.0-fItop)*a*i+fItop*ap*ip); + +// ICaL +// a variant updated by jakub, using a changed activation curve +// it computes both ICaL in subspace and myoplasm (_i) + +// calculate ICaL, ICaNa, ICaK +real dss=1.0763*exp(-1.0070*exp(-0.0829*(v))); // magyar +if(v >31.4978) dss = 1; // activation cannot be greater than 1 +real td= 0.6+1.0/(exp(-0.05*(v+6.0))+exp(0.09*(v+14.0))); +real dd=(dss-d)/td; // Rush-Larsen +real fss=1.0/(1.0+exp((v+19.58)/3.696)); +real tff=7.0+1.0/(0.0045*exp(-(v+20.0)/10.0)+0.0045*exp((v+20.0)/10.0)); +real tfs=1000.0+1.0/(0.000035*exp(-(v+5.0)/4.0)+0.000035*exp((v+5.0)/6.0)); +real Aff=0.6; +real Afs=1.0-Aff; +real dff=(fss-ff)/tff; // Rush-Larsen +real dfs=(fss-fs)/tfs; // Rush-Larsen +real f=Aff*ff+Afs*fs; +real fcass=fss; +real tfcaf=7.0+1.0/(0.04*exp(-(v-4.0)/7.0)+0.04*exp((v-4.0)/7.0)); +real tfcas=100.0+1.0/(0.00012*exp(-v/3.0)+0.00012*exp(v/7.0)); + +real Afcaf=0.3+0.6/(1.0+exp((v-10.0)/10.0)); +real Afcas=1.0-Afcaf; +real dfcaf=(fcass-fcaf)/tfcaf; // Rush-Larsen +real dfcas=(fcass-fcas)/tfcas; // Rush-Larsen +real fca=Afcaf*fcaf+Afcas*fcas; + +real tjca = 75; +//real tjca = 72.5; +real jcass = 1.0/(1.0+exp((v+18.08)/(2.7916))); +real djca=(jcass-jca)/tjca; // Rush-Larsen +real tffp=2.5*tff; +real dffp=(fss-ffp)/tffp; // Rush-Larsen +real fp=Aff*ffp+Afs*fs; +real tfcafp=2.5*tfcaf; +real dfcafp=(fcass-fcafp)/tfcafp; // Rush-Larsen +real fcap=Afcaf*fcafp+Afcas*fcas; + +// SS nca +real Kmn=0.002; +real k2n=500.0; +real km2n=jca*1; +real anca=1.0/(k2n/km2n+pow((1.0+Kmn/cass),4.0)); +real dnca=anca*k2n-nca*km2n; // Euler + +// myoplasmic nca +real anca_i = 1.0/(k2n/km2n+pow((1.0+Kmn/cai),4.0)); +real dnca_i = anca_i*k2n-nca_i*km2n; // Euler + +// SS driving force +//real clo = 150; // Extracellular Cl [mM] +//real cli = 24.0; // Intracellular Cl [mM] +real Io = 0.5*(nao + ko + clo + 4*cao)/1000; // ionic strength outside. /1000 is for things being in micromolar +real Ii = 0.5*(nass + kss + cli + 4*cass)/1000; // ionic strength outside. /1000 is for things being in micromolar +//real Ii = 0.5*(nass + kss + clss + 4*cass)/1000; // (dynCl) ionic strength outside. /1000 is for things being in micromolar +// The ionic strength is too high for basic DebHuc. We'll use Davies +real dielConstant = 74; // water at 37 degrees. +real temp = 310; // body temp in kelvins. +real constA = 1.82*pow(10,6)*pow((dielConstant*temp),(-1.5)); + +real gamma_cai = exp(-constA * 4 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii)); +real gamma_cao = exp(-constA * 4 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io)); +real gamma_nai = exp(-constA * 1 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii)); +real gamma_nao = exp(-constA * 1 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io)); +real gamma_ki = exp(-constA * 1 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii)); +real gamma_kao = exp(-constA * 1 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io)); + +real PhiCaL_ss = 4.0*vffrt*(gamma_cai*cass*exp(2.0*vfrt)-gamma_cao*cao)/(exp(2.0*vfrt)-1.0); +real PhiCaNa_ss = 1.0*vffrt*(gamma_nai*nass*exp(1.0*vfrt)-gamma_nao*nao)/(exp(1.0*vfrt)-1.0); +real PhiCaK_ss = 1.0*vffrt*(gamma_ki*kss*exp(1.0*vfrt)-gamma_kao*ko)/(exp(1.0*vfrt)-1.0); + +// Myo driving force +Io = 0.5*(nao + ko + clo + 4*cao)/1000; // ionic strength outside. /1000 is for things being in micromolar +Ii = 0.5*(nai + ki + cli + 4*cai)/1000; // ionic strength outside. /1000 is for things being in micromolar +// The ionic strength is too high for basic DebHuc. We'll use Davies + +gamma_cai = exp(-constA * 4 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii)); +gamma_cao = exp(-constA * 4 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io)); +gamma_nai = exp(-constA * 1 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii)); +gamma_nao = exp(-constA * 1 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io)); +gamma_ki = exp(-constA * 1 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii)); +gamma_kao = exp(-constA * 1 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io)); + +real gammaCaoMyo = gamma_cao; +real gammaCaiMyo = gamma_cai; + +real PhiCaL_i = 4.0*vffrt*(gamma_cai*cai*exp(2.0*vfrt)-gamma_cao*cao)/(exp(2.0*vfrt)-1.0); +real PhiCaNa_i = 1.0*vffrt*(gamma_nai*nai*exp(1.0*vfrt)-gamma_nao*nao)/(exp(1.0*vfrt)-1.0); +real PhiCaK_i = 1.0*vffrt*(gamma_ki*ki*exp(1.0*vfrt)-gamma_kao*ko)/(exp(1.0*vfrt)-1.0); +// The rest +real PCa=8.3757e-05 * ICaL_Multiplier; +if (celltype==EPI) + PCa=PCa*1.2; +else if (celltype==MID) + //PCa=PCa*2; + PCa = PCa*1.8; + +real PCap=1.1*PCa; +real PCaNa=0.00125*PCa; +real PCaK=3.574e-4*PCa; +real PCaNap=0.00125*PCap; +real PCaKp=3.574e-4*PCap; + +real ICaL_ss=(1.0-fICaLp)*PCa*PhiCaL_ss*d*(f*(1.0-nca)+jca*fca*nca)+fICaLp*PCap*PhiCaL_ss*d*(fp*(1.0-nca)+jca*fcap*nca); +real ICaNa_ss=(1.0-fICaLp)*PCaNa*PhiCaNa_ss*d*(f*(1.0-nca)+jca*fca*nca)+fICaLp*PCaNap*PhiCaNa_ss*d*(fp*(1.0-nca)+jca*fcap*nca); +real ICaK_ss=(1.0-fICaLp)*PCaK*PhiCaK_ss*d*(f*(1.0-nca)+jca*fca*nca)+fICaLp*PCaKp*PhiCaK_ss*d*(fp*(1.0-nca)+jca*fcap*nca); + +real ICaL_i=(1.0-fICaLp)*PCa*PhiCaL_i*d*(f*(1.0-nca_i)+jca*fca*nca_i)+fICaLp*PCap*PhiCaL_i*d*(fp*(1.0-nca_i)+jca*fcap*nca_i); +real ICaNa_i=(1.0-fICaLp)*PCaNa*PhiCaNa_i*d*(f*(1.0-nca_i)+jca*fca*nca_i)+fICaLp*PCaNap*PhiCaNa_i*d*(fp*(1.0-nca_i)+jca*fcap*nca_i); +real ICaK_i=(1.0-fICaLp)*PCaK*PhiCaK_i*d*(f*(1.0-nca_i)+jca*fca*nca_i)+fICaLp*PCaKp*PhiCaK_i*d*(fp*(1.0-nca_i)+jca*fcap*nca_i); + +// And we weight ICaL (in ss) and ICaL_i +ICaL_i = ICaL_i * (1-ICaL_fractionSS); +ICaNa_i = ICaNa_i * (1-ICaL_fractionSS); +ICaK_i = ICaK_i * (1-ICaL_fractionSS); +ICaL_ss = ICaL_ss * ICaL_fractionSS; +ICaNa_ss = ICaNa_ss * ICaL_fractionSS; +ICaK_ss = ICaK_ss * ICaL_fractionSS; + +real ICaL = ICaL_ss + ICaL_i; +real ICaNa = ICaNa_ss + ICaNa_i; +real ICaK = ICaK_ss + ICaK_i; +real ICaL_tot = ICaL + ICaNa + ICaK; + +// Ikr +// Variant based on Lu-Vandenberg +// Extracting state vector +// IMPORTANT: For reason of backward compatibility of naming of an older version of a MM IKr +// c3 in code is c0 in article diagram, c2 is c1, c1 is c2 +real c0 = ikr_c0; +real c1 = ikr_c1; +real c2 = ikr_c2; +//real c0 = ikr_c2; +//real c1 = ikr_c1; +//real c2 = ikr_c0; +real o = ikr_o; +real I = ikr_i; +real b = 0; // no channels blocked in via the mechanism of specific MM states + +// transition rates +// from c0 to c1 in l-v model, +real alpha = 0.1161 * exp(0.2990 * vfrt); +// from c1 to c0 in l-v/ +real beta = 0.2442 * exp(-1.604 * vfrt); + +// from c1 to c2 in l-v/ +real alpha1 = 1.25 * 0.1235; +// from c2 to c1 in l-v/ +real beta1 = 0.1911; + +// from c2 to o/ c1 to o +real alpha2 =0.0578 * exp(0.9710 * vfrt); +// from o to c2/ +real beta2 = 0.349e-3* exp(-1.062 * vfrt); + +// from o to i +real alphai = 0.2533 * exp(0.5953 * vfrt); +// from i to o +real betai = 1.25* 0.0522 * exp(-0.8209 * vfrt); + +// from c2 to i (from c1 in orig) +real alphac2ToI = 0.52e-4 * exp(1.525 * vfrt); +// from i to c2 +real betaItoC2 = 0.85e-8 * exp(-1.842 * vfrt); +betaItoC2 = (beta2 * betai * alphac2ToI)/(alpha2 * alphai); +// transitions themselves +// for reason of backward compatibility of naming of an older version of a +// MM IKr, c3 in code is c0 in article diagram, c2 is c1, c1 is c2. + +real dc0 = c1 * beta - c0 * alpha; // Euler +real dc1 = c0 * alpha + c2*beta1 - c1*(beta+alpha1); // Euler +real dc2 = c1 * alpha1 + o*beta2 + I*betaItoC2 - c2 * (beta1 + alpha2 + alphac2ToI); // Euler +real delta_o = c2 * alpha2 + I*betai - o*(beta2+alphai); // Euler +real di = c2*alphac2ToI + o*alphai - I*(betaItoC2 + betai); // Euler + +real GKr = 0.0321 * sqrt(ko/5) * IKr_Multiplier; // 1st element compensates for change to ko (sqrt(5/5.4)* 0.0362) +if (celltype==EPI) + GKr=GKr*1.3; +else if (celltype==MID) + GKr=GKr*0.8; + +real IKr = GKr * o * (v-EK); + +// calculate IKs +real xs1ss=1.0/(1.0+exp((-(v+11.60))/8.932)); +real txs1=817.3+1.0/(2.326e-4*exp((v+48.28)/17.80)+0.001292*exp((-(v+210.0))/230.0)); +real dxs1=(xs1ss-xs1)/txs1; // Rush-Larsen +real xs2ss=xs1ss; +real txs2=1.0/(0.01*exp((v-50.0)/20.0)+0.0193*exp((-(v+66.54))/31.0)); +real dxs2=(xs2ss-xs2)/txs2; // Rush-Larsen +real KsCa=1.0+0.6/(1.0+pow((3.8e-5/cai),1.4)); +real GKs= 0.0011*IKs_Multiplier; +if (celltype==EPI) + GKs=GKs*1.4; +real IKs = GKs*KsCa*xs1*xs2*(v-EKs); + +// IK1 +real aK1 = 4.094/(1+exp(0.1217*(v-EK-49.934))); +real bK1 = (15.72*exp(0.0674*(v-EK-3.257))+exp(0.0618*(v-EK-594.31)))/(1+exp(-0.1629*(v-EK+14.207))); +real K1ss = aK1/(aK1+bK1); +real GK1=IK1_Multiplier * 0.6992; // 0.7266; // * sqrt(5/5.4)) +if (celltype==EPI) + GK1=GK1*1.2; +else if (celltype==MID) + GK1=GK1*1.3; +real IK1=GK1*sqrt(ko/5)*K1ss*(v-EK); + +// IKCa +real fIKCass = 0.8; +real kdikca = 6.05e-4; +real ikcan = 3.5; +real GKCa = 0.003 * IKCa_Multiplier; +real IKCa_ss = GKCa * fIKCass * pow(cass,ikcan) / (pow(cass,ikcan) + pow(kdikca,ikcan)) * (v-EK); +real IKCa_i = GKCa * (1.0-fIKCass) * pow(cai,ikcan) / (pow(cai,ikcan) + pow(kdikca,ikcan)) * (v-EK); +real IKCa = IKCa_ss + IKCa_i; + +// INaCa +real zca = 2.0; +real kna1=15.0; +real kna2=5.0; +real kna3=88.12; +real kasymm=12.5; +real wna=6.0e4; +real wca=6.0e4; +real wnaca=5.0e3; +real kcaon=1.5e6; +real kcaoff=5.0e3; +real qna=0.5224; +real qca=0.1670; +real hca=exp((qca*v*F)/(R*T)); +real hna=exp((qna*v*F)/(R*T)); + +real h1=1+nai/kna3*(1+hna); +real h2=(nai*hna)/(kna3*h1); +real h3=1.0/h1; +real h4=1.0+nai/kna1*(1+nai/kna2); +real h5=nai*nai/(h4*kna1*kna2); +real h6=1.0/h4; +real h7=1.0+nao/kna3*(1.0+1.0/hna); +real h8=nao/(kna3*hna*h7); +real h9=1.0/h7; +real h10=kasymm+1.0+nao/kna1*(1.0+nao/kna2); +real h11=nao*nao/(h10*kna1*kna2); +real h12=1.0/h10; +real k1=h12*cao*kcaon; +real k2=kcaoff; +real k3p=h9*wca; +real k3pp=h8*wnaca; +real k3=k3p+k3pp; +real k4p=h3*wca/hca; +real k4pp=h2*wnaca; +real k4=k4p+k4pp; +real k5=kcaoff; +real k6=h6*cai*kcaon; +real k7=h5*h2*wna; +real k8=h8*h11*wna; +real x1=k2*k4*(k7+k6)+k5*k7*(k2+k3); +real x2=k1*k7*(k4+k5)+k4*k6*(k1+k8); +real x3=k1*k3*(k7+k6)+k8*k6*(k2+k3); +real x4=k2*k8*(k4+k5)+k3*k5*(k1+k8); +real E1=x1/(x1+x2+x3+x4); +real E2=x2/(x1+x2+x3+x4); +real E3=x3/(x1+x2+x3+x4); +real E4=x4/(x1+x2+x3+x4); +real KmCaAct=150.0e-6; +real allo=1.0/(1.0+pow((KmCaAct/cai),2.0)); +real zna=1.0; +real JncxNa=3.0*(E4*k7-E1*k8)+E3*k4pp-E2*k3pp; +real JncxCa=E2*k2-E1*k1; +real Gncx= 0.0034 * INaCa_Multiplier; +if (celltype==EPI) + Gncx=Gncx*1.1; +else if (celltype==MID) + Gncx=Gncx*1.4; +real INaCa_i = (1-INaCa_fractionSS)*Gncx*allo*(zna*JncxNa+zca*JncxCa); + +// calculate INaCa_ss +h1=1+nass/kna3*(1+hna); +h2=(nass*hna)/(kna3*h1); +h3=1.0/h1; +h4=1.0+nass/kna1*(1+nass/kna2); +h5=nass*nass/(h4*kna1*kna2); +h6=1.0/h4; +h7=1.0+nao/kna3*(1.0+1.0/hna); +h8=nao/(kna3*hna*h7); +h9=1.0/h7; +h10=kasymm+1.0+nao/kna1*(1+nao/kna2); +h11=nao*nao/(h10*kna1*kna2); +h12=1.0/h10; +k1=h12*cao*kcaon; +k2=kcaoff; +k3p=h9*wca; +k3pp=h8*wnaca; +k3=k3p+k3pp; +k4p=h3*wca/hca; +k4pp=h2*wnaca; +k4=k4p+k4pp; +k5=kcaoff; +k6=h6*cass*kcaon; +k7=h5*h2*wna; +k8=h8*h11*wna; +x1=k2*k4*(k7+k6)+k5*k7*(k2+k3); +x2=k1*k7*(k4+k5)+k4*k6*(k1+k8); +x3=k1*k3*(k7+k6)+k8*k6*(k2+k3); +x4=k2*k8*(k4+k5)+k3*k5*(k1+k8); +E1=x1/(x1+x2+x3+x4); +E2=x2/(x1+x2+x3+x4); +E3=x3/(x1+x2+x3+x4); +E4=x4/(x1+x2+x3+x4); +KmCaAct=150.0e-6 ; +allo=1.0/(1.0+pow((KmCaAct/cass),2.0)); +JncxNa=3.0*(E4*k7-E1*k8)+E3*k4pp-E2*k3pp; +JncxCa=E2*k2-E1*k1; +real INaCa_ss = INaCa_fractionSS*Gncx*allo*(zna*JncxNa+zca*JncxCa); + +// calculate INaK +real k1p=949.5; +real k1m=182.4; +real k2p=687.2; +real k2m=39.4; +k3p=1899.0; +real k3m=79300.0; +k4p=639.0; +real k4m=40.0; +real Knai0=9.073; +real Knao0=27.78; +real delta=-0.1550; +real Knai=Knai0*exp((delta*v*F)/(3.0*R*T)); +real Knao=Knao0*exp(((1.0-delta)*v*F)/(3.0*R*T)); +real Kki=0.5; +real Kko=0.3582; +real MgADP=0.05; +real MgATP=9.8; +real Kmgatp=1.698e-7; +real H=1.0e-7; +real eP=4.2; +real Khp=1.698e-7; +real Knap=224.0; +real Kxkur=292.0; +real P=eP/(1.0+H/Khp+nai/Knap+ki/Kxkur); +real a1=(k1p*pow((nai/Knai),3.0))/(pow((1.0+nai/Knai),3.0)+pow((1.0+ki/Kki),2.0)-1.0); +real b1=k1m*MgADP; +real a2=k2p; +real b2=(k2m*pow((nao/Knao),3.0))/(pow((1.0+nao/Knao),3.0)+pow((1.0+ko/Kko),2.0)-1.0); +real a3=(k3p*pow((ko/Kko),2.0))/(pow((1.0+nao/Knao),3.0)+pow((1.0+ko/Kko),2.0)-1.0); +real b3=(k3m*P*H)/(1.0+MgATP/Kmgatp); +real a4=(k4p*MgATP/Kmgatp)/(1.0+MgATP/Kmgatp); +real b4=(k4m*pow((ki/Kki),2.0))/(pow((1.0+nai/Knai),3.0)+pow((1.0+ki/Kki),2.0)-1.0); +x1=a4*a1*a2+b2*b4*b3+a2*b4*b3+b3*a1*a2; +x2=b2*b1*b4+a1*a2*a3+a3*b1*b4+a2*a3*b4; +x3=a2*a3*a4+b3*b2*b1+b2*b1*a4+a3*a4*b1; +x4=b4*b3*b2+a3*a4*a1+b2*a4*a1+b3*b2*a1; +E1=x1/(x1+x2+x3+x4); +E2=x2/(x1+x2+x3+x4); +E3=x3/(x1+x2+x3+x4); +E4=x4/(x1+x2+x3+x4); +real zk=1.0; +real JnakNa=3.0*(E1*a3-E2*b3); +real JnakK=2.0*(E4*b1-E3*a1); +real Pnak= 15.4509; +if (celltype==EPI) + Pnak=Pnak*0.9; +else if (celltype==MID) + Pnak=Pnak*0.7; + +real INaK = Pnak*(zna*JnakNa+zk*JnakK)*INaK_Multiplier; + +// Minor/background currents +// calculate IKb +real xkb=1.0/(1.0+exp(-(v-10.8968)/(23.9871))); +real GKb=0.0189*IKb_Multiplier; + +if (IKCa_Multiplier > 0.0) + GKb = GKb*0.9; +if (celltype==EPI) + GKb=GKb*0.6; +real IKb = GKb*xkb*(v-EK); + +// calculate INab +real PNab=1.9239e-09*INab_Multiplier; +real INab=PNab*vffrt*(nai*exp(vfrt)-nao)/(exp(vfrt)-1.0); + +// calculate ICab +real PCab=5.9194e-08*ICab_Multiplier; +real ICab=PCab*4.0*vffrt*(gammaCaiMyo*cai*exp(2.0*vfrt)-gammaCaoMyo*cao)/(exp(2.0*vfrt)-1.0); + +// calculate IpCa +real GpCa=5e-04*IpCa_Multiplier; +real IpCa=GpCa*cai/(0.0005+cai); + +// Chloride +// I_ClCa: Ca-activated Cl Current, I_Clbk: background Cl Current + +real ECl = (R*T/F)*log(cli/clo); // [mV] +//real EClss = (R*T/F)*log(clss/clo); // [mV] dynCl + +real Fjunc = 1; +real Fsl = 1-Fjunc; // fraction in SS and in myoplasm - as per literature, I(Ca)Cl is in junctional subspace + +real GClCa = ICaCl_Multiplier * 0.2843; // [mS/uF] +real GClB = IClb_Multiplier * 1.98e-3; // [mS/uF] +real KdClCa = 0.1; // [mM] + +//real I_ClCa_junc = Fjunc*GClCa/(1+KdClCa/cass)*(v-EClss); // dynCl +//real I_ClCa_sl = Fsl*GClCa/(1+KdClCa/cai)*(v-ECl); // dynCl +real I_ClCa_junc = Fjunc*GClCa/(1+KdClCa/cass)*(v-ECl); +real I_ClCa_sl = Fsl*GClCa/(1+KdClCa/cai)*(v-ECl); + +real I_ClCa = I_ClCa_junc+I_ClCa_sl; +real I_Clbk = GClB*(v-ECl); + +// Calcium handling +// calculate ryanodione receptor calcium induced calcium release from the jsr +real fJrelp=(1.0/(1.0+KmCaMK/CaMKa)); + +// Jrel +real jsrMidpoint = 1.7; + +real bt=4.75; +real a_rel=0.5*bt; +real Jrel_inf=a_rel*(-ICaL)/(1.0+pow((jsrMidpoint/cajsr),8.0)); +if (celltype==MID) + Jrel_inf=Jrel_inf*1.7; +real tau_rel=bt/(1.0+0.0123/cajsr); + +if (tau_rel<0.001) + tau_rel=0.001; + +real dJrelnp=(Jrel_inf-Jrel_np)/tau_rel; // Rush-Larsen + +real btp=1.25*bt; +real a_relp=0.5*btp; +real Jrel_infp=a_relp*(-ICaL)/(1.0+pow((jsrMidpoint/cajsr),8.0)); +if (celltype==MID) + Jrel_infp=Jrel_infp*1.7; +real tau_relp=btp/(1.0+0.0123/cajsr); + +if (tau_relp<0.001) + tau_relp=0.001; + +tau_relp = tau_relp*taurelp_Multiplier; +real dJrelp=(Jrel_infp-Jrel_p)/tau_relp; // Rush-Larsen +real Jrel = Jrel_Multiplier * 1.5378 * ((1.0-fJrelp)*Jrel_np+fJrelp*Jrel_p); + +real fJupp=(1.0/(1.0+KmCaMK/CaMKa)); + +// Jup +// calculate serca pump, ca uptake flux +// camkFactor = 2.4; +// gjup = 0.00696; +// Jupnp=Jup_Multiplier * gjup*cai/(cai+0.001); +// Jupp=Jup_Multiplier * camkFactor*gjup*cai/(cai + 8.2500e-04); +// if celltype==1 +// Jupnp=Jupnp*1.3; +// Jupp=Jupp*1.3; +// end +// +// +// Jleak=Jup_Multiplier * 0.00629 * cansr/15.0; +// Jup=(1.0-fJupp)*Jupnp+fJupp*Jupp-Jleak; + +// calculate serca pump, ca uptake flux +real Jupnp = Jup_Multiplier * 0.005425*cai/(cai+0.00092); +real Jupp = Jup_Multiplier * 2.75*0.005425*cai/(cai+0.00092-0.00017); +if (celltype==EPI) { + Jupnp=Jupnp*1.3; + Jupp=Jupp*1.3; +} +real Jleak=Jup_Multiplier* 0.0048825*cansr/15.0; +real Jup=(1.0-fJupp)*Jupnp+fJupp*Jupp-Jleak; + +//calculate tranlocation flux +real Jtr=(cansr-cajsr)/60; + +// I_katp current (fkatp current) +//real I_katp = (fkatp * gkatp * akik * bkik * (v - EK)); + +// stimulus current +real Istim = calc_I_stim; + +//update the membrane voltage +//real dv=-(INa+INaL+Ito+ICaL+ICaNa+ICaK+IKr+IKs+IK1+INaCa_i+INaCa_ss+INaK+INab+IKb+IpCa+ICab+I_ClCa+I_Clbk+Istim); // Euler +real dv=-(INa+INaL+Ito+ICaL+ICaNa+ICaK+IKr+IKs+IK1+IKCa+INaCa_i+INaCa_ss+INaK+INab+IKb+IpCa+ICab+I_ClCa+I_Clbk+Istim); // Euler +//real dv=-(INa+INaL+Ito+ICaL+ICaNa+ICaK+IKr+IKs+IK1+INaCa_i+INaCa_ss+INaK+INab+IKb+IpCa+ICab+I_ClCa+I_Clbk+I_katp+Istim); // Euler + +// calculate diffusion fluxes +real JdiffNa=(nass-nai)/2.0; +real JdiffK=(kss-ki)/2.0; +real Jdiff=(cass-cai)/0.2; +//real JdiffCl=(clss-cli)/2.0; // dynCl + +// calcium buffer constants +real cmdnmax= 0.05; +if (celltype==EPI) + cmdnmax=cmdnmax*1.3; +real kmcmdn=0.00238; +//real trpnmax=0.07; +real kmtrpn=0.0005; +real BSRmax=0.047; +real KmBSR = 0.00087; +real BSLmax=1.124; +real KmBSL = 0.0087; +real csqnmax=10.0; +real kmcsqn=0.8; + +// update intracellular concentrations, using buffers for cai, cass, cajsr +real dnai=-(ICaNa_i+INa+INaL+3.0*INaCa_i+3.0*INaK+INab)*Acap/(F*vmyo)+JdiffNa*vss/vmyo; // Euler +real dnass=-(ICaNa_ss+3.0*INaCa_ss)*Acap/(F*vss)-JdiffNa; // Euler + +real dki=-(ICaK_i+Ito+IKr+IKs+IK1+IKb+Istim-2.0*INaK)*Acap/(F*vmyo)+JdiffK*vss/vmyo; // Euler +real dkss=-(ICaK_ss)*Acap/(F*vss)-JdiffK; // Euler + +//real Bcai=1.0/(1.0+cmdnmax*kmcmdn/pow((kmcmdn+cai),2.0)+trpnmax*kmtrpn/pow((kmtrpn+cai),2.0)); // dynCl +real Bcai=1.0/(1.0+cmdnmax*kmcmdn/pow((kmcmdn+cai),2.0)); +//real dcai=Bcai*(-(ICaL_i + IpCa+ICab-2.0*INaCa_i)*Acap/(2.0*F*vmyo)-Jup*vnsr/vmyo+Jdiff*vss/vmyo); // dynCl +real dcai=Bcai*(-(ICaL_i + IpCa+ICab-2.0*INaCa_i)*Acap/(2.0*F*vmyo)-Jup*vnsr/vmyo+Jdiff*vss/vmyo-dCa_TRPN*trpnmax); + +real Bcass=1.0/(1.0+BSRmax*KmBSR/pow((KmBSR+cass),2.0)+BSLmax*KmBSL/pow((KmBSL+cass),2.0)); +real dcass=Bcass*(-(ICaL_ss-2.0*INaCa_ss)*Acap/(2.0*F*vss)+Jrel*vjsr/vss-Jdiff); + +real dcansr=Jup-Jtr*vjsr/vnsr; + +real Bcajsr=1.0/(1.0+csqnmax*kmcsqn/pow((kmcsqn+cajsr),2.0)); +real dcajsr=Bcajsr*(Jtr-Jrel); + +//real dcli = - (I_Clbk + I_ClCa_sl)*Acap/(-1*F*vmyo)+JdiffCl*vss/vmyo; // dynCl +//real dclss = - I_ClCa_junc*Acap/(-1*F*vss)-JdiffCl; // dynCl + +// Right-hand side +rDY_[0] = dv; +rDY_[1] = dnai; +rDY_[2] = dnass; +rDY_[3] = dki; +rDY_[4] = dkss; +rDY_[5] = dcai; +rDY_[6] = dcass; +rDY_[7] = dcansr; +rDY_[8] = dcajsr; +rDY_[9] = dm; +rDY_[10] = dhp; +rDY_[11] = dh; +rDY_[12] = dj; +rDY_[13] = djp; +rDY_[14] = dmL; +rDY_[15] = dhL; +rDY_[16] = dhLp; +rDY_[17] = da; +rDY_[18] = diF; +rDY_[19] = diS; +rDY_[20] = dap; +rDY_[21] = diFp; +rDY_[22] = diSp; +rDY_[23] = dd; +rDY_[24] = dff; +rDY_[25] = dfs; +rDY_[26] = dfcaf; +rDY_[27] = dfcas; +rDY_[28] = djca; +rDY_[29] = dnca; +rDY_[30] = dnca_i; +rDY_[31] = dffp; +rDY_[32] = dfcafp; +rDY_[33] = dxs1; +rDY_[34] = dxs2; +rDY_[35] = dJrelnp; +rDY_[36] = dCaMKt; +rDY_[37] = dc0; +rDY_[38] = dc1; +rDY_[39] = dc2; +rDY_[40] = delta_o; +rDY_[41] = di; +rDY_[42] = dJrelp; +// ----------------------- +// Land-Niederer +rDY_[43] = dXS; +rDY_[44] = dXW; +rDY_[45] = dCa_TRPN; +rDY_[46] = dTmBlocked; +rDY_[47] = dZETAS; +rDY_[48] = dZETAW; diff --git a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.cu b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.cu new file mode 100644 index 00000000..5b0f64c7 --- /dev/null +++ b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.cu @@ -0,0 +1,974 @@ +#include "ToRORd_Land_mixed_endo_mid_epi.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) { + int threadID = blockDim.x * blockIdx.x + threadIdx.x; + + if (threadID < num_volumes) { + for (int i = 0; i < NEQ; i++) { + // Default initial conditions for ENDO cell (from original Matlab script) + *((real * )((char *) sv + pitch * 0) + threadID) = -8.863699e+01; + *((real * )((char *) sv + pitch * 1) + threadID) = 1.189734e+01; + *((real * )((char *) sv + pitch * 2) + threadID) = 1.189766e+01; + *((real * )((char *) sv + pitch * 3) + threadID) = 1.412345e+02; + *((real * )((char *) sv + pitch * 4) + threadID) = 1.412344e+02; + *((real * )((char *) sv + pitch * 5) + threadID) = 7.267473e-05; + *((real * )((char *) sv + pitch * 6) + threadID) = 6.337870e-05; + *((real * )((char *) sv + pitch * 7) + threadID) = 1.532653e+00; + *((real * )((char *) sv + pitch * 8) + threadID) = 1.533946e+00; + *((real * )((char *) sv + pitch * 9) + threadID) = 8.280078e-04; + *((real * )((char *) sv + pitch * 10) + threadID) = 6.665272e-01; + *((real * )((char *) sv + pitch * 11) + threadID) = 8.260208e-01; + *((real * )((char *) sv + pitch * 12) + threadID) = 8.260560e-01; + *((real * )((char *) sv + pitch * 13) + threadID) = 8.258509e-01; + *((real * )((char *) sv + pitch * 14) + threadID) = 1.668686e-04; + *((real * )((char *) sv + pitch * 15) + threadID) = 5.228306e-01; + *((real * )((char *) sv + pitch * 16) + threadID) = 2.859696e-01; + *((real * )((char *) sv + pitch * 17) + threadID) = 9.591370e-04; + *((real * )((char *) sv + pitch * 18) + threadID) = 9.996012e-01; + *((real * )((char *) sv + pitch * 19) + threadID) = 5.934016e-01; + *((real * )((char *) sv + pitch * 20) + threadID) = 4.886961e-04; + *((real * )((char *) sv + pitch * 21) + threadID) = 9.996011e-01; + *((real * )((char *) sv + pitch * 22) + threadID) = 6.546687e-01; + *((real * )((char *) sv + pitch * 23) + threadID) = 9.500075e-32; + *((real * )((char *) sv + pitch * 24) + threadID) = 1.000000e+00; + *((real * )((char *) sv + pitch * 25) + threadID) = 9.392580e-01; + *((real * )((char *) sv + pitch * 26) + threadID) = 1.000000e+00; + *((real * )((char *) sv + pitch * 27) + threadID) = 9.998984e-01; + *((real * )((char *) sv + pitch * 28) + threadID) = 9.999783e-01; + *((real * )((char *) sv + pitch * 29) + threadID) = 4.448162e-04; + *((real * )((char *) sv + pitch * 30) + threadID) = 7.550725e-04; + *((real * )((char *) sv + pitch * 31) + threadID) = 1.000000e+00; + *((real * )((char *) sv + pitch * 32) + threadID) = 1.000000e+00; + *((real * )((char *) sv + pitch * 33) + threadID) = 2.424047e-01; + *((real * )((char *) sv + pitch * 34) + threadID) = 1.795377e-04; + *((real * )((char *) sv + pitch * 35) + threadID) = -6.883086e-25; + *((real * )((char *) sv + pitch * 36) + threadID) = 1.117498e-02; + *((real * )((char *) sv + pitch * 37) + threadID) = 9.980366e-01; + *((real * )((char *) sv + pitch * 38) + threadID) = 8.588018e-04; + *((real * )((char *) sv + pitch * 39) + threadID) = 7.097447e-04; + *((real * )((char *) sv + pitch * 40) + threadID) = 3.812617e-04; + *((real * )((char *) sv + pitch * 41) + threadID) = 1.357116e-05; + *((real * )((char *) sv + pitch * 42) + threadID) = 2.302525e-23; + *((real * )((char *) sv + pitch * 43) + threadID) = 1.561941e-04; + *((real * )((char *) sv + pitch * 44) + threadID) = 2.351289e-04; + *((real * )((char *) sv + pitch * 45) + threadID) = 8.077631e-03; + *((real * )((char *) sv + pitch * 46) + threadID) = 9.993734e-01; + *((real * )((char *) sv + pitch * 47) + threadID) = 0.000000e+00; + *((real * )((char *) sv + pitch * 48) + threadID) = 0.000000e+00; + } + + if(use_adpt_dt) { + *((real *)((char *)sv + pitch * 49) + threadID) = min_dt; // dt + *((real *)((char *)sv + pitch * 50) + threadID) = 0.0; // time_new + *((real *)((char *)sv + pitch * 51) + threadID) = 0.0; // previous dt + } + } +} + +__global__ void kernel_set_model_initial_conditions_endo_mid_epi(real *sv, int num_volumes, size_t pitch, bool use_adpt_dt, real min_dt,\ + real *initial_endo, real *initial_epi, real *initial_mid, real *transmurality) { + int threadID = blockDim.x * blockIdx.x + threadIdx.x; + + if (threadID < num_volumes) { + for (int i = 0; i < NEQ; i++) { + if (transmurality[threadID] == ENDO) + *((real * )((char *) sv + pitch * i) + threadID) = initial_endo[i]; + else if (transmurality[threadID] == EPI) + *((real * )((char *) sv + pitch * i) + threadID) = initial_epi[i]; + else + *((real * )((char *) sv + pitch * i) + threadID) = initial_mid[i]; + } + + if(use_adpt_dt) { + *((real *)((char *)sv + pitch * (NEQ+1)) + threadID) = min_dt; // dt + *((real *)((char *)sv + pitch * (NEQ+2)) + threadID) = 0.0; // time_new + *((real *)((char *)sv + pitch * (NEQ+3)) + threadID) = 0.0; // previous dt + } + } +} + +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 timestep to solve the ODEs\n"); + } else { + log_info("Using Fixed timestep 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)); + + // Get initial condition from extra_data + real *initial_conditions_endo = NULL; + real *initial_conditions_epi = NULL; + real *initial_conditions_mid = NULL; + real *transmurality = NULL; + real *initial_conditions_endo_device = NULL; + real *initial_conditions_epi_device = NULL; + real *initial_conditions_mid_device = NULL; + real *transmurality_device = NULL; + + if(solver->ode_extra_data) { + struct extra_data_for_torord_land *extra_data = (struct extra_data_for_torord_land*)solver->ode_extra_data; + initial_conditions_endo = extra_data->initial_ss_endo; + initial_conditions_epi = extra_data->initial_ss_epi; + initial_conditions_mid = extra_data->initial_ss_mid; + transmurality = extra_data->transmurality; + check_cuda_error(cudaMalloc((void **)&initial_conditions_endo_device, sizeof(real)*NEQ)); + check_cuda_error(cudaMemcpy(initial_conditions_endo_device, initial_conditions_endo, sizeof(real)*NEQ, cudaMemcpyHostToDevice)); + check_cuda_error(cudaMalloc((void **)&initial_conditions_epi_device, sizeof(real)*NEQ)); + check_cuda_error(cudaMemcpy(initial_conditions_epi_device, initial_conditions_epi, sizeof(real)*NEQ, cudaMemcpyHostToDevice)); + check_cuda_error(cudaMalloc((void **)&initial_conditions_mid_device, sizeof(real)*NEQ)); + check_cuda_error(cudaMemcpy(initial_conditions_mid_device, initial_conditions_mid, sizeof(real)*NEQ, cudaMemcpyHostToDevice)); + check_cuda_error(cudaMalloc((void **)&transmurality_device, sizeof(real)*num_volumes)); + check_cuda_error(cudaMemcpy(transmurality_device, transmurality, sizeof(real)*num_volumes, cudaMemcpyHostToDevice)); + } + else { + log_info("[INFO] You should supply a mask function to tag the cells when using this mixed model!\n"); + log_info("[INFO] Considering all cells ENDO!\n"); + } + + if (solver->ode_extra_data) { + kernel_set_model_initial_conditions_endo_mid_epi<<>>(solver->sv, num_volumes, pitch_h, use_adpt_dt, solver->min_dt,\ + initial_conditions_endo_device, initial_conditions_epi_device, initial_conditions_mid_device,\ + transmurality_device); + } + else { + kernel_set_model_initial_conditions<<>>(solver->sv, num_volumes, pitch_h, use_adpt_dt, solver->min_dt); + } + + check_cuda_error(cudaPeekAtLastError()); + cudaDeviceSynchronize(); + + check_cuda_error(cudaFree(initial_conditions_endo_device)); + check_cuda_error(cudaFree(initial_conditions_epi_device)); + check_cuda_error(cudaFree(initial_conditions_mid_device)); + check_cuda_error(cudaFree(transmurality_device)); + + 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 = NULL; + 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 if exists + uint32_t num_volumes = ode_solver->original_num_cells; + real *transmurality = NULL; + real *transmurality_device = NULL; + int num_extra_parameters = 20; + real extra_par[num_extra_parameters]; + real *extra_par_device = NULL; + if(ode_solver->ode_extra_data) { + struct extra_data_for_torord_land *extra_data = (struct extra_data_for_torord_land*)ode_solver->ode_extra_data; + extra_par[0] = extra_data->INa_Multiplier; + extra_par[1] = extra_data->INaL_Multiplier; + extra_par[2] = extra_data->INaCa_Multiplier; + extra_par[3] = extra_data->INaK_Multiplier; + extra_par[4] = extra_data->INab_Multiplier; + extra_par[5] = extra_data->Ito_Multiplier; + extra_par[6] = extra_data->IKr_Multiplier; + extra_par[7] = extra_data->IKs_Multiplier; + extra_par[8] = extra_data->IK1_Multiplier; + extra_par[9] = extra_data->IKb_Multiplier; + extra_par[10] = extra_data->IKCa_Multiplier; + extra_par[11] = extra_data->ICaL_Multiplier; + extra_par[12] = extra_data->ICab_Multiplier; + extra_par[13] = extra_data->IpCa_Multiplier; + extra_par[14] = extra_data->ICaCl_Multiplier; + extra_par[15] = extra_data->IClb_Multiplier; + extra_par[16] = extra_data->Jrel_Multiplier; + extra_par[17] = extra_data->Jup_Multiplier; + extra_par[18] = extra_data->aCaMK_Multiplier; + extra_par[19] = extra_data->taurelp_Multiplier; + transmurality = extra_data->transmurality; + + check_cuda_error(cudaMalloc((void **)&transmurality_device, sizeof(real)*num_volumes)); + check_cuda_error(cudaMemcpy(transmurality_device, transmurality, sizeof(real)*num_volumes, cudaMemcpyHostToDevice)); + + check_cuda_error(cudaMalloc((void **)&extra_par_device, sizeof(real)*num_extra_parameters)); + check_cuda_error(cudaMemcpy(extra_par_device, extra_par, sizeof(real)*num_extra_parameters, cudaMemcpyHostToDevice)); + } + else { + // Default: initialize all current modifiers + for (uint32_t i = 0; i < num_extra_parameters; i++) { + if (i == 9) + extra_par[i] = 0.0; + else + extra_par[i] = 1.0; + } + + check_cuda_error(cudaMalloc((void **)&extra_par_device, sizeof(real)*num_extra_parameters)); + check_cuda_error(cudaMemcpy(extra_par_device, extra_par, sizeof(real)*num_extra_parameters, cudaMemcpyHostToDevice)); + } + + // Transmurality mapping is defined on 'extra_data' function + if (ode_solver->ode_extra_data) { + solve_endo_mid_epi_gpu<<>>(current_t, dt, sv, stims_currents_device, cells_to_solve_device, transmurality_device, extra_par_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); + } + // No transmurality: all cells ENDO + else { + solve_gpu<<>>(current_t, dt, sv, stims_currents_device, cells_to_solve_device, extra_par_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); + } + + check_cuda_error(cudaPeekAtLastError()); + + if (stims_currents_device) check_cuda_error(cudaFree(stims_currents_device)); + if (cells_to_solve_device) check_cuda_error(cudaFree(cells_to_solve_device)); + if (transmurality_device) check_cuda_error(cudaFree(transmurality_device)); + if (extra_par_device) check_cuda_error(cudaFree(extra_par_device)); +} + +__global__ void solve_gpu(real cur_time, real dt, real *sv, real *stim_currents, uint32_t *cells_to_solve, real *extra_params,\ + uint32_t num_cells_to_solve, int num_steps, size_t pitch, bool use_adpt, real abstol, real reltol, real max_dt) { + const real TOLERANCE = 1e-8; + int threadID = blockDim.x * blockIdx.x + threadIdx.x; + int sv_id; + + // 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]; + real a[NEQ], b[NEQ]; + + for(int n = 0; n < num_steps; ++n) { + + RHS_RL_gpu(a, b, sv, rDY, stim_currents[threadID], 0.0, extra_params, sv_id, dt, pitch, false); + + // Solve variables based on its type: + // Non-linear = Euler + // Hodkin-Huxley = Rush-Larsen || Euler (if 'a' coefficient is too small) + SOLVE_EQUATION_EULER_GPU(0); // v + SOLVE_EQUATION_EULER_GPU(1); // nai + SOLVE_EQUATION_EULER_GPU(2); // nass + SOLVE_EQUATION_EULER_GPU(3); // ki + SOLVE_EQUATION_EULER_GPU(4); // kss + SOLVE_EQUATION_EULER_GPU(5); // cai + SOLVE_EQUATION_EULER_GPU(6); // cass + SOLVE_EQUATION_EULER_GPU(7); // cansr + SOLVE_EQUATION_EULER_GPU(8); // cajsr + SOLVE_EQUATION_RUSH_LARSEN_GPU(9); // m + SOLVE_EQUATION_RUSH_LARSEN_GPU(10); // hp + SOLVE_EQUATION_RUSH_LARSEN_GPU(11); // h + SOLVE_EQUATION_RUSH_LARSEN_GPU(12); // j + SOLVE_EQUATION_RUSH_LARSEN_GPU(13); // jp + SOLVE_EQUATION_RUSH_LARSEN_GPU(14); // mL + SOLVE_EQUATION_RUSH_LARSEN_GPU(15); // hL + SOLVE_EQUATION_RUSH_LARSEN_GPU(16); // hLp + SOLVE_EQUATION_RUSH_LARSEN_GPU(17); // a + SOLVE_EQUATION_RUSH_LARSEN_GPU(18); // iF + SOLVE_EQUATION_RUSH_LARSEN_GPU(19); // iS + SOLVE_EQUATION_RUSH_LARSEN_GPU(20); // ap + SOLVE_EQUATION_RUSH_LARSEN_GPU(21); // iFp + SOLVE_EQUATION_RUSH_LARSEN_GPU(22); // iSp + SOLVE_EQUATION_RUSH_LARSEN_GPU(23); // d + SOLVE_EQUATION_RUSH_LARSEN_GPU(24); // ff + SOLVE_EQUATION_RUSH_LARSEN_GPU(25); // fs + SOLVE_EQUATION_RUSH_LARSEN_GPU(26); // fcaf + SOLVE_EQUATION_RUSH_LARSEN_GPU(27); // fcas + SOLVE_EQUATION_RUSH_LARSEN_GPU(28); // jca + SOLVE_EQUATION_EULER_GPU(29); // nca + SOLVE_EQUATION_EULER_GPU(30); // nca_i + SOLVE_EQUATION_RUSH_LARSEN_GPU(31); // ffp + SOLVE_EQUATION_RUSH_LARSEN_GPU(32); // fcafp + SOLVE_EQUATION_RUSH_LARSEN_GPU(33); // xs1 + SOLVE_EQUATION_RUSH_LARSEN_GPU(34); // xs2 + SOLVE_EQUATION_RUSH_LARSEN_GPU(35); // Jrel_np + SOLVE_EQUATION_EULER_GPU(36); // CaMKt + SOLVE_EQUATION_EULER_GPU(37); // ikr_c0 + SOLVE_EQUATION_EULER_GPU(38); // ikr_c1 + SOLVE_EQUATION_EULER_GPU(39); // ikr_c2 + SOLVE_EQUATION_EULER_GPU(40); // ikr_o + SOLVE_EQUATION_EULER_GPU(41); // ikr_i + SOLVE_EQUATION_RUSH_LARSEN_GPU(42); // Jrel_p + // --------------------------------------------------- + // Land-Niederer + SOLVE_EQUATION_EULER_GPU(43); // XS + SOLVE_EQUATION_EULER_GPU(44); // XW + SOLVE_EQUATION_EULER_GPU(45); // Ca_TRPN + SOLVE_EQUATION_EULER_GPU(46); // TmBlocked + SOLVE_EQUATION_EULER_GPU(47); // ZETAS + SOLVE_EQUATION_EULER_GPU(48); // ZETAW + } + } else { + solve_forward_euler_gpu_adpt(sv, stim_currents[threadID], 0.0, extra_params, cur_time + max_dt, sv_id, pitch, abstol, reltol, dt, max_dt); + } + } +} + +__global__ void solve_endo_mid_epi_gpu(real cur_time, real dt, real *sv, real *stim_currents, uint32_t *cells_to_solve, real *transmurality, real *extra_params,\ + uint32_t num_cells_to_solve, int num_steps, size_t pitch, bool use_adpt, real abstol, real reltol, real max_dt) { + const real TOLERANCE = 1e-8; + int threadID = blockDim.x * blockIdx.x + threadIdx.x; + int sv_id; + + // 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]; + real a[NEQ], b[NEQ]; + + for(int n = 0; n < num_steps; ++n) { + + RHS_RL_gpu(a, b, sv, rDY, stim_currents[threadID], transmurality[threadID], extra_params, sv_id, dt, pitch, false); + + // Solve variables based on its type: + // Non-linear = Euler + // Hodkin-Huxley = Rush-Larsen || Euler (if 'a' coefficient is too small) + SOLVE_EQUATION_EULER_GPU(0); // v + SOLVE_EQUATION_EULER_GPU(1); // nai + SOLVE_EQUATION_EULER_GPU(2); // nass + SOLVE_EQUATION_EULER_GPU(3); // ki + SOLVE_EQUATION_EULER_GPU(4); // kss + SOLVE_EQUATION_EULER_GPU(5); // cai + SOLVE_EQUATION_EULER_GPU(6); // cass + SOLVE_EQUATION_EULER_GPU(7); // cansr + SOLVE_EQUATION_EULER_GPU(8); // cajsr + SOLVE_EQUATION_RUSH_LARSEN_GPU(9); // m + SOLVE_EQUATION_RUSH_LARSEN_GPU(10); // hp + SOLVE_EQUATION_RUSH_LARSEN_GPU(11); // h + SOLVE_EQUATION_RUSH_LARSEN_GPU(12); // j + SOLVE_EQUATION_RUSH_LARSEN_GPU(13); // jp + SOLVE_EQUATION_RUSH_LARSEN_GPU(14); // mL + SOLVE_EQUATION_RUSH_LARSEN_GPU(15); // hL + SOLVE_EQUATION_RUSH_LARSEN_GPU(16); // hLp + SOLVE_EQUATION_RUSH_LARSEN_GPU(17); // a + SOLVE_EQUATION_RUSH_LARSEN_GPU(18); // iF + SOLVE_EQUATION_RUSH_LARSEN_GPU(19); // iS + SOLVE_EQUATION_RUSH_LARSEN_GPU(20); // ap + SOLVE_EQUATION_RUSH_LARSEN_GPU(21); // iFp + SOLVE_EQUATION_RUSH_LARSEN_GPU(22); // iSp + SOLVE_EQUATION_RUSH_LARSEN_GPU(23); // d + SOLVE_EQUATION_RUSH_LARSEN_GPU(24); // ff + SOLVE_EQUATION_RUSH_LARSEN_GPU(25); // fs + SOLVE_EQUATION_RUSH_LARSEN_GPU(26); // fcaf + SOLVE_EQUATION_RUSH_LARSEN_GPU(27); // fcas + SOLVE_EQUATION_RUSH_LARSEN_GPU(28); // jca + SOLVE_EQUATION_EULER_GPU(29); // nca + SOLVE_EQUATION_EULER_GPU(30); // nca_i + SOLVE_EQUATION_RUSH_LARSEN_GPU(31); // ffp + SOLVE_EQUATION_RUSH_LARSEN_GPU(32); // fcafp + SOLVE_EQUATION_RUSH_LARSEN_GPU(33); // xs1 + SOLVE_EQUATION_RUSH_LARSEN_GPU(34); // xs2 + SOLVE_EQUATION_RUSH_LARSEN_GPU(35); // Jrel_np + SOLVE_EQUATION_EULER_GPU(36); // CaMKt + SOLVE_EQUATION_EULER_GPU(37); // ikr_c0 + SOLVE_EQUATION_EULER_GPU(38); // ikr_c1 + SOLVE_EQUATION_EULER_GPU(39); // ikr_c2 + SOLVE_EQUATION_EULER_GPU(40); // ikr_o + SOLVE_EQUATION_EULER_GPU(41); // ikr_i + SOLVE_EQUATION_RUSH_LARSEN_GPU(42); // Jrel_p + // --------------------------------------------------- + // Land-Niederer + SOLVE_EQUATION_EULER_GPU(43); // XS + SOLVE_EQUATION_EULER_GPU(44); // XW + SOLVE_EQUATION_EULER_GPU(45); // Ca_TRPN + SOLVE_EQUATION_EULER_GPU(46); // TmBlocked + SOLVE_EQUATION_EULER_GPU(47); // ZETAS + SOLVE_EQUATION_EULER_GPU(48); // ZETAW + } + } else { + solve_forward_euler_gpu_adpt(sv, stim_currents[threadID], transmurality[threadID], extra_params, cur_time + max_dt, sv_id, pitch, abstol, reltol, dt, max_dt); + } + } +} + +inline __device__ void solve_forward_euler_gpu_adpt(real *sv, real stim_curr, real mapping, real *extra_params, real final_time, int thread_id, size_t pitch, real abstol, real reltol, real min_dt, real max_dt) { + + #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, mapping, extra_params, thread_id, dt, pitch, true); + 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, mapping, extra_params, thread_id, dt, pitch, true); + 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; +} + +inline __device__ void RHS_gpu(real *sv, real *rDY_, real stim_current, \ + real mapping, real *extra_params, int threadID_, real dt, size_t pitch, bool use_adpt_dt) { + + // Current modifiers + real INa_Multiplier = extra_params[0]; + real INaL_Multiplier = extra_params[1]; + real INaCa_Multiplier = extra_params[2]; + real INaK_Multiplier = extra_params[3]; + real INab_Multiplier = extra_params[4]; + real Ito_Multiplier = extra_params[5]; + real IKr_Multiplier = extra_params[6]; + real IKs_Multiplier = extra_params[7]; + real IK1_Multiplier = extra_params[8]; + real IKb_Multiplier = extra_params[9]; + real IKCa_Multiplier = extra_params[10]; + real ICaL_Multiplier = extra_params[11]; + real ICab_Multiplier = extra_params[12]; + real IpCa_Multiplier = extra_params[13]; + real ICaCl_Multiplier = extra_params[14]; + real IClb_Multiplier = extra_params[15]; + real Jrel_Multiplier = extra_params[16]; + real Jup_Multiplier = extra_params[17]; + real aCaMK_Multiplier = extra_params[18]; + real taurelp_Multiplier = extra_params[19]; + + // Get the celltype for the current cell + real celltype = mapping; + + // Get the stimulus current from the current cell + real calc_I_stim = stim_current; + + // Constant variables + const real cli = 24; // Intracellular Cl [mM] + const real clo = 150; // Extracellular Cl [mM] + + // State variables + real v; + real nai; + real nass; + real ki; + real kss; + real cai; + real cass; + real cansr; + real cajsr; + real m; + real hp; + real h; + real j; + real jp; + real mL; + real hL; + real hLp; + real a; + real iF; + real iS; + real ap; + real iFp; + real iSp; + + // ical + real d; + real ff; + real fs; + real fcaf; + real fcas; + real jca; + real nca; + real nca_i; + real ffp; + real fcafp; + + real xs1; + real xs2; + real Jrel_np; + real CaMKt; + + // new MM ICaL states + real ikr_c0; + real ikr_c1; + real ikr_c2; + real ikr_o; + real ikr_i; + real Jrel_p; + + real XS; + real XW; + real Ca_TRPN; + real TmBlocked; + real ZETAS; + real ZETAW; + + if (use_adpt_dt) { + v = sv[0]; + nai = sv[1]; + nass = sv[2]; + ki = sv[3]; + kss = sv[4]; + cai = sv[5]; + cass = sv[6]; + cansr = sv[7]; + cajsr = sv[8]; + m = sv[9]; + hp = sv[10]; + h = sv[11]; + j = sv[12]; + jp = sv[13]; + mL = sv[14]; + hL = sv[15]; + hLp = sv[16]; + a = sv[17]; + iF = sv[18]; + iS = sv[19]; + ap = sv[20]; + iFp = sv[21]; + iSp = sv[22]; + + // ical + d = sv[23]; + ff = sv[24]; + fs = sv[25]; + fcaf = sv[26]; + fcas = sv[27]; + jca = sv[28]; + nca = sv[29]; + nca_i = sv[30]; + ffp = sv[31]; + fcafp = sv[32]; + + xs1 = sv[33]; + xs2 = sv[34]; + Jrel_np = sv[35]; + CaMKt = sv[36]; + + // new MM ICaL states + ikr_c0 = sv[37]; + ikr_c1 = sv[38]; + ikr_c2 = sv[39]; + ikr_o = sv[40]; + ikr_i = sv[41]; + Jrel_p = sv[42]; + + XS = (sv[43] > 0.0) ? sv[43] : 0.0; + XW = (sv[44] > 0.0) ? sv[44] : 0.0; + Ca_TRPN = (sv[45] > 0.0) ? sv[45] : 0.0; + TmBlocked = sv[46]; + ZETAS = sv[47]; + ZETAW = sv[48]; + } else { + v = *((real *)((char *)sv + pitch * 0) + threadID_); + nai = *((real *)((char *)sv + pitch * 1) + threadID_); + nass = *((real *)((char *)sv + pitch * 2) + threadID_); + ki = *((real *)((char *)sv + pitch * 3) + threadID_); + kss = *((real *)((char *)sv + pitch * 4) + threadID_); + cai = *((real *)((char *)sv + pitch * 5) + threadID_); + cass = *((real *)((char *)sv + pitch * 6) + threadID_); + cansr = *((real *)((char *)sv + pitch * 7) + threadID_); + cajsr = *((real *)((char *)sv + pitch * 8) + threadID_); + m = *((real *)((char *)sv + pitch * 9) + threadID_); + hp = *((real *)((char *)sv + pitch * 10) + threadID_); + h = *((real *)((char *)sv + pitch * 11) + threadID_); + j = *((real *)((char *)sv + pitch * 12) + threadID_); + jp = *((real *)((char *)sv + pitch * 13) + threadID_); + mL = *((real *)((char *)sv + pitch * 14) + threadID_); + hL = *((real *)((char *)sv + pitch * 15) + threadID_); + hLp = *((real *)((char *)sv + pitch * 16) + threadID_); + a = *((real *)((char *)sv + pitch * 17) + threadID_); + iF = *((real *)((char *)sv + pitch * 18) + threadID_); + iS = *((real *)((char *)sv + pitch * 19) + threadID_); + ap = *((real *)((char *)sv + pitch * 20) + threadID_); + iFp = *((real *)((char *)sv + pitch * 21) + threadID_); + iSp = *((real *)((char *)sv + pitch * 22) + threadID_); + + // ical + d = *((real *)((char *)sv + pitch * 23) + threadID_); + ff = *((real *)((char *)sv + pitch * 24) + threadID_); + fs = *((real *)((char *)sv + pitch * 25) + threadID_); + fcaf = *((real *)((char *)sv + pitch * 26) + threadID_); + fcas = *((real *)((char *)sv + pitch * 27) + threadID_); + jca = *((real *)((char *)sv + pitch * 28) + threadID_); + nca = *((real *)((char *)sv + pitch * 29) + threadID_); + nca_i = *((real *)((char *)sv + pitch * 30) + threadID_); + ffp = *((real *)((char *)sv + pitch * 31) + threadID_); + fcafp = *((real *)((char *)sv + pitch * 32) + threadID_); + + xs1 = *((real *)((char *)sv + pitch * 33) + threadID_); + xs2 = *((real *)((char *)sv + pitch * 34) + threadID_); + Jrel_np = *((real *)((char *)sv + pitch * 35) + threadID_); + CaMKt = *((real *)((char *)sv + pitch * 36) + threadID_); + + // new MM ICaL states + ikr_c0 = *((real *)((char *)sv + pitch * 37) + threadID_); + ikr_c1 = *((real *)((char *)sv + pitch * 38) + threadID_); + ikr_c2 = *((real *)((char *)sv + pitch * 39) + threadID_); + ikr_o = *((real *)((char *)sv + pitch * 40) + threadID_); + ikr_i = *((real *)((char *)sv + pitch * 41) + threadID_); + Jrel_p = *((real *)((char *)sv + pitch * 42) + threadID_); + + // Land-Niederer model + XS = fmaxf(0,(*((real *)((char *)sv + pitch * 43) + threadID_))); + XW = fmaxf(0,(*((real *)((char *)sv + pitch * 44) + threadID_))); + Ca_TRPN = fmaxf(0,(*((real *)((char *)sv + pitch * 45) + threadID_))); + TmBlocked = (*((real *)((char *)sv + pitch * 46) + threadID_)); + ZETAS = (*((real *)((char *)sv + pitch * 47) + threadID_)); + ZETAW = (*((real *)((char *)sv + pitch * 48) + threadID_)); + } + + #include "ToRORd_Land_mixed_endo_mid_epi.common.c" +} + +inline __device__ void RHS_RL_gpu(real *a_, real *b_, real *sv, real *rDY_, real stim_current, \ + real mapping, real *extra_params, int threadID_, real dt, size_t pitch, bool use_adpt_dt) { + + // Current modifiers + real INa_Multiplier = extra_params[0]; + real INaL_Multiplier = extra_params[1]; + real INaCa_Multiplier = extra_params[2]; + real INaK_Multiplier = extra_params[3]; + real INab_Multiplier = extra_params[4]; + real Ito_Multiplier = extra_params[5]; + real IKr_Multiplier = extra_params[6]; + real IKs_Multiplier = extra_params[7]; + real IK1_Multiplier = extra_params[8]; + real IKb_Multiplier = extra_params[9]; + real IKCa_Multiplier = extra_params[10]; + real ICaL_Multiplier = extra_params[11]; + real ICab_Multiplier = extra_params[12]; + real IpCa_Multiplier = extra_params[13]; + real ICaCl_Multiplier = extra_params[14]; + real IClb_Multiplier = extra_params[15]; + real Jrel_Multiplier = extra_params[16]; + real Jup_Multiplier = extra_params[17]; + real aCaMK_Multiplier = extra_params[18]; + real taurelp_Multiplier = extra_params[19]; + + // Get the celltype for the current cell + real celltype = mapping; + + // Get the stimulus current from the current cell + real calc_I_stim = stim_current; + + // Constant variables + const real cli = 24; // Intracellular Cl [mM] + const real clo = 150; // Extracellular Cl [mM] + + // State variables + real v; + real nai; + real nass; + real ki; + real kss; + real cai; + real cass; + real cansr; + real cajsr; + real m; + real hp; + real h; + real j; + real jp; + real mL; + real hL; + real hLp; + real a; + real iF; + real iS; + real ap; + real iFp; + real iSp; + + // ical + real d; + real ff; + real fs; + real fcaf; + real fcas; + real jca; + real nca; + real nca_i; + real ffp; + real fcafp; + + real xs1; + real xs2; + real Jrel_np; + real CaMKt; + + // new MM ICaL states + real ikr_c0; + real ikr_c1; + real ikr_c2; + real ikr_o; + real ikr_i; + real Jrel_p; + + real XS; + real XW; + real Ca_TRPN; + real TmBlocked; + real ZETAS; + real ZETAW; + + if (use_adpt_dt) { + v = sv[0]; + nai = sv[1]; + nass = sv[2]; + ki = sv[3]; + kss = sv[4]; + cai = sv[5]; + cass = sv[6]; + cansr = sv[7]; + cajsr = sv[8]; + m = sv[9]; + hp = sv[10]; + h = sv[11]; + j = sv[12]; + jp = sv[13]; + mL = sv[14]; + hL = sv[15]; + hLp = sv[16]; + a = sv[17]; + iF = sv[18]; + iS = sv[19]; + ap = sv[20]; + iFp = sv[21]; + iSp = sv[22]; + + // ical + d = sv[23]; + ff = sv[24]; + fs = sv[25]; + fcaf = sv[26]; + fcas = sv[27]; + jca = sv[28]; + nca = sv[29]; + nca_i = sv[30]; + ffp = sv[31]; + fcafp = sv[32]; + + xs1 = sv[33]; + xs2 = sv[34]; + Jrel_np = sv[35]; + CaMKt = sv[36]; + + // new MM ICaL states + ikr_c0 = sv[37]; + ikr_c1 = sv[38]; + ikr_c2 = sv[39]; + ikr_o = sv[40]; + ikr_i = sv[41]; + Jrel_p = sv[42]; + + XS = (sv[43] > 0.0) ? sv[43] : 0.0; + XW = (sv[44] > 0.0) ? sv[44] : 0.0; + Ca_TRPN = (sv[45] > 0.0) ? sv[45] : 0.0; + TmBlocked = sv[46]; + ZETAS = sv[47]; + ZETAW = sv[48]; + } else { + v = *((real *)((char *)sv + pitch * 0) + threadID_); + nai = *((real *)((char *)sv + pitch * 1) + threadID_); + nass = *((real *)((char *)sv + pitch * 2) + threadID_); + ki = *((real *)((char *)sv + pitch * 3) + threadID_); + kss = *((real *)((char *)sv + pitch * 4) + threadID_); + cai = *((real *)((char *)sv + pitch * 5) + threadID_); + cass = *((real *)((char *)sv + pitch * 6) + threadID_); + cansr = *((real *)((char *)sv + pitch * 7) + threadID_); + cajsr = *((real *)((char *)sv + pitch * 8) + threadID_); + m = *((real *)((char *)sv + pitch * 9) + threadID_); + hp = *((real *)((char *)sv + pitch * 10) + threadID_); + h = *((real *)((char *)sv + pitch * 11) + threadID_); + j = *((real *)((char *)sv + pitch * 12) + threadID_); + jp = *((real *)((char *)sv + pitch * 13) + threadID_); + mL = *((real *)((char *)sv + pitch * 14) + threadID_); + hL = *((real *)((char *)sv + pitch * 15) + threadID_); + hLp = *((real *)((char *)sv + pitch * 16) + threadID_); + a = *((real *)((char *)sv + pitch * 17) + threadID_); + iF = *((real *)((char *)sv + pitch * 18) + threadID_); + iS = *((real *)((char *)sv + pitch * 19) + threadID_); + ap = *((real *)((char *)sv + pitch * 20) + threadID_); + iFp = *((real *)((char *)sv + pitch * 21) + threadID_); + iSp = *((real *)((char *)sv + pitch * 22) + threadID_); + + // ical + d = *((real *)((char *)sv + pitch * 23) + threadID_); + ff = *((real *)((char *)sv + pitch * 24) + threadID_); + fs = *((real *)((char *)sv + pitch * 25) + threadID_); + fcaf = *((real *)((char *)sv + pitch * 26) + threadID_); + fcas = *((real *)((char *)sv + pitch * 27) + threadID_); + jca = *((real *)((char *)sv + pitch * 28) + threadID_); + nca = *((real *)((char *)sv + pitch * 29) + threadID_); + nca_i = *((real *)((char *)sv + pitch * 30) + threadID_); + ffp = *((real *)((char *)sv + pitch * 31) + threadID_); + fcafp = *((real *)((char *)sv + pitch * 32) + threadID_); + + xs1 = *((real *)((char *)sv + pitch * 33) + threadID_); + xs2 = *((real *)((char *)sv + pitch * 34) + threadID_); + Jrel_np = *((real *)((char *)sv + pitch * 35) + threadID_); + CaMKt = *((real *)((char *)sv + pitch * 36) + threadID_); + + // new MM ICaL states + ikr_c0 = *((real *)((char *)sv + pitch * 37) + threadID_); + ikr_c1 = *((real *)((char *)sv + pitch * 38) + threadID_); + ikr_c2 = *((real *)((char *)sv + pitch * 39) + threadID_); + ikr_o = *((real *)((char *)sv + pitch * 40) + threadID_); + ikr_i = *((real *)((char *)sv + pitch * 41) + threadID_); + Jrel_p = *((real *)((char *)sv + pitch * 42) + threadID_); + + // Land-Niederer model + XS = fmaxf(0,(*((real *)((char *)sv + pitch * 43) + threadID_))); + XW = fmaxf(0,(*((real *)((char *)sv + pitch * 44) + threadID_))); + Ca_TRPN = fmaxf(0,(*((real *)((char *)sv + pitch * 45) + threadID_))); + TmBlocked = (*((real *)((char *)sv + pitch * 46) + threadID_)); + ZETAS = (*((real *)((char *)sv + pitch * 47) + threadID_)); + ZETAW = (*((real *)((char *)sv + pitch * 48) + threadID_)); + } + + #include "ToRORd_Land_mixed_endo_mid_epi_RL.common.c" +} diff --git a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.h b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.h new file mode 100644 index 00000000..a2a874e3 --- /dev/null +++ b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.h @@ -0,0 +1,107 @@ +#ifndef MONOALG3D_MODEL_TORORD_LAND_MIXED_ENDO_MID_EPI_H +#define MONOALG3D_MODEL_TORORD_LAND_MIXED_ENDO_MID_EPI_H + +// (Margara, F., Wang, Z.J., Levrero-Florencio, F., Santiago, A., Vázquez, M., Bueno-Orovio, A., +// and Rodriguez, B. (2020). In-silico human electro-mechanical ventricular modelling and simulation for +// drug-induced pro-arrhythmia and inotropic risk assessment. Progress in Biophysics and Molecular Biology). + + +#include "../model_common.h" +#include "../../extra_data_library/helper_functions.h" +#include + +#define NEQ 49 +#define INITIAL_V (-8.863699e+01) + +#define ENDO 0.0 +#define MID 1.0 +#define EPI 2.0 + +// CPU macros +#define SOLVE_EQUATION_CONSTANT_CPU(id) sv[id] = rY[id] + +#define SOLVE_EQUATION_EULER_CPU(id) sv[id] = dt * rDY[id] + rY[id] + +#define SOLVE_EQUATION_RUSH_LARSEN_CPU(id) sv[id] = (abs(a[id]) < TOLERANCE) ? rY[id] + dt * (rY[id] * a[id] + b[id]) : \ + exp(a[id] * dt)*(rY[id] + (b[id] / a[id])) - (b[id] / a[id] ) + +#define SOLVE_EQUATION_ADAPT_RUSH_LARSEN_CONSTANT_CPU(id) sv[id] = rDY[id]; + +#define SOLVE_EQUATION_ADAPT_RUSH_LARSEN_EULER_CPU(id) edos_old_aux_[id] = sv[id]; \ + edos_new_euler_[id] = _k1__[id] * *dt + edos_old_aux_[id]; \ + sv[id] = edos_new_euler_[id] + +#define SOLVE_EQUATION_ADAPT_RUSH_LARSEN_RL_CPU(id) edos_old_aux_[id] = sv[id]; \ + edos_new_euler_[id] = (fabs(a_[id]) < abs_tol) ? edos_old_aux_[id] + (edos_old_aux_[id] * a_[id] + b_[id])*(*dt) : exp(a_[id]*(*dt))*(edos_old_aux_[id] + (b_[id] / a_[id])) - (b_[id] / a_[id]); \ + sv[id] = edos_new_euler_[id] + +#define SOLVE_ERROR_ADAPT_RUSH_LARSEN_EULER_CPU(id) _k2__[id] = rDY[id]; \ + f = (_k1__[id] + _k2__[id]) * 0.5; \ + y_2nd_order = edos_old_aux_[id] + (*dt) * f; \ + auxError = (fabs(y_2nd_order) < abs_tol) ? fabs(edos_new_euler_[id] - abs_tol) : fabs( (y_2nd_order - edos_new_euler_[id])/(y_2nd_order) ); \ + greatestError = (auxError > greatestError) ? auxError : greatestError + +#define SOLVE_ERROR_ADAPT_RUSH_LARSEN_RL_CPU(id) _k2__[id] = rDY[id]; \ + as = (a_[id] + a_new[id]) * 0.5; \ + bs = (b_[id] + b_new[id]) * 0.5; \ + y_2nd_order = (fabs(as) < abs_tol) ? edos_old_aux_[id] + (*dt) * (edos_old_aux_[id]*as + bs) : exp(as*(*dt))*(edos_old_aux_[id] + (bs/as)) - (bs/as); \ + auxError = (fabs(y_2nd_order) < abs_tol) ? fabs(edos_new_euler_[id] - abs_tol) : fabs( (y_2nd_order - edos_new_euler_[id])/(y_2nd_order) ); \ + greatestError = (auxError > greatestError) ? auxError : greatestError + +// GPU macros +#define SOLVE_EQUATION_CONSTANT_GPU(id) *((real *)((char *)sv + pitch * id) + sv_id) = *((real *)((char *)sv + pitch * id) + sv_id) + +#define SOLVE_EQUATION_EULER_GPU(id) *((real *)((char *)sv + pitch * id) + sv_id) = *((real *)((char *)sv + pitch * id) + sv_id) + dt * rDY[id] + +#define SOLVE_EQUATION_RUSH_LARSEN_GPU(id) *((real *)((char *)sv + pitch * id) + sv_id) = (abs(a[id]) < TOLERANCE) ? *((real *)((char *)sv + pitch * id) + sv_id) + dt * ( *((real *)((char *)sv + pitch * id) + sv_id) * a[id] + b[id]) : \ + exp(a[id] * dt)*(*((real *)((char *)sv + pitch * id) + sv_id) + (b[id] / a[id])) - (b[id] / a[id]) + +#define SOLVE_EQUATION_ADAPT_RUSH_LARSEN_CONSTANT_GPU(id) sv_local[id] = rDY[id]; + +#define SOLVE_EQUATION_ADAPT_RUSH_LARSEN_EULER_GPU(id) edos_old_aux_[id] = sv_local[id]; \ + edos_new_euler_[id] = _k1__[id] * dt + edos_old_aux_[id]; \ + sv_local[id] = edos_new_euler_[id] + +#define SOLVE_EQUATION_ADAPT_RUSH_LARSEN_RL_GPU(id) edos_old_aux_[id] = sv_local[id]; \ + edos_new_euler_[id] = (a_[id] < abstol) ? edos_old_aux_[id] + (edos_old_aux_[id] * a_[id] + b_[id])*(dt) : exp(a_[id]*(dt))*(edos_old_aux_[id] + (b_[id] / a_[id])) - (b_[id] / a_[id]); \ + sv_local[id] = edos_new_euler_[id] + +#define SOLVE_ERROR_ADAPT_RUSH_LARSEN_EULER_GPU(id) _k2__[id] = rDY[id]; \ + f = (_k1__[id] + _k2__[id]) * 0.5; \ + y_2nd_order = edos_old_aux_[id] + (dt) * f; \ + auxError = (fabs(y_2nd_order) < abstol) ? fabs(edos_new_euler_[id] - abstol) : fabs( (y_2nd_order - edos_new_euler_[id])/(y_2nd_order) ); \ + greatestError = (auxError > greatestError) ? auxError : greatestError + +#define SOLVE_ERROR_ADAPT_RUSH_LARSEN_RL_GPU(id) _k2__[id] = rDY[id]; \ + as = (a_[id] + a_new[id]) * 0.5; \ + bs = (b_[id] + b_new[id]) * 0.5; \ + y_2nd_order = (fabs(as) < abstol) ? edos_old_aux_[id] + (dt) * (edos_old_aux_[id]*as + bs) : exp(as*(dt))*(edos_old_aux_[id] + (bs/as)) - (bs/as); \ + auxError = (fabs(y_2nd_order) < abstol) ? fabs(edos_new_euler_[id] - abstol) : fabs( (y_2nd_order - edos_new_euler_[id])/(y_2nd_order) ); \ + greatestError = (auxError > greatestError) ? auxError : greatestError + +#ifdef __CUDACC__ + +#include "../../gpu_utils/gpu_utils.h" + +__global__ void kernel_set_model_initial_conditions(real *sv, int num_volumes, size_t pitch, bool use_adpt_dt, real min_dt); +__global__ void kernel_set_model_initial_conditions_endo_mid_epi(real *sv, int num_volumes, size_t pitch, bool use_adpt_dt, real min_dt,\ + real *initial_endo, real *initial_epi, real *initial_mid, real *transmurality); + +__global__ void solve_gpu(real cur_time, real dt, real *sv, real *stim_currents, uint32_t *cells_to_solve, real *extra_params,\ + uint32_t num_cells_to_solve, int num_steps, size_t pitch, bool use_adpt, real abstol, real reltol, real max_dt); +__global__ void solve_endo_mid_epi_gpu(real cur_time, real dt, real *sv, real *stim_currents, uint32_t *cells_to_solve, real *transmurality, real *extra_params,\ + uint32_t num_cells_to_solve, int num_steps, size_t pitch, bool use_adpt, real abstol, real reltol, real max_dt); + +inline __device__ void RHS_gpu(real *sv, real *rDY_, real stim_current, real transmurality, real *extra_params, int threadID_, real dt, size_t pitch, bool use_adpt_dt); +inline __device__ void RHS_RL_gpu(real *a_, real *b_, real *sv, real *rDY_, real stim_current, real transmurality, real *extra_params, int threadID_, real dt, size_t pitch, bool use_adpt_dt); +inline __device__ void solve_forward_euler_gpu_adpt(real *sv, real stim_curr, real transmurality, real *extra_params, real final_time, int thread_id, size_t pitch, real abstol, real reltol, real min_dt, real max_dt); + +#endif + +void RHS_cpu(const real *sv, real *rDY_, real stim_current, real dt, real transmurality, real const *extra_params); +void RHS_RL_cpu (real *a_, real *b_, const real *sv, real *rDY_, real stim_current, real dt, real transmurality, real const *extra_params); +void solve_forward_euler_cpu_adpt(real *sv, real stim_curr, real transmurality, real final_time, int sv_id, struct ode_solver *solver, real const *extra_params); +void solve_model_ode_cpu(real dt, real *sv, real stim_current, real transmurality, real const *extra_params); + +#endif //MONOALG3D_MODEL_TORORD_LAND_MIXED_ENDO_MID_EPI_H + diff --git a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi_RL.common.c b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi_RL.common.c new file mode 100644 index 00000000..66cb582a --- /dev/null +++ b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi_RL.common.c @@ -0,0 +1,886 @@ +// Avoid NaN errors +//if (v >= 100.0) v = 100.0; +//if (v <= -100.0) v = -100.0; + +//real aux, aux2, aux3; + +// Changeable parameters +real nao = 140.0; +real cao = 1.8; +real ko = 5.0; +// Localization of ICaL and NCX: the fraction in junctional subspace +real ICaL_fractionSS = 0.8; +real INaCa_fractionSS = 0.35; +// Additional parameters +real ca50 = 0.805; +real trpnmax = 0.07; + +// INPUT CODE: +int mode = 0; // 0 = "intact", 1 = "skinned" +real lambda = 1.0; +real lambda_rate = 0.0; + +// EC parameters +real perm50 = 0.35; +real TRPN_n = 2; +real koff = 0.1; +real dr = 0.25; +real wfrac = 0.5; +real TOT_A = 25; +real ktm_unblock = 0.021; +real beta_1 = -2.4; +real beta_0 = 2.3; +real gamma = 0.0085; +real gamma_wu = 0.615; +real phi = 2.23; + +real nperm; +real Tref; +real nu; +real mu; +if (mode == 1) { + nperm = 2.2; + //ca50=2.5; + Tref = 40.5; + nu = 1; + mu = 1; +} +else { + nperm = 2.036; + //ca50=0.805; + //ca50 = 0.5; + //ca50 = 0.7; + Tref = 120; + nu = 7; + mu = 3; +} + +real k_ws = 0.004; +real k_uw = 0.026; + +real lambda_min = 0.87; +real lambda_max = 1.2; + +//real dydt[6] = {0,0,0,0,0,0}; + +k_ws = k_ws*mu; +k_uw = k_uw*nu; + +real cdw = phi*k_uw*(1-dr)*(1-wfrac)/((1-dr)*wfrac); +real cds = phi*k_ws*(1-dr)*wfrac/dr; +real k_wu = k_uw*(1/wfrac-1)-k_ws; +real k_su = k_ws*(1/dr-1)*wfrac; +real A = (0.25*TOT_A)/((1-dr)*wfrac+dr)*(dr/0.25); + +//real lambda0 = (lambda_max < lambda) ? lambda_max : lambda; +real lambda0 = fminf(lambda_max, lambda); +//real lambda_aux = (lambda_min < lambda0) ? lambda_min : lambda0; +//aux = 1+beta_0*(lambda0+lambda_aux-(1+lambda_min)); +//real Lfac = (aux > 0) ? aux : 0; +real Lfac = fmaxf(0,1+beta_0*(lambda0+fminf(lambda_min,lambda0)-(1+lambda_min))); + +real XU = (1-TmBlocked)-XW-XS; // unattached available xb = all - tm blocked - already prepowerstroke - already post-poststroke - no overlap +real xb_ws = k_ws*XW; +real xb_uw = k_uw*XU; +real xb_wu = k_wu*XW; +real xb_su = k_su*XS; + +//aux = ((ZETAS>0)*ZETAS > (ZETAS<-1)*(-ZETAS-1)) ? (ZETAS>0)*ZETAS : (ZETAS<-1)*(-ZETAS-1); +real gamma_rate=gamma*fmaxf((ZETAS>0)*ZETAS,(ZETAS<-1)*(-ZETAS-1)); +real xb_su_gamma=gamma_rate*XS; +real gamma_rate_w=gamma_wu*abs(ZETAW); // weak xbs don't like being strained +real xb_wu_gamma=gamma_rate_w*XW; + +//dydt[0] = xb_ws-xb_su-xb_su_gamma; +//dydt[1] = xb_uw-xb_wu-xb_ws-xb_wu_gamma; +real dXS = xb_ws-xb_su-xb_su_gamma; +real dXW = xb_uw-xb_wu-xb_ws-xb_wu_gamma; + +//aux = (lambda-1 < 0.2) ? (lambda-1) : 0.2; +//ca50 = ca50+beta_1*aux; +ca50=ca50+beta_1*fminf(0.2,lambda-1); +//dydt[2] = koff*(pow(((cai*1000)/ca50),TRPN_n)*(1-Ca_TRPN)-Ca_TRPN); // untouched +real dCa_TRPN = koff*(pow(((cai*1000)/ca50),TRPN_n)*(1-Ca_TRPN)-Ca_TRPN); + +real XSSS = dr*0.5; +real XWSS = (1-dr)*wfrac*0.5; +real ktm_block = ktm_unblock*(pow(perm50,nperm))*0.5/(0.5-XSSS-XWSS); + +//aux = pow(Ca_TRPN,-(nperm/2)); +//aux2 = pow(Ca_TRPN,(nperm/2)); +//aux3 = (100 < aux) ? 100 : aux; +//dydt[3] = ktm_block*aux3*XU-ktm_unblock*aux2*TmBlocked; +//real dTmBlocked = ktm_block*aux3*XU-ktm_unblock*aux2*TmBlocked; +real dTmBlocked = ktm_block*fminf(100, pow(Ca_TRPN,-(nperm/2)))*XU - ktm_unblock*( pow(Ca_TRPN,(nperm/2)))*TmBlocked; + +// velocity dependence -- assumes distortion resets on W->S +//dydt[4] = A*lambda_rate-cds*ZETAS; // - gamma_rate * ZETAS; +//dydt[5] = A*lambda_rate-cdw*ZETAW; // - gamma_rate_w * ZETAW; +real dZETAS = A*lambda_rate-cds*ZETAS; // - gamma_rate * ZETAS; +real dZETAW = A*lambda_rate-cdw*ZETAW; // - gamma_rate_w * ZETAW; + +// Active Force +real Ta = Lfac*(Tref/dr)*((ZETAS+1)*XS+(ZETAW)*XW); + +// physical constants +real R=8314.0; +real T=310.0; +real F=96485.0; + +// cell geometry +real L=0.01; +real rad=0.0011; +real vcell=1000*3.14*rad*rad*L; +real Ageo=2*3.14*rad*rad+2*3.14*rad*L; +real Acap=2*Ageo; +real vmyo=0.68*vcell; +real vnsr=0.0552*vcell; +real vjsr=0.0048*vcell; +real vss=0.02*vcell; + +real fkatp = 0.0; +real gkatp = 4.3195; + +// CaMK constants +real KmCaMK=0.15; +real aCaMK=0.05*aCaMK_Multiplier; +real bCaMK=0.00068; +real CaMKo=0.05; +real KmCaM=0.0015; +// update CaMK +real CaMKb=CaMKo*(1.0-CaMKt)/(1.0+KmCaM/cass); +real CaMKa=CaMKb+CaMKt; +real dCaMKt=aCaMK*CaMKb*(CaMKb+CaMKt)-bCaMK*CaMKt; // Euler + +// reversal potentials +real ENa=(R*T/F)*log(nao/nai); +real EK=(R*T/F)*log(ko/ki); +real PKNa=0.01833; +real EKs=(R*T/F)*log((ko+PKNa*nao)/(ki+PKNa*nai)); + +// convenient shorthand calculations +real vffrt=v*F*F/(R*T); +real vfrt=v*F/(R*T); +real frt = F/(R*T); + +real K_o_n = 5.0; +real A_atp = 2.0; +real K_atp = 0.25; +real akik = pow((ko / K_o_n), 0.24); +real bkik = (1.0 / (1.0 + pow((A_atp / K_atp), 2.0))); + +real fINap=(1.0/(1.0+KmCaMK/CaMKa)); +real fINaLp=(1.0/(1.0+KmCaMK/CaMKa)); +real fItop=(1.0/(1.0+KmCaMK/CaMKa)); +real fICaLp=(1.0/(1.0+KmCaMK/CaMKa)); + +// INa formulations +// The Grandi implementation updated with INa phosphorylation. +// m gate +real mss = 1 / (pow(1 + exp( -(56.86 + v) / 9.03 ),2)); +real taum = 0.1292 * exp(-pow((v+45.79)/15.54,2)) + 0.06487 * exp(-pow((v-4.823)/51.12,2)); +real dm = (mss - m) / taum; // Rush-Larsen + +// h gate +real ah = (v >= -40) ? (0) : (0.057 * exp( -(v + 80) / 6.8 )); +real bh = (v >= -40) ? (0.77 / (0.13*(1 + exp( -(v + 10.66) / 11.1 )))) : ((2.7 * exp( 0.079 * v) + 3.1*pow(10,5) * exp(0.3485 * v))); +real tauh = 1 / (ah + bh); +real hss = 1 / (pow(1 + exp( (v + 71.55)/7.43 ),2)); +real dh = (hss - h) / tauh; // Rush-Larsen +// j gate +real aj = (v >= -40) ? (0) : (((-2.5428 * pow(10,4)*exp(0.2444*v) - 6.948*pow(10,-6) * exp(-0.04391*v)) * (v + 37.78)) / (1 + exp( 0.311 * (v + 79.23) ))); +real bj = (v >= -40) ? ((0.6 * exp( 0.057 * v)) / (1 + exp( -0.1 * (v + 32) ))) : ((0.02424 * exp( -0.01052 * v )) / (1 + exp( -0.1378 * (v + 40.14) ))); +real tauj = 1 / (aj + bj); +real jss = 1 / pow((1 + exp( (v + 71.55)/7.43 )),2); +real dj = (jss - j) / tauj; // Rush-Larsen + +// h phosphorylated +real hssp = 1 / pow((1 + exp( (v + 71.55 + 6)/7.43 )),2); +real dhp = (hssp - hp) / tauh; // Rush-Larsen +// j phosphorylated +real taujp = 1.46 * tauj; +real djp = (jss - jp) / taujp; // Rush-Larsen +real GNa = 11.7802; +real INa = INa_Multiplier*GNa*(v-ENa)*pow(m,3.0)*((1.0-fINap)*h*j+fINap*hp*jp); + +// INaL +// calculate INaL +real mLss=1.0/(1.0+exp((-(v+42.85))/5.264)); +real tm = 0.1292 * exp(-pow(((v+45.79)/15.54),2)) + 0.06487 * exp(-pow(((v-4.823)/51.12),2)); +real tmL=tm; +real dmL=(mLss-mL)/tmL; // Rush-Larsen +real hLss=1.0/(1.0+exp((v+87.61)/7.488)); +real thL=200.0; +real dhL=(hLss-hL)/thL; // Rush-Larsen +real hLssp=1.0/(1.0+exp((v+93.81)/7.488)); +real thLp=3.0*thL; +real dhLp=(hLssp-hLp)/thLp; // Rush-Larsen +real GNaL=0.0279*INaL_Multiplier; +if (celltype==EPI) GNaL=GNaL*0.6; +real INaL = GNaL*(v-ENa)*mL*((1.0-fINaLp)*hL+fINaLp*hLp); + +// ITo +// calculate Ito +real ass=1.0/(1.0+exp((-(v-14.34))/14.82)); +real ta=1.0515/(1.0/(1.2089*(1.0+exp(-(v-18.4099)/29.3814)))+3.5/(1.0+exp((v+100.0)/29.3814))); +real da=(ass-a)/ta; // Rush-Larsen +real iss=1.0/(1.0+exp((v+43.94)/5.711)); +real delta_epi = (celltype == EPI) ? 1.0-(0.95/(1.0+exp((v+70.0)/5.0))) : 1.0; +real tiF=4.562+1/(0.3933*exp((-(v+100.0))/100.0)+0.08004*exp((v+50.0)/16.59)); +real tiS=23.62+1/(0.001416*exp((-(v+96.52))/59.05)+1.780e-8*exp((v+114.1)/8.079)); +tiF=tiF*delta_epi; +tiS=tiS*delta_epi; +real AiF=1.0/(1.0+exp((v-213.6)/151.2)); +real AiS=1.0-AiF; +real diF=(iss-iF)/tiF; // Rush-Larsen +real diS=(iss-iS)/tiS; // Rush-Larsen +real i=AiF*iF+AiS*iS; +real assp=1.0/(1.0+exp((-(v-24.34))/14.82)); +real dap=(assp-ap)/ta; // Rush-Larsen +real dti_develop=1.354+1.0e-4/(exp((v-167.4)/15.89)+exp(-(v-12.23)/0.2154)); +real dti_recover=1.0-0.5/(1.0+exp((v+70.0)/20.0)); +real tiFp=dti_develop*dti_recover*tiF; +real tiSp=dti_develop*dti_recover*tiS; +real diFp=(iss-iFp)/tiFp; // Rush-Larsen +real diSp=(iss-iSp)/tiSp; // Rush-Larsen +real ip=AiF*iFp+AiS*iSp; +real Gto=0.16*Ito_Multiplier; +Gto = (celltype == EPI || celltype == MID) ? Gto*2.0 : Gto; +real Ito = Gto*(v-EK)*((1.0-fItop)*a*i+fItop*ap*ip); + +// ICaL +// a variant updated by jakub, using a changed activation curve +// it computes both ICaL in subspace and myoplasm (_i) + +// calculate ICaL, ICaNa, ICaK +real dss=1.0763*exp(-1.0070*exp(-0.0829*(v))); // magyar +if(v >31.4978) dss = 1; // activation cannot be greater than 1 +real td= 0.6+1.0/(exp(-0.05*(v+6.0))+exp(0.09*(v+14.0))); +real dd=(dss-d)/td; // Rush-Larsen +real fss=1.0/(1.0+exp((v+19.58)/3.696)); +real tff=7.0+1.0/(0.0045*exp(-(v+20.0)/10.0)+0.0045*exp((v+20.0)/10.0)); +real tfs=1000.0+1.0/(0.000035*exp(-(v+5.0)/4.0)+0.000035*exp((v+5.0)/6.0)); +real Aff=0.6; +real Afs=1.0-Aff; +real dff=(fss-ff)/tff; // Rush-Larsen +real dfs=(fss-fs)/tfs; // Rush-Larsen +real f=Aff*ff+Afs*fs; +real fcass=fss; +real tfcaf=7.0+1.0/(0.04*exp(-(v-4.0)/7.0)+0.04*exp((v-4.0)/7.0)); +real tfcas=100.0+1.0/(0.00012*exp(-v/3.0)+0.00012*exp(v/7.0)); + +real Afcaf=0.3+0.6/(1.0+exp((v-10.0)/10.0)); +real Afcas=1.0-Afcaf; +real dfcaf=(fcass-fcaf)/tfcaf; // Rush-Larsen +real dfcas=(fcass-fcas)/tfcas; // Rush-Larsen +real fca=Afcaf*fcaf+Afcas*fcas; + +real tjca = 75; +//real tjca = 72.5; +real jcass = 1.0/(1.0+exp((v+18.08)/(2.7916))); +real djca=(jcass-jca)/tjca; // Rush-Larsen +real tffp=2.5*tff; +real dffp=(fss-ffp)/tffp; // Rush-Larsen +real fp=Aff*ffp+Afs*fs; +real tfcafp=2.5*tfcaf; +real dfcafp=(fcass-fcafp)/tfcafp; // Rush-Larsen +real fcap=Afcaf*fcafp+Afcas*fcas; + +// SS nca +real Kmn=0.002; +real k2n=500.0; +real km2n=jca*1; +real anca=1.0/(k2n/km2n+pow((1.0+Kmn/cass),4.0)); +real dnca=anca*k2n-nca*km2n; // Euler + +// myoplasmic nca +real anca_i = 1.0/(k2n/km2n+pow((1.0+Kmn/cai),4.0)); +real dnca_i = anca_i*k2n-nca_i*km2n; // Euler + +// SS driving force +//real clo = 150; // Extracellular Cl [mM] +//real cli = 24.0; // Intracellular Cl [mM] +real Io = 0.5*(nao + ko + clo + 4*cao)/1000; // ionic strength outside. /1000 is for things being in micromolar +real Ii = 0.5*(nass + kss + cli + 4*cass)/1000; // ionic strength outside. /1000 is for things being in micromolar +//real Ii = 0.5*(nass + kss + clss + 4*cass)/1000; // (dynCl) ionic strength outside. /1000 is for things being in micromolar +// The ionic strength is too high for basic DebHuc. We'll use Davies +real dielConstant = 74; // water at 37 degrees. +real temp = 310; // body temp in kelvins. +real constA = 1.82*pow(10,6)*pow((dielConstant*temp),(-1.5)); + +real gamma_cai = exp(-constA * 4 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii)); +real gamma_cao = exp(-constA * 4 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io)); +real gamma_nai = exp(-constA * 1 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii)); +real gamma_nao = exp(-constA * 1 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io)); +real gamma_ki = exp(-constA * 1 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii)); +real gamma_kao = exp(-constA * 1 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io)); + +real PhiCaL_ss = 4.0*vffrt*(gamma_cai*cass*exp(2.0*vfrt)-gamma_cao*cao)/(exp(2.0*vfrt)-1.0); +real PhiCaNa_ss = 1.0*vffrt*(gamma_nai*nass*exp(1.0*vfrt)-gamma_nao*nao)/(exp(1.0*vfrt)-1.0); +real PhiCaK_ss = 1.0*vffrt*(gamma_ki*kss*exp(1.0*vfrt)-gamma_kao*ko)/(exp(1.0*vfrt)-1.0); + +// Myo driving force +Io = 0.5*(nao + ko + clo + 4*cao)/1000; // ionic strength outside. /1000 is for things being in micromolar +Ii = 0.5*(nai + ki + cli + 4*cai)/1000; // ionic strength outside. /1000 is for things being in micromolar +// The ionic strength is too high for basic DebHuc. We'll use Davies + +gamma_cai = exp(-constA * 4 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii)); +gamma_cao = exp(-constA * 4 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io)); +gamma_nai = exp(-constA * 1 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii)); +gamma_nao = exp(-constA * 1 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io)); +gamma_ki = exp(-constA * 1 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii)); +gamma_kao = exp(-constA * 1 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io)); + +real gammaCaoMyo = gamma_cao; +real gammaCaiMyo = gamma_cai; + +real PhiCaL_i = 4.0*vffrt*(gamma_cai*cai*exp(2.0*vfrt)-gamma_cao*cao)/(exp(2.0*vfrt)-1.0); +real PhiCaNa_i = 1.0*vffrt*(gamma_nai*nai*exp(1.0*vfrt)-gamma_nao*nao)/(exp(1.0*vfrt)-1.0); +real PhiCaK_i = 1.0*vffrt*(gamma_ki*ki*exp(1.0*vfrt)-gamma_kao*ko)/(exp(1.0*vfrt)-1.0); +// The rest +real PCa=8.3757e-05 * ICaL_Multiplier; +if (celltype==EPI) + PCa=PCa*1.2; +else if (celltype==MID) + //PCa=PCa*2; + PCa = PCa*1.8; + +real PCap=1.1*PCa; +real PCaNa=0.00125*PCa; +real PCaK=3.574e-4*PCa; +real PCaNap=0.00125*PCap; +real PCaKp=3.574e-4*PCap; + +real ICaL_ss=(1.0-fICaLp)*PCa*PhiCaL_ss*d*(f*(1.0-nca)+jca*fca*nca)+fICaLp*PCap*PhiCaL_ss*d*(fp*(1.0-nca)+jca*fcap*nca); +real ICaNa_ss=(1.0-fICaLp)*PCaNa*PhiCaNa_ss*d*(f*(1.0-nca)+jca*fca*nca)+fICaLp*PCaNap*PhiCaNa_ss*d*(fp*(1.0-nca)+jca*fcap*nca); +real ICaK_ss=(1.0-fICaLp)*PCaK*PhiCaK_ss*d*(f*(1.0-nca)+jca*fca*nca)+fICaLp*PCaKp*PhiCaK_ss*d*(fp*(1.0-nca)+jca*fcap*nca); + +real ICaL_i=(1.0-fICaLp)*PCa*PhiCaL_i*d*(f*(1.0-nca_i)+jca*fca*nca_i)+fICaLp*PCap*PhiCaL_i*d*(fp*(1.0-nca_i)+jca*fcap*nca_i); +real ICaNa_i=(1.0-fICaLp)*PCaNa*PhiCaNa_i*d*(f*(1.0-nca_i)+jca*fca*nca_i)+fICaLp*PCaNap*PhiCaNa_i*d*(fp*(1.0-nca_i)+jca*fcap*nca_i); +real ICaK_i=(1.0-fICaLp)*PCaK*PhiCaK_i*d*(f*(1.0-nca_i)+jca*fca*nca_i)+fICaLp*PCaKp*PhiCaK_i*d*(fp*(1.0-nca_i)+jca*fcap*nca_i); + +// And we weight ICaL (in ss) and ICaL_i +ICaL_i = ICaL_i * (1-ICaL_fractionSS); +ICaNa_i = ICaNa_i * (1-ICaL_fractionSS); +ICaK_i = ICaK_i * (1-ICaL_fractionSS); +ICaL_ss = ICaL_ss * ICaL_fractionSS; +ICaNa_ss = ICaNa_ss * ICaL_fractionSS; +ICaK_ss = ICaK_ss * ICaL_fractionSS; + +real ICaL = ICaL_ss + ICaL_i; +real ICaNa = ICaNa_ss + ICaNa_i; +real ICaK = ICaK_ss + ICaK_i; +real ICaL_tot = ICaL + ICaNa + ICaK; + +// Ikr +// Variant based on Lu-Vandenberg +// Extracting state vector +// IMPORTANT: For reason of backward compatibility of naming of an older version of a MM IKr +// c3 in code is c0 in article diagram, c2 is c1, c1 is c2 +real c0 = ikr_c0; +real c1 = ikr_c1; +real c2 = ikr_c2; +//real c0 = ikr_c2; +//real c1 = ikr_c1; +//real c2 = ikr_c0; +real o = ikr_o; +real I = ikr_i; +real b = 0; // no channels blocked in via the mechanism of specific MM states + +// transition rates +// from c0 to c1 in l-v model, +real alpha = 0.1161 * exp(0.2990 * vfrt); +// from c1 to c0 in l-v/ +real beta = 0.2442 * exp(-1.604 * vfrt); + +// from c1 to c2 in l-v/ +real alpha1 = 1.25 * 0.1235; +// from c2 to c1 in l-v/ +real beta1 = 0.1911; + +// from c2 to o/ c1 to o +real alpha2 =0.0578 * exp(0.9710 * vfrt); +// from o to c2/ +real beta2 = 0.349e-3* exp(-1.062 * vfrt); + +// from o to i +real alphai = 0.2533 * exp(0.5953 * vfrt); +// from i to o +real betai = 1.25* 0.0522 * exp(-0.8209 * vfrt); + +// from c2 to i (from c1 in orig) +real alphac2ToI = 0.52e-4 * exp(1.525 * vfrt); +// from i to c2 +real betaItoC2 = 0.85e-8 * exp(-1.842 * vfrt); +betaItoC2 = (beta2 * betai * alphac2ToI)/(alpha2 * alphai); +// transitions themselves +// for reason of backward compatibility of naming of an older version of a +// MM IKr, c3 in code is c0 in article diagram, c2 is c1, c1 is c2. + +real dc0 = c1 * beta - c0 * alpha; // Euler +real dc1 = c0 * alpha + c2*beta1 - c1*(beta+alpha1); // Euler +real dc2 = c1 * alpha1 + o*beta2 + I*betaItoC2 - c2 * (beta1 + alpha2 + alphac2ToI); // Euler +real delta_o = c2 * alpha2 + I*betai - o*(beta2+alphai); // Euler +real di = c2*alphac2ToI + o*alphai - I*(betaItoC2 + betai); // Euler + +real GKr = 0.0321 * sqrt(ko/5) * IKr_Multiplier; // 1st element compensates for change to ko (sqrt(5/5.4)* 0.0362) +if (celltype==EPI) + GKr=GKr*1.3; +else if (celltype==MID) + GKr=GKr*0.8; + +real IKr = GKr * o * (v-EK); + +// calculate IKs +real xs1ss=1.0/(1.0+exp((-(v+11.60))/8.932)); +real txs1=817.3+1.0/(2.326e-4*exp((v+48.28)/17.80)+0.001292*exp((-(v+210.0))/230.0)); +real dxs1=(xs1ss-xs1)/txs1; // Rush-Larsen +real xs2ss=xs1ss; +real txs2=1.0/(0.01*exp((v-50.0)/20.0)+0.0193*exp((-(v+66.54))/31.0)); +real dxs2=(xs2ss-xs2)/txs2; // Rush-Larsen +real KsCa=1.0+0.6/(1.0+pow((3.8e-5/cai),1.4)); +real GKs= 0.0011*IKs_Multiplier; +if (celltype==EPI) + GKs=GKs*1.4; +real IKs = GKs*KsCa*xs1*xs2*(v-EKs); + +// IK1 +real aK1 = 4.094/(1+exp(0.1217*(v-EK-49.934))); +real bK1 = (15.72*exp(0.0674*(v-EK-3.257))+exp(0.0618*(v-EK-594.31)))/(1+exp(-0.1629*(v-EK+14.207))); +real K1ss = aK1/(aK1+bK1); +real GK1=IK1_Multiplier * 0.6992; // 0.7266; // * sqrt(5/5.4)) +if (celltype==EPI) + GK1=GK1*1.2; +else if (celltype==MID) + GK1=GK1*1.3; +real IK1=GK1*sqrt(ko/5)*K1ss*(v-EK); + +// IKCa +real fIKCass = 0.8; +real kdikca = 6.05e-4; +real ikcan = 3.5; +real GKCa = 0.003 * IKCa_Multiplier; +real IKCa_ss = GKCa * fIKCass * pow(cass,ikcan) / (pow(cass,ikcan) + pow(kdikca,ikcan)) * (v-EK); +real IKCa_i = GKCa * (1.0-fIKCass) * pow(cai,ikcan) / (pow(cai,ikcan) + pow(kdikca,ikcan)) * (v-EK); +real IKCa = IKCa_ss + IKCa_i; + +// INaCa +real zca = 2.0; +real kna1=15.0; +real kna2=5.0; +real kna3=88.12; +real kasymm=12.5; +real wna=6.0e4; +real wca=6.0e4; +real wnaca=5.0e3; +real kcaon=1.5e6; +real kcaoff=5.0e3; +real qna=0.5224; +real qca=0.1670; +real hca=exp((qca*v*F)/(R*T)); +real hna=exp((qna*v*F)/(R*T)); + +real h1=1+nai/kna3*(1+hna); +real h2=(nai*hna)/(kna3*h1); +real h3=1.0/h1; +real h4=1.0+nai/kna1*(1+nai/kna2); +real h5=nai*nai/(h4*kna1*kna2); +real h6=1.0/h4; +real h7=1.0+nao/kna3*(1.0+1.0/hna); +real h8=nao/(kna3*hna*h7); +real h9=1.0/h7; +real h10=kasymm+1.0+nao/kna1*(1.0+nao/kna2); +real h11=nao*nao/(h10*kna1*kna2); +real h12=1.0/h10; +real k1=h12*cao*kcaon; +real k2=kcaoff; +real k3p=h9*wca; +real k3pp=h8*wnaca; +real k3=k3p+k3pp; +real k4p=h3*wca/hca; +real k4pp=h2*wnaca; +real k4=k4p+k4pp; +real k5=kcaoff; +real k6=h6*cai*kcaon; +real k7=h5*h2*wna; +real k8=h8*h11*wna; +real x1=k2*k4*(k7+k6)+k5*k7*(k2+k3); +real x2=k1*k7*(k4+k5)+k4*k6*(k1+k8); +real x3=k1*k3*(k7+k6)+k8*k6*(k2+k3); +real x4=k2*k8*(k4+k5)+k3*k5*(k1+k8); +real E1=x1/(x1+x2+x3+x4); +real E2=x2/(x1+x2+x3+x4); +real E3=x3/(x1+x2+x3+x4); +real E4=x4/(x1+x2+x3+x4); +real KmCaAct=150.0e-6; +real allo=1.0/(1.0+pow((KmCaAct/cai),2.0)); +real zna=1.0; +real JncxNa=3.0*(E4*k7-E1*k8)+E3*k4pp-E2*k3pp; +real JncxCa=E2*k2-E1*k1; +real Gncx= 0.0034 * INaCa_Multiplier; +if (celltype==EPI) + Gncx=Gncx*1.1; +else if (celltype==MID) + Gncx=Gncx*1.4; +real INaCa_i = (1-INaCa_fractionSS)*Gncx*allo*(zna*JncxNa+zca*JncxCa); + +// calculate INaCa_ss +h1=1+nass/kna3*(1+hna); +h2=(nass*hna)/(kna3*h1); +h3=1.0/h1; +h4=1.0+nass/kna1*(1+nass/kna2); +h5=nass*nass/(h4*kna1*kna2); +h6=1.0/h4; +h7=1.0+nao/kna3*(1.0+1.0/hna); +h8=nao/(kna3*hna*h7); +h9=1.0/h7; +h10=kasymm+1.0+nao/kna1*(1+nao/kna2); +h11=nao*nao/(h10*kna1*kna2); +h12=1.0/h10; +k1=h12*cao*kcaon; +k2=kcaoff; +k3p=h9*wca; +k3pp=h8*wnaca; +k3=k3p+k3pp; +k4p=h3*wca/hca; +k4pp=h2*wnaca; +k4=k4p+k4pp; +k5=kcaoff; +k6=h6*cass*kcaon; +k7=h5*h2*wna; +k8=h8*h11*wna; +x1=k2*k4*(k7+k6)+k5*k7*(k2+k3); +x2=k1*k7*(k4+k5)+k4*k6*(k1+k8); +x3=k1*k3*(k7+k6)+k8*k6*(k2+k3); +x4=k2*k8*(k4+k5)+k3*k5*(k1+k8); +E1=x1/(x1+x2+x3+x4); +E2=x2/(x1+x2+x3+x4); +E3=x3/(x1+x2+x3+x4); +E4=x4/(x1+x2+x3+x4); +KmCaAct=150.0e-6 ; +allo=1.0/(1.0+pow((KmCaAct/cass),2.0)); +JncxNa=3.0*(E4*k7-E1*k8)+E3*k4pp-E2*k3pp; +JncxCa=E2*k2-E1*k1; +real INaCa_ss = INaCa_fractionSS*Gncx*allo*(zna*JncxNa+zca*JncxCa); + +// calculate INaK +real k1p=949.5; +real k1m=182.4; +real k2p=687.2; +real k2m=39.4; +k3p=1899.0; +real k3m=79300.0; +k4p=639.0; +real k4m=40.0; +real Knai0=9.073; +real Knao0=27.78; +real delta=-0.1550; +real Knai=Knai0*exp((delta*v*F)/(3.0*R*T)); +real Knao=Knao0*exp(((1.0-delta)*v*F)/(3.0*R*T)); +real Kki=0.5; +real Kko=0.3582; +real MgADP=0.05; +real MgATP=9.8; +real Kmgatp=1.698e-7; +real H=1.0e-7; +real eP=4.2; +real Khp=1.698e-7; +real Knap=224.0; +real Kxkur=292.0; +real P=eP/(1.0+H/Khp+nai/Knap+ki/Kxkur); +real a1=(k1p*pow((nai/Knai),3.0))/(pow((1.0+nai/Knai),3.0)+pow((1.0+ki/Kki),2.0)-1.0); +real b1=k1m*MgADP; +real a2=k2p; +real b2=(k2m*pow((nao/Knao),3.0))/(pow((1.0+nao/Knao),3.0)+pow((1.0+ko/Kko),2.0)-1.0); +real a3=(k3p*pow((ko/Kko),2.0))/(pow((1.0+nao/Knao),3.0)+pow((1.0+ko/Kko),2.0)-1.0); +real b3=(k3m*P*H)/(1.0+MgATP/Kmgatp); +real a4=(k4p*MgATP/Kmgatp)/(1.0+MgATP/Kmgatp); +real b4=(k4m*pow((ki/Kki),2.0))/(pow((1.0+nai/Knai),3.0)+pow((1.0+ki/Kki),2.0)-1.0); +x1=a4*a1*a2+b2*b4*b3+a2*b4*b3+b3*a1*a2; +x2=b2*b1*b4+a1*a2*a3+a3*b1*b4+a2*a3*b4; +x3=a2*a3*a4+b3*b2*b1+b2*b1*a4+a3*a4*b1; +x4=b4*b3*b2+a3*a4*a1+b2*a4*a1+b3*b2*a1; +E1=x1/(x1+x2+x3+x4); +E2=x2/(x1+x2+x3+x4); +E3=x3/(x1+x2+x3+x4); +E4=x4/(x1+x2+x3+x4); +real zk=1.0; +real JnakNa=3.0*(E1*a3-E2*b3); +real JnakK=2.0*(E4*b1-E3*a1); +real Pnak= 15.4509; +if (celltype==EPI) + Pnak=Pnak*0.9; +else if (celltype==MID) + Pnak=Pnak*0.7; + +real INaK = Pnak*(zna*JnakNa+zk*JnakK)*INaK_Multiplier; + +// Minor/background currents +// calculate IKb +real xkb=1.0/(1.0+exp(-(v-10.8968)/(23.9871))); +real GKb=0.0189*IKb_Multiplier; + +if (IKCa_Multiplier > 0.0) + GKb = GKb*0.9; +if (celltype==EPI) + GKb=GKb*0.6; +real IKb = GKb*xkb*(v-EK); + +// calculate INab +real PNab=1.9239e-09*INab_Multiplier; +real INab=PNab*vffrt*(nai*exp(vfrt)-nao)/(exp(vfrt)-1.0); + +// calculate ICab +real PCab=5.9194e-08*ICab_Multiplier; +real ICab=PCab*4.0*vffrt*(gammaCaiMyo*cai*exp(2.0*vfrt)-gammaCaoMyo*cao)/(exp(2.0*vfrt)-1.0); + +// calculate IpCa +real GpCa=5e-04*IpCa_Multiplier; +real IpCa=GpCa*cai/(0.0005+cai); + +// Chloride +// I_ClCa: Ca-activated Cl Current, I_Clbk: background Cl Current + +real ECl = (R*T/F)*log(cli/clo); // [mV] +//real EClss = (R*T/F)*log(clss/clo); // [mV] dynCl + +real Fjunc = 1; +real Fsl = 1-Fjunc; // fraction in SS and in myoplasm - as per literature, I(Ca)Cl is in junctional subspace + +real GClCa = ICaCl_Multiplier * 0.2843; // [mS/uF] +real GClB = IClb_Multiplier * 1.98e-3; // [mS/uF] +real KdClCa = 0.1; // [mM] + +//real I_ClCa_junc = Fjunc*GClCa/(1+KdClCa/cass)*(v-EClss); // dynCl +//real I_ClCa_sl = Fsl*GClCa/(1+KdClCa/cai)*(v-ECl); // dynCl +real I_ClCa_junc = Fjunc*GClCa/(1+KdClCa/cass)*(v-ECl); +real I_ClCa_sl = Fsl*GClCa/(1+KdClCa/cai)*(v-ECl); + +real I_ClCa = I_ClCa_junc+I_ClCa_sl; +real I_Clbk = GClB*(v-ECl); + +// Calcium handling +// calculate ryanodione receptor calcium induced calcium release from the jsr +real fJrelp=(1.0/(1.0+KmCaMK/CaMKa)); + +// Jrel +real jsrMidpoint = 1.7; + +real bt=4.75; +real a_rel=0.5*bt; +real Jrel_inf=a_rel*(-ICaL)/(1.0+pow((jsrMidpoint/cajsr),8.0)); +if (celltype==MID) + Jrel_inf=Jrel_inf*1.7; +real tau_rel=bt/(1.0+0.0123/cajsr); + +if (tau_rel<0.001) + tau_rel=0.001; + +real dJrelnp=(Jrel_inf-Jrel_np)/tau_rel; // Rush-Larsen + +real btp=1.25*bt; +real a_relp=0.5*btp; +real Jrel_infp=a_relp*(-ICaL)/(1.0+pow((jsrMidpoint/cajsr),8.0)); +if (celltype==MID) + Jrel_infp=Jrel_infp*1.7; +real tau_relp=btp/(1.0+0.0123/cajsr); + +if (tau_relp<0.001) + tau_relp=0.001; + +tau_relp = tau_relp*taurelp_Multiplier; +real dJrelp=(Jrel_infp-Jrel_p)/tau_relp; // Rush-Larsen +real Jrel = 1.5378 * ((1.0-fJrelp)*Jrel_np+fJrelp*Jrel_p)*Jrel_Multiplier; + +real fJupp=(1.0/(1.0+KmCaMK/CaMKa)); + +// Jup +// calculate serca pump, ca uptake flux +// camkFactor = 2.4; +// gjup = 0.00696; +// Jupnp=Jup_Multiplier * gjup*cai/(cai+0.001); +// Jupp=Jup_Multiplier * camkFactor*gjup*cai/(cai + 8.2500e-04); +// if celltype==1 +// Jupnp=Jupnp*1.3; +// Jupp=Jupp*1.3; +// end +// +// +// Jleak=Jup_Multiplier * 0.00629 * cansr/15.0; +// Jup=(1.0-fJupp)*Jupnp+fJupp*Jupp-Jleak; + +// calculate serca pump, ca uptake flux +real Jupnp = Jup_Multiplier * 0.005425*cai/(cai+0.00092); +real Jupp = Jup_Multiplier * 2.75*0.005425*cai/(cai+0.00092-0.00017); +if (celltype==EPI) { + Jupnp=Jupnp*1.3; + Jupp=Jupp*1.3; +} +real Jleak=Jup_Multiplier* 0.0048825*cansr/15.0; +real Jup=(1.0-fJupp)*Jupnp+fJupp*Jupp-Jleak; + +//calculate tranlocation flux +real Jtr=(cansr-cajsr)/60; + +// I_katp current (fkatp current) +//real I_katp = (fkatp * gkatp * akik * bkik * (v - EK)); + +// stimulus current +real Istim = calc_I_stim; + +//update the membrane voltage +//real dv=-(INa+INaL+Ito+ICaL+ICaNa+ICaK+IKr+IKs+IK1+INaCa_i+INaCa_ss+INaK+INab+IKb+IpCa+ICab+I_ClCa+I_Clbk+Istim); // Euler +real dv=-(INa+INaL+Ito+ICaL+ICaNa+ICaK+IKr+IKs+IK1+IKCa+INaCa_i+INaCa_ss+INaK+INab+IKb+IpCa+ICab+I_ClCa+I_Clbk+Istim); // Euler +//real dv=-(INa+INaL+Ito+ICaL+ICaNa+ICaK+IKr+IKs+IK1+INaCa_i+INaCa_ss+INaK+INab+IKb+IpCa+ICab+I_ClCa+I_Clbk+I_katp+Istim); // Euler + +// calculate diffusion fluxes +real JdiffNa=(nass-nai)/2.0; +real JdiffK=(kss-ki)/2.0; +real Jdiff=(cass-cai)/0.2; +//real JdiffCl=(clss-cli)/2.0; // dynCl + +// calcium buffer constants +real cmdnmax= 0.05; +if (celltype==EPI) + cmdnmax=cmdnmax*1.3; +real kmcmdn=0.00238; +//real trpnmax=0.07; +real kmtrpn=0.0005; +real BSRmax=0.047; +real KmBSR = 0.00087; +real BSLmax=1.124; +real KmBSL = 0.0087; +real csqnmax=10.0; +real kmcsqn=0.8; + +// update intracellular concentrations, using buffers for cai, cass, cajsr +real dnai=-(ICaNa_i+INa+INaL+3.0*INaCa_i+3.0*INaK+INab)*Acap/(F*vmyo)+JdiffNa*vss/vmyo; // Euler +real dnass=-(ICaNa_ss+3.0*INaCa_ss)*Acap/(F*vss)-JdiffNa; // Euler + +real dki=-(ICaK_i+Ito+IKr+IKs+IK1+IKb+Istim-2.0*INaK)*Acap/(F*vmyo)+JdiffK*vss/vmyo; // Euler +real dkss=-(ICaK_ss)*Acap/(F*vss)-JdiffK; // Euler + +//real Bcai=1.0/(1.0+cmdnmax*kmcmdn/pow((kmcmdn+cai),2.0)+trpnmax*kmtrpn/pow((kmtrpn+cai),2.0)); // dynCl +real Bcai=1.0/(1.0+cmdnmax*kmcmdn/pow((kmcmdn+cai),2.0)); +//real dcai=Bcai*(-(ICaL_i + IpCa+ICab-2.0*INaCa_i)*Acap/(2.0*F*vmyo)-Jup*vnsr/vmyo+Jdiff*vss/vmyo); // dynCl +real dcai=Bcai*(-(ICaL_i + IpCa+ICab-2.0*INaCa_i)*Acap/(2.0*F*vmyo)-Jup*vnsr/vmyo+Jdiff*vss/vmyo-dCa_TRPN*trpnmax); + +real Bcass=1.0/(1.0+BSRmax*KmBSR/pow((KmBSR+cass),2.0)+BSLmax*KmBSL/pow((KmBSL+cass),2.0)); +real dcass=Bcass*(-(ICaL_ss-2.0*INaCa_ss)*Acap/(2.0*F*vss)+Jrel*vjsr/vss-Jdiff); + +real dcansr=Jup-Jtr*vjsr/vnsr; + +real Bcajsr=1.0/(1.0+csqnmax*kmcsqn/pow((kmcsqn+cajsr),2.0)); +real dcajsr=Bcajsr*(Jtr-Jrel); + +//real dcli = - (I_Clbk + I_ClCa_sl)*Acap/(-1*F*vmyo)+JdiffCl*vss/vmyo; // dynCl +//real dclss = - I_ClCa_junc*Acap/(-1*F*vss)-JdiffCl; // dynCl + +// Compute 'a' coefficients for the Hodkin-Huxley variables +a_[ 9] = -1.0 / tm; +a_[10] = -1.0 / tauh; +a_[11] = -1.0 / tauh; +a_[12] = -1.0 / tauj; +a_[13] = -1.0 / taujp; +a_[14] = -1.0 / tmL; +a_[15] = -1.0 / thL; +a_[16] = -1.0 / thLp; +a_[17] = -1.0 / ta; +a_[18] = -1.0 / tiF; +a_[19] = -1.0 / tiS; +a_[20] = -1.0 / ta; +a_[21] = -1.0 / tiFp; +a_[22] = -1.0 / tiSp; +a_[23] = -1.0 / td; +a_[24] = -1.0 / tff; +a_[25] = -1.0 / tfs; +a_[26] = -1.0 / tfcaf; +a_[27] = -1.0 / tfcas; +a_[28] = -1.0 / tjca; +a_[31] = -1.0 / tffp; +a_[32] = -1.0 / tfcafp; +a_[33] = -1.0 / txs1; +a_[34] = -1.0 / txs2; +a_[35] = -1.0 / tau_relp; +a_[42] = -1.0 / tau_rel; + +// Compute 'b' coefficients for the Hodkin-Huxley variables +b_[ 9] = mss / tm; +b_[10] = hssp / tauh; +b_[11] = hss / tauh; +b_[12] = jss / tauj; +b_[13] = jss / taujp; +b_[14] = mLss / tmL; +b_[15] = hLss / thL; +b_[16] = hLssp / thLp; +b_[17] = ass / ta; +b_[18] = iss / tiF; +b_[19] = iss / tiS; +b_[20] = assp / ta; +b_[21] = iss / tiFp; +b_[22] = iss / tiSp; +b_[23] = dss / td; +b_[24] = fss / tff; +b_[25] = fss / tfs; +b_[26] = fcass / tfcaf; +b_[27] = fcass / tfcas; +b_[28] = jcass / tjca; +b_[31] = fss / tffp; +b_[32] = fcass / tfcafp; +b_[33] = xs1ss / txs1; +b_[34] = xs2ss / txs2; +b_[35] = Jrel_infp / tau_relp; +b_[42] = Jrel_inf / tau_rel; + +// Right-hand side +rDY_[0] = dv; +rDY_[1] = dnai; +rDY_[2] = dnass; +rDY_[3] = dki; +rDY_[4] = dkss; +rDY_[5] = dcai; +rDY_[6] = dcass; +rDY_[7] = dcansr; +rDY_[8] = dcajsr; +rDY_[9] = dm; +rDY_[10] = dhp; +rDY_[11] = dh; +rDY_[12] = dj; +rDY_[13] = djp; +rDY_[14] = dmL; +rDY_[15] = dhL; +rDY_[16] = dhLp; +rDY_[17] = da; +rDY_[18] = diF; +rDY_[19] = diS; +rDY_[20] = dap; +rDY_[21] = diFp; +rDY_[22] = diSp; +rDY_[23] = dd; +rDY_[24] = dff; +rDY_[25] = dfs; +rDY_[26] = dfcaf; +rDY_[27] = dfcas; +rDY_[28] = djca; +rDY_[29] = dnca; +rDY_[30] = dnca_i; +rDY_[31] = dffp; +rDY_[32] = dfcafp; +rDY_[33] = dxs1; +rDY_[34] = dxs2; +rDY_[35] = dJrelnp; +rDY_[36] = dCaMKt; +rDY_[37] = dc0; +rDY_[38] = dc1; +rDY_[39] = dc2; +rDY_[40] = delta_o; +rDY_[41] = di; +rDY_[42] = dJrelp; +// ----------------------- +// Land-Niederer +rDY_[43] = dXS; +rDY_[44] = dXW; +rDY_[45] = dCa_TRPN; +rDY_[46] = dTmBlocked; +rDY_[47] = dZETAS; +rDY_[48] = dZETAW; diff --git a/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi.c b/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi.c index 1b1afe25..ff2e3448 100644 --- a/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi.c +++ b/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi.c @@ -228,12 +228,12 @@ SOLVE_MODEL_ODES(solve_model_odes_cpu) { if(adpt) { if (ode_solver->ode_extra_data) { - //solve_forward_euler_cpu_adpt(sv + (sv_id * NEQ), stim_currents[i], transmurality[i], current_t + dt, sv_id, ode_solver, extra_par); - solve_rush_larsen_cpu_adpt(sv + (sv_id * NEQ), stim_currents[i], transmurality[i], current_t + dt, sv_id, ode_solver, extra_par); + solve_forward_euler_cpu_adpt(sv + (sv_id * NEQ), stim_currents[i], transmurality[i], current_t + dt, sv_id, ode_solver, extra_par); + //solve_rush_larsen_cpu_adpt(sv + (sv_id * NEQ), stim_currents[i], transmurality[i], current_t + dt, sv_id, ode_solver, extra_par); } else { - //solve_forward_euler_cpu_adpt(sv + (sv_id * NEQ), stim_currents[i], 0.0, current_t + dt, sv_id, ode_solver, extra_par); - solve_rush_larsen_cpu_adpt(sv + (sv_id * NEQ), stim_currents[i], 0.0, current_t + dt, sv_id, ode_solver, extra_par); + solve_forward_euler_cpu_adpt(sv + (sv_id * NEQ), stim_currents[i], 0.0, current_t + dt, sv_id, ode_solver, extra_par); + //solve_rush_larsen_cpu_adpt(sv + (sv_id * NEQ), stim_currents[i], 0.0, current_t + dt, sv_id, ode_solver, extra_par); } } else { diff --git a/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi.cu b/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi.cu index b45321bb..1cb8a96e 100644 --- a/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi.cu +++ b/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi.cu @@ -377,8 +377,8 @@ __global__ void solve_gpu(real cur_time, real dt, real *sv, real *stim_currents, SOLVE_EQUATION_RUSH_LARSEN_GPU(42); // Jrel_p } } else { - //solve_forward_euler_gpu_adpt(sv, stim_currents[threadID], 0.0, extra_params, cur_time + max_dt, sv_id, pitch, abstol, reltol, dt, max_dt); - solve_rush_larsen_gpu_adpt(sv, stim_currents[threadID], 0.0, extra_params, cur_time + max_dt, sv_id, pitch, abstol, reltol, dt, max_dt); + solve_forward_euler_gpu_adpt(sv, stim_currents[threadID], 0.0, extra_params, cur_time + max_dt, sv_id, pitch, abstol, reltol, dt, max_dt); + //solve_rush_larsen_gpu_adpt(sv, stim_currents[threadID], 0.0, extra_params, cur_time + max_dt, sv_id, pitch, abstol, reltol, dt, max_dt); } } } diff --git a/src/models_library/ToROrd/build.sh b/src/models_library/ToROrd/build.sh index e069da97..0da7b92e 100644 --- a/src/models_library/ToROrd/build.sh +++ b/src/models_library/ToROrd/build.sh @@ -18,3 +18,10 @@ MODEL_FILE_GPU="ToRORd_dynCl_mixed_endo_mid_epi.cu" COMMON_HEADERS="ToRORd_dynCl_mixed_endo_mid_epi.h" COMPILE_MODEL_LIB "ToRORd_dynCl_mixed_endo_mid_epi" "$MODEL_FILE_CPU" "$MODEL_FILE_GPU" "$COMMON_HEADERS" + +############## ToRORd Land Mixed ENDO_MID_EPI ############################## +MODEL_FILE_CPU="ToRORd_Land_mixed_endo_mid_epi.c" +MODEL_FILE_GPU="ToRORd_Land_mixed_endo_mid_epi.cu" +COMMON_HEADERS="ToRORd_Land_mixed_endo_mid_epi.h" + +COMPILE_MODEL_LIB "ToRORd_Land_mixed_endo_mid_epi" "$MODEL_FILE_CPU" "$MODEL_FILE_GPU" "$COMMON_HEADERS" 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.c b/src/save_mesh_library/save_mesh.c index 9d44896b..7ebaec6e 100644 --- a/src/save_mesh_library/save_mesh.c +++ b/src/save_mesh_library/save_mesh.c @@ -156,72 +156,45 @@ SAVE_MESH(save_one_cell_state_variables) { check_cuda_error(cudaMemcpy2D(cell_sv, sizeof(real), ode_solver->sv + params->cell_sv_position, ode_solver->pitch, sizeof(real), ode_solver->model_data.number_of_ode_equations, cudaMemcpyDeviceToHost)); - fprintf(params->file, "%lf %lf ", time_info->current_t, cell_sv[0]); - for(int i = 2; i < 11; i++) { - fprintf(params->file, " %lf ", cell_sv[i]); + // All state variables + for (uint32_t i = 0; i < ode_solver->model_data.number_of_ode_equations; i++) { + fprintf(params->file, "%g, ", cell_sv[i]); } + fprintf(params->file, "\n"); - fprintf(params->file, " %lf ", cell_sv[13]); - fprintf(params->file, " %lf ", cell_sv[11]); - fprintf(params->file, " %lf ", cell_sv[12]); - - for(int i = 14; i < 30; i++) { - fprintf(params->file, " %lf ", cell_sv[i]); - } - - fprintf(params->file, " %lf ", cell_sv[32]); - fprintf(params->file, " %lf ", cell_sv[33]); - - fprintf(params->file, " %lf ", cell_sv[30]); - fprintf(params->file, " %lf ", cell_sv[31]); - - fprintf(params->file, " %lf ", cell_sv[39]); - fprintf(params->file, " %lf ", cell_sv[40]); - fprintf(params->file, " %lf ", cell_sv[41]); - fprintf(params->file, " %lf ", cell_sv[1]); - - fprintf(params->file, " %lf ", cell_sv[34]); - fprintf(params->file, " %lf ", cell_sv[35]); - fprintf(params->file, " %lf ", cell_sv[36]); - fprintf(params->file, " %lf ", cell_sv[37]); - fprintf(params->file, " %lf ", cell_sv[38]); - fprintf(params->file, " %lf\n", cell_sv[42]); free(cell_sv); #endif } else { real *cell_sv = &ode_solver->sv[params->cell_sv_position * ode_solver->model_data.number_of_ode_equations]; - fprintf(params->file, "%lf %lf ", time_info->current_t, cell_sv[0]); - for(int i = 2; i < 11; i++) { - fprintf(params->file, " %lf ", cell_sv[i]); - } + // Time and transmembrane potential + //fprintf(params->file, "%g %g\n", time_info->current_t, cell_sv[0]); + + // Time, Cai and Vm + //fprintf(params->file, "%g %g %g\n", time_info->current_t, cell_sv[0], cell_sv[5]); - fprintf(params->file, " %lf ", cell_sv[13]); - fprintf(params->file, " %lf ", cell_sv[11]); - fprintf(params->file, " %lf ", cell_sv[12]); + // Only transmembrane potential + //fprintf(params->file, "%g\n", cell_sv[0]); - for(int i = 14; i < 30; i++) { - fprintf(params->file, " %lf ", cell_sv[i]); + // All state variables + for (uint32_t i = 0; i < ode_solver->model_data.number_of_ode_equations; i++) { + fprintf(params->file, "%g, ", cell_sv[i]); } - - fprintf(params->file, " %lf ", cell_sv[32]); - fprintf(params->file, " %lf ", cell_sv[33]); - - fprintf(params->file, " %lf ", cell_sv[30]); - fprintf(params->file, " %lf ", cell_sv[31]); - - fprintf(params->file, " %lf ", cell_sv[39]); - fprintf(params->file, " %lf ", cell_sv[40]); - fprintf(params->file, " %lf ", cell_sv[41]); - fprintf(params->file, " %lf ", cell_sv[1]); - - fprintf(params->file, " %lf ", cell_sv[34]); - fprintf(params->file, " %lf ", cell_sv[35]); - fprintf(params->file, " %lf ", cell_sv[36]); - fprintf(params->file, " %lf ", cell_sv[37]); - fprintf(params->file, " %lf ", cell_sv[38]); - fprintf(params->file, " %lf\n", cell_sv[42]); + fprintf(params->file, "\n"); + + // Time, Cai and Vm at certain timestep + //if (time_info->current_t >= 15200) { + // fprintf(params->file, "%.3lf %g %g\n", time_info->current_t, cell_sv[0], cell_sv[5]); + //} + + // All state-variables at certain timestep + //if (time_info->current_t >= 15200) { + // for (uint32_t i = 0; i < ode_solver->model_data.number_of_ode_equations; i++) { + // fprintf(params->file, "%g, ", cell_sv[i]); + // } + // fprintf(params->file, "\n"); + //} } } 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) {