diff --git a/scripts/tuneCV/main.cpp b/scripts/tuneCV/main.cpp index 658be3b8..8ff9440c 100644 --- a/scripts/tuneCV/main.cpp +++ b/scripts/tuneCV/main.cpp @@ -11,10 +11,6 @@ #include #include #include -#include -#include -#include -#include #include #include #include @@ -23,57 +19,87 @@ #include #include #include +#include using namespace std; +// Change your MonoAlg3D path here: +// ---------------------------------------------------------- +const char MONOALG_PATH[500] = "/home/berg/Github/MonoAlg3D_C"; +// ---------------------------------------------------------- + const double TOLERANCE = 1.0e-03; double calculate_conduction_velocity_from_simulation () { - string filename = "outputs/cable/purkinje_activation_time_map_pulse_it_0.vtp"; + string filename = "outputs/cable/tissue_activation_time_map_pulse_it_0.vtu"; // Read all the data from the file - vtkSmartPointer reader = vtkSmartPointer::New(); + vtkSmartPointer reader = vtkSmartPointer::New(); reader->SetFileName(filename.c_str()); reader->Update(); - vtkPolyData* polydata = reader->GetOutput(); - uint32_t num_points = polydata->GetNumberOfPoints(); - uint32_t num_lines = polydata->GetNumberOfLines(); + vtkUnstructuredGrid* unstructuredGrid = reader->GetOutput(); + uint32_t num_points = unstructuredGrid->GetNumberOfPoints(); + uint32_t num_lines = unstructuredGrid->GetNumberOfCells(); + + vtkSmartPointer cellLocator = vtkSmartPointer::New(); + cellLocator->SetDataSet(unstructuredGrid); + cellLocator->BuildLocator(); + // Points of interest + double point_x_0[3] = {7500, 250, 250}; + double point_x_1[3] = {12500, 250, 250}; + + // Find (closest points) Cell indexes for CV computation + vtkIdType cellId_x_0; // the cell id of the cell containing the closest point + vtkIdType cellId_x_1; // the cell id of the cell containing the closest point + + double closestPoint_x_0[3]; + double closestPoint_x_1[3]; + + int subId; // not needed + double closestPointDist2; // not needed + + cellLocator->FindClosestPoint(point_x_0, closestPoint_x_0, cellId_x_0, subId, closestPointDist2); + cellLocator->FindClosestPoint(point_x_1, closestPoint_x_1, cellId_x_1, subId, closestPointDist2); + + double delta_s_x = sqrt(pow(closestPoint_x_0[0]-closestPoint_x_1[0], 2)); // only changes over x-axis + + cout << delta_s_x << endl; + // Read points scalar values string array_name = "Scalars_"; - vtkSmartPointer array = vtkFloatArray::SafeDownCast(polydata->GetPointData()->GetArray(array_name.c_str())); - - double cv = 0.0; - if(array) - { - // Cell indexes for CV computation - uint32_t middle_cell = num_points/2; - uint32_t prev_cell = middle_cell-25; - uint32_t next_cell = middle_cell+25; + vtkSmartPointer array = vtkFloatArray::SafeDownCast(unstructuredGrid->GetCellData()->GetArray(array_name.c_str())); + + double cv_x = -1.0; + + if(array) { + double delta_lat_x = (array->GetValue(cellId_x_1) - array->GetValue(cellId_x_0)); // ms - double delta_s = 5000.0; - double delta_t = (array->GetValue(next_cell) - array->GetValue(prev_cell)); - cv = (delta_s / delta_t)*0.001; // {m/s} + cout << delta_lat_x << endl; + cout << array->GetValue(cellId_x_0) << endl; + cout << array->GetValue(cellId_x_1) << endl; + + cv_x = (delta_s_x / delta_lat_x)*0.001; // {m/s} } - else - { - cerr << "[!] ERROR! No Scalar_value found for the points!" << endl; + else { + cerr << "[!] ERROR! No 'Scalar_value' found for the points!" << endl; exit(EXIT_FAILURE); } - return cv; + return cv_x; } // TODO: Maybe pass a pre-configured config file as an input parameter with the cellular model setup that the user will use -void write_configuration_file (const double sigma) -{ - FILE *file = fopen("/home/berg/Github/MonoAlg3D_C/scripts/tuneCV/configs/cable.ini","w+"); +void write_configuration_file (const double sigma) { + char filename[500]; + sprintf(filename, "%s/scripts/tuneCV/configs/cable.ini", MONOALG_PATH); + FILE *file = fopen(filename,"w+"); fprintf(file,"[main]\n"); fprintf(file,"num_threads=4\n"); - fprintf(file,"dt_pde=0.02\n"); + fprintf(file,"dt_pde=0.01\n"); fprintf(file,"simulation_time=150.0\n"); fprintf(file,"abort_on_no_activity=false\n"); fprintf(file,"use_adaptivity=false\n"); @@ -82,80 +108,83 @@ void write_configuration_file (const double sigma) fprintf(file,"[update_monodomain]\n"); fprintf(file,"main_function=update_monodomain_default\n"); - fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libdefault_update_monodomain.so\n"); + fprintf(file,"library_file=%s/shared_libs/libdefault_update_monodomain.so\n", MONOALG_PATH); fprintf(file,"\n"); + // For saving the LATs in a format that can be read for calculating the CVs fprintf(file,"[save_result]\n"); fprintf(file,"print_rate=1\n"); - fprintf(file,"output_dir=/home/berg/Github/MonoAlg3D_C/scripts/tuneCV/outputs/cable\n"); + fprintf(file,"output_dir=%s/scripts/tuneCV/outputs/cable\n", MONOALG_PATH); fprintf(file,"save_pvd=true\n"); fprintf(file,"file_prefix=V\n"); fprintf(file,"save_activation_time=true\n"); - fprintf(file,"save_apd=true\n"); - fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libdefault_save_mesh_purkinje.so\n"); - fprintf(file,"main_function=save_purkinje_with_activation_times\n"); - fprintf(file,"init_function=init_save_purkinje_with_activation_times\n"); - fprintf(file,"end_function=end_save_purkinje_with_activation_times\n"); + fprintf(file,"activation_threshold=-50.0\n"); + fprintf(file,"save_apd=false\n"); + fprintf(file,"library_file=%s/shared_libs/libdefault_save_mesh_purkinje.so\n", MONOALG_PATH); + fprintf(file,"main_function=save_tissue_with_activation_times\n"); + fprintf(file,"init_function=init_save_tissue_with_activation_times\n"); + fprintf(file,"end_function=end_save_tissue_with_activation_times\n"); fprintf(file,"remove_older_simulation=true\n"); fprintf(file,"\n"); fprintf(file,"[assembly_matrix]\n"); fprintf(file,"init_function=set_initial_conditions_fvm\n"); - fprintf(file,"sigma_purkinje=%g\n",sigma); - fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libpurkinje_matrix_assembly.so\n"); - fprintf(file,"main_function=purkinje_fibers_assembly_matrix\n"); + fprintf(file,"sigma_x=%g\n",sigma); + fprintf(file,"sigma_y=%g\n",sigma); + fprintf(file,"sigma_z=%g\n",sigma); + fprintf(file,"library_file=%s/shared_libs/libdefault_matrix_assembly.so\n", MONOALG_PATH); + fprintf(file,"main_function=homogeneous_sigma_assembly_matrix\n"); fprintf(file,"\n"); fprintf(file,"[linear_system_solver]\n"); fprintf(file,"tolerance=1e-16\n"); fprintf(file,"use_preconditioner=yes\n"); fprintf(file,"max_iterations=500\n"); - fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libdefault_linear_system_solver.so\n"); + fprintf(file,"library_file=%s/shared_libs/libdefault_linear_system_solver.so\n", MONOALG_PATH); fprintf(file,"main_function=conjugate_gradient\n"); fprintf(file,"\n"); - fprintf(file,"[purkinje]\n"); - fprintf(file,"name=Simple Purkinje\n"); - //fprintf(file,"dx=100.0\n"); - fprintf(file,"dx=400.0\n"); - fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libdefault_purkinje.so\n"); - fprintf(file,"main_function=initialize_purkinje_with_custom_mesh\n"); - //fprintf(file,"network_file=/home/berg/Github/MonoAlg3D_C/networks/simple_cable_10cm.vtk\n"); - fprintf(file,"network_file=/home/berg/Github/MonoAlg3D_C/networks/simple_cable_40cm.vtk\n"); + fprintf(file,"[domain]\n"); + fprintf(file,"name=Simple Cable\n"); + fprintf(file,"start_dx=500.0\n"); + fprintf(file,"start_dy=500.0\n"); + fprintf(file,"start_dz=500.0\n"); + fprintf(file,"cable_length=20000.0\n"); + fprintf(file,"library_file=%s/shared_libs/libdefault_domains.so\n", MONOALG_PATH); + fprintf(file,"main_function=initialize_grid_with_cable_mesh\n"); fprintf(file,"\n"); fprintf(file,"[ode_solver]\n"); - //fprintf(file,"adaptive=true\n"); - fprintf(file,"dt=0.02\n"); + fprintf(file,"dt=0.01\n"); fprintf(file,"use_gpu=no\n"); fprintf(file,"gpu_id=0\n"); - //fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libtrovato_2019.so\n"); - //fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libToRORd_fkatp_endo.so\n"); - fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libten_tusscher_3_endo.so\n"); + fprintf(file,"library_file=%s/shared_libs/libToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.so\n", MONOALG_PATH); fprintf(file,"\n"); - fprintf(file,"[purkinje_stim_his]\n"); + fprintf(file,"[stim_benchmark]\n"); fprintf(file,"start = 0.0\n"); - fprintf(file,"duration = 5.0\n"); - //fprintf(file,"current = -40.0\n"); - fprintf(file,"current = -53.0\n"); - fprintf(file,"id_limit = 25\n"); - fprintf(file,"main_function=stim_if_id_less_than\n"); - fprintf(file,"library_file=/home/berg/Github/MonoAlg3D_C/shared_libs/libdefault_stimuli.so\n"); + fprintf(file,"duration = 2.0\n"); + fprintf(file,"current = -50.0\n"); + fprintf(file, "min_x = 0.0\n"); + fprintf(file, "max_x = 500.0\n"); + fprintf(file, "min_y = 0.0\n"); + fprintf(file, "max_y = 500.0\n"); + fprintf(file, "min_z = 0.0\n"); + fprintf(file, "max_z = 3000.0\n"); + fprintf(file,"main_function=stim_x_y_z_limits\n"); + fprintf(file,"library_file=%s/shared_libs/libdefault_stimuli.so\n", MONOALG_PATH); fprintf(file,"\n"); fclose(file); } -double calculate_error (const double cv, const double target_cv) -{ +double calculate_error (const double cv, const double target_cv) { return sqrt(pow(cv,2)-pow(target_cv,2))/target_cv; } -int main (int argc, char *argv[]) -{ - if (argc-1 != 1) - { +int main (int argc, char *argv[]) { + if (argc-1 != 1) { + cerr << "=============================================================================" << endl; cerr << "Usage:> " << argv[0] << " " << endl; cerr << "=============================================================================" << endl; @@ -175,12 +204,13 @@ int main (int argc, char *argv[]) double target_cv = atof(argv[1]); double sigma = 0.0002; - do - { + do { write_configuration_file(sigma); // Run the simulation - system("/home/berg/Github/MonoAlg3D_C/bin/MonoAlg3D -c /home/berg/Github/MonoAlg3D_C/scripts/tuneCV/configs/cable.ini"); + char command[500]; + sprintf(command,"%s/bin/MonoAlg3D -c %s/scripts/tuneCV/configs/cable.ini", MONOALG_PATH, MONOALG_PATH); + system(command); cv = calculate_conduction_velocity_from_simulation(); factor = pow(target_cv/cv,2); diff --git a/scripts/tuneCVbenchmark/.gitignore b/scripts/tuneCVbenchmark/.gitignore new file mode 100644 index 00000000..c7f2fd3a --- /dev/null +++ b/scripts/tuneCVbenchmark/.gitignore @@ -0,0 +1,5 @@ +/tables +build/* +outputs/* +bin/* +configs/* diff --git a/scripts/tuneCVbenchmark/CMakeLists.txt b/scripts/tuneCVbenchmark/CMakeLists.txt new file mode 100644 index 00000000..6928af75 --- /dev/null +++ b/scripts/tuneCVbenchmark/CMakeLists.txt @@ -0,0 +1,12 @@ +cmake_minimum_required(VERSION 2.8) + +PROJECT(tuneCV) + +find_package(VTK REQUIRED) +include(${VTK_USE_FILE}) + +set( CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/bin ) + +add_executable(tuneCV main.cpp ) + +target_link_libraries(tuneCV ${VTK_LIBRARIES}) diff --git a/scripts/tuneCVbenchmark/clean_project.sh b/scripts/tuneCVbenchmark/clean_project.sh new file mode 100755 index 00000000..d25ced98 --- /dev/null +++ b/scripts/tuneCVbenchmark/clean_project.sh @@ -0,0 +1,15 @@ +#!/bin/bash + +BUILD_DIR="build" + +if [[ "$#" -eq 1 ]]; then + BUILD_DIR=$1 +fi + +if [[ ! -d "${BUILD_DIR}" ]]; then + echo "Directory ${BUILD_DIR} does not exist" + exit +fi + +rm -fr ${BUILD_DIR}/* +rm -fr bin/* diff --git a/scripts/tuneCVbenchmark/main.cpp b/scripts/tuneCVbenchmark/main.cpp new file mode 100644 index 00000000..d558b11a --- /dev/null +++ b/scripts/tuneCVbenchmark/main.cpp @@ -0,0 +1,225 @@ +// Author: Lucas Berg + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +const double TOLERANCE = 1.0e-02; // 1 cm/s + +double calculate_conduction_velocity_from_cable_simulation () +{ + string filename = "outputs/cable/tissue_activation_time_map_pulse_it_0.vtu"; + + // Read all the data from the file + vtkSmartPointer reader = vtkSmartPointer::New(); + reader->SetFileName(filename.c_str()); + reader->Update(); + + vtkUnstructuredGrid* unstructuredGrid = reader->GetOutput(); + uint32_t num_points = unstructuredGrid->GetNumberOfPoints(); + uint32_t num_lines = unstructuredGrid->GetNumberOfCells(); + + vtkSmartPointer cellLocator = vtkSmartPointer::New(); + cellLocator->SetDataSet(unstructuredGrid); + cellLocator->BuildLocator(); + + // Points of interest + double point_x_0[3] = {4000, 250, 250}; + double point_x_1[3] = {18000, 250, 250}; + + // Find (closest points) Cell indexes for CV computation + vtkIdType cellId_x_0; // the cell id of the cell containing the closest point + vtkIdType cellId_x_1; // the cell id of the cell containing the closest point + + double closestPoint_x_0[3]; + double closestPoint_x_1[3]; + + int subId; // not needed + double closestPointDist2; // not needed + + cellLocator->FindClosestPoint(point_x_0, closestPoint_x_0, cellId_x_0, subId, closestPointDist2); + cellLocator->FindClosestPoint(point_x_1, closestPoint_x_1, cellId_x_1, subId, closestPointDist2); + + double delta_s_x = sqrt(pow(closestPoint_x_0[0]-closestPoint_x_1[0], 2)); // only changes over x-axis + + cout << delta_s_x << endl; + + // Read points scalar values + string array_name = "Scalars_"; + vtkSmartPointer array = vtkFloatArray::SafeDownCast(unstructuredGrid->GetCellData()->GetArray(array_name.c_str())); + + double cv_x = -1.0; + + if(array) + { + + double delta_lat_x = (array->GetValue(cellId_x_1) - array->GetValue(cellId_x_0)); // ms + + cout << delta_lat_x << endl; + cout << array->GetValue(cellId_x_0) << endl; + cout << array->GetValue(cellId_x_1) << endl; + + cv_x = (delta_s_x / delta_lat_x)*0.001; // {m/s} + } + else + { + cerr << "[!] ERROR! No Scalar_value found for the points!" << endl; + exit(EXIT_FAILURE); + } + + return cv_x; +} + +// TODO: Maybe pass a pre-configured config file as an input parameter with the cellular model setup that the user will use +void write_configuration_file (const double sigma) +{ + FILE *file = fopen("/home/jenny/MonoAlg3D_C/scripts/tuneCVbenchmark/configs/cable.ini","w+"); + + fprintf(file,"[main]\n"); + fprintf(file,"num_threads=6\n"); + fprintf(file,"dt_pde=0.01\n"); + fprintf(file,"simulation_time=100.0\n"); + fprintf(file,"abort_on_no_activity=false\n"); + fprintf(file,"use_adaptivity=false\n"); + fprintf(file,"quiet=true\n"); + fprintf(file,"\n"); + + fprintf(file,"[update_monodomain]\n"); + fprintf(file,"main_function=update_monodomain_default\n"); + fprintf(file,"library_file=/home/jenny/MonoAlg3D_C/shared_libs/libdefault_update_monodomain.so\n"); + fprintf(file,"\n"); + + // For saving the LATs in a format that can be read for calculating the CVs + fprintf(file,"[save_result]\n"); + fprintf(file,"print_rate=1\n"); + fprintf(file,"output_dir=/home/jenny/MonoAlg3D_C/scripts/tuneCVbenchmark/outputs/cable\n"); + fprintf(file,"save_pvd=true\n"); + fprintf(file,"file_prefix=V\n"); + fprintf(file,"save_activation_time=true\n"); + fprintf(file,"save_apd=false\n"); + fprintf(file,"library_file=/home/jenny/MonoAlg3D_C/shared_libs/libdefault_save_mesh_purkinje.so\n"); + fprintf(file,"main_function=save_tissue_with_activation_times\n"); + fprintf(file,"init_function=init_save_tissue_with_activation_times\n"); + fprintf(file,"end_function=end_save_tissue_with_activation_times\n"); + fprintf(file,"remove_older_simulation=true\n"); + fprintf(file,"\n"); + + fprintf(file,"[assembly_matrix]\n"); + fprintf(file,"init_function=set_initial_conditions_fvm\n"); + fprintf(file,"sigma_x=%g\n",sigma); + fprintf(file,"sigma_y=%g\n",sigma); + fprintf(file,"sigma_z=%g\n",sigma); + fprintf(file,"library_file=/home/jenny/MonoAlg3D_C/shared_libs/libdefault_matrix_assembly.so\n"); + fprintf(file,"main_function=homogeneous_sigma_assembly_matrix\n"); + fprintf(file,"\n"); + + fprintf(file,"[linear_system_solver]\n"); + fprintf(file,"tolerance=1e-15\n"); + fprintf(file,"use_preconditioner=no\n"); + fprintf(file,"use_gpu=yes\n"); + fprintf(file,"max_iterations=200\n"); + fprintf(file,"library_file=/home/jenny/MonoAlg3D_C/shared_libs/libdefault_linear_system_solver.so\n"); + fprintf(file,"init_function=init_conjugate_gradient\n"); + fprintf(file,"end_function=end_conjugate_gradient\n"); + fprintf(file,"main_function=conjugate_gradient\n"); + fprintf(file,"\n"); + + fprintf(file,"[domain]\n"); + fprintf(file,"name=Simple Cable\n"); + fprintf(file,"start_dx=500.0\n"); + fprintf(file,"start_dy=500.0\n"); + fprintf(file,"start_dz=500.0\n"); + fprintf(file,"cable_length=20000.0\n"); + fprintf(file,"library_file=/home/jenny/MonoAlg3D_C/shared_libs/libdefault_domains.so\n"); + fprintf(file,"main_function=initialize_grid_with_cable_mesh\n"); + fprintf(file,"\n"); + + fprintf(file,"[ode_solver]\n"); + fprintf(file,"dt=0.01\n"); + fprintf(file,"use_gpu=yes\n"); + fprintf(file,"gpu_id=0\n"); + fprintf(file,"library_file=/home/jenny/MonoAlg3D_C/shared_libs/libten_tusscher_tt3_mixed_endo_mid_epi.so\n"); + fprintf(file,"\n"); + + fprintf(file,"[stim_benchmark]\n"); + fprintf(file,"start = 0.0\n"); + fprintf(file,"duration = 2.0\n"); + fprintf(file,"current = -50.0\n"); + fprintf(file, "min_x = 0.0\n"); + fprintf(file, "max_x = 500.0\n"); + fprintf(file, "min_y = 0.0\n"); + fprintf(file, "max_y = 500.0\n"); + fprintf(file, "min_z = 0.0\n"); + fprintf(file, "max_z = 3000.0\n"); + fprintf(file,"main_function=stim_x_y_z_limits\n"); + fprintf(file,"library_file=/home/jenny/MonoAlg3D_C/shared_libs/libdefault_stimuli.so\n"); + fprintf(file,"\n"); + + fclose(file); +} + +int main (int argc, char *argv[]) +{ + if (argc-1 != 1) + { + cerr << "=============================================================================" << endl; + cerr << "Usage:> " << argv[0] << " " << endl; + cerr << "=============================================================================" << endl; + cerr << " = Target conduction velocity in m/s" << endl; + cerr << "=============================================================================" << endl; + cerr << "Example:" << endl; + cerr << argv[0] << " 0.67 (Longitudinal normal direction ventricle)" << endl; + cerr << argv[0] << " 0.33 (Transversal normal direction ventricle)" << endl; + cerr << argv[0] << " 0.17 (Sheet normal direction ventricle)" << endl; + cerr << argv[0] << " 1.90 (Purkinje fiber)" << endl; + cerr << "=============================================================================" << endl; + + exit(EXIT_FAILURE); + } + + double cv, factor; + double target_cv = atof(argv[1]); + double sigma = 0.0002; + + do + { + write_configuration_file(sigma); + + // Run the simulation + system("/home/jenny/MonoAlg3D_C/bin/MonoAlg3D -c /home/jenny/MonoAlg3D_C/scripts/tuneCVbenchmark/configs/cable.ini"); + + cv = calculate_conduction_velocity_from_cable_simulation(); + factor = pow(target_cv/cv,2); + sigma = sigma*factor; + + printf("\n|| Target CV = %g m/s || Computed CV = %g m/s || Factor = %g || Adjusted sigma = %g mS/um ||\n\n",target_cv,cv,factor,sigma); + + }while ( fabs(cv-target_cv) > TOLERANCE ); + + printf("\n[+] Target conductivity = %g mS/um\n",sigma); + + return 0; +} diff --git a/scripts/tuneCVbenchmark/recompile_project.sh b/scripts/tuneCVbenchmark/recompile_project.sh new file mode 100755 index 00000000..5a19dab7 --- /dev/null +++ b/scripts/tuneCVbenchmark/recompile_project.sh @@ -0,0 +1,21 @@ +#!/bin/bash + +BUILD_DIR="build" +BUILD_TYPE="Release" + +if [[ "$#" -eq 1 ]]; then + BUILD_DIR=$1 +fi + +if [[ "$#" -eq 2 ]]; then + BUILD_DIR=$1 + BUILD_TYPE=$2 +fi + + +if [[ ! -d "${BUILD_DIR}" ]]; then + echo "Directory ${BUILD_DIR} does not exist. Creating." + mkdir ${BUILD_DIR} +fi + +cd ${BUILD_DIR}; cmake -DCMAKE_BUILD_TYPE=${BUILD_TYPE} ..; make diff --git a/src/extra_data_library/helper_functions.c b/src/extra_data_library/helper_functions.c index 729490fb..b6377868 100644 --- a/src/extra_data_library/helper_functions.c +++ b/src/extra_data_library/helper_functions.c @@ -756,6 +756,447 @@ struct extra_data_for_torord_land * set_common_torord_Land_data (struct config * return extra_data; } +struct extra_data_for_torord_land_twave * set_common_torord_Land_twave_data (struct config *config, uint32_t num_cells) { + const uint32_t neq = 49; + struct extra_data_for_torord_land_twave *extra_data = MALLOC_ONE_TYPE(struct extra_data_for_torord_land_twave); + + 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); + extra_data->sf_IKs = MALLOC_ARRAY_OF_TYPE(real, num_cells); + extra_data->basetoapex = MALLOC_ARRAY_OF_TYPE(real, num_cells); + + return extra_data; +} + +struct extra_data_for_torord_gksgkrtjca_twave * set_common_torord_gksgkrtjca_twave_data (struct config *config, uint32_t num_cells) { + const uint32_t neq = 43; + struct extra_data_for_torord_gksgkrtjca_twave *extra_data = MALLOC_ONE_TYPE(struct extra_data_for_torord_gksgkrtjca_twave); + + 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); + + // Set the initial conditions (celltype = ENDO) + extra_data->initial_ss_endo[0] = -8.890585e+01; + extra_data->initial_ss_endo[1] = 1.107642e-02; + extra_data->initial_ss_endo[2] = 6.504164e-05; + extra_data->initial_ss_endo[3] = 1.210818e+01; + extra_data->initial_ss_endo[4] = 1.210851e+01; + extra_data->initial_ss_endo[5] = 1.426206e+02; + extra_data->initial_ss_endo[6] = 1.426205e+02; + extra_data->initial_ss_endo[7] = 1.530373e+00; + extra_data->initial_ss_endo[8] = 1.528032e+00; + extra_data->initial_ss_endo[9] = 7.455488e-05; + extra_data->initial_ss_endo[10] = 7.814592e-04; + extra_data->initial_ss_endo[11] = 8.313839e-01; + extra_data->initial_ss_endo[12] = 8.311938e-01; + extra_data->initial_ss_endo[13] = 6.752873e-01; + extra_data->initial_ss_endo[14] = 8.308255e-01; + extra_data->initial_ss_endo[15] = 1.585610e-04; + extra_data->initial_ss_endo[16] = 5.294475e-01; + extra_data->initial_ss_endo[17] = 2.896996e-01; + extra_data->initial_ss_endo[18] = 9.419166e-04; + extra_data->initial_ss_endo[19] = 9.996194e-01; + extra_data->initial_ss_endo[20] = 5.938602e-01; + extra_data->initial_ss_endo[21] = 4.799180e-04; + extra_data->initial_ss_endo[22] = 9.996194e-01; + extra_data->initial_ss_endo[23] = 6.543754e-01; + extra_data->initial_ss_endo[24] = -2.898677e-33; + extra_data->initial_ss_endo[25] = 1.000000e+00; + extra_data->initial_ss_endo[26] = 9.389659e-01; + extra_data->initial_ss_endo[27] = 1.000000e+00; + extra_data->initial_ss_endo[28] = 9.999003e-01; + extra_data->initial_ss_endo[29] = 9.999773e-01; + extra_data->initial_ss_endo[30] = 1.000000e+00; + extra_data->initial_ss_endo[31] = 1.000000e+00; + extra_data->initial_ss_endo[32] = 4.920606e-04; + extra_data->initial_ss_endo[33] = 8.337021e-04; + extra_data->initial_ss_endo[34] = 6.962775e-04; + extra_data->initial_ss_endo[35] = 8.425453e-04; + extra_data->initial_ss_endo[36] = 9.980807e-01; + extra_data->initial_ss_endo[37] = 1.289824e-05; + extra_data->initial_ss_endo[38] = 3.675442e-04; + extra_data->initial_ss_endo[39] = 2.471690e-01; + extra_data->initial_ss_endo[40] = 1.742987e-04; + extra_data->initial_ss_endo[41] = 5.421027e-24; + extra_data->initial_ss_endo[42] = 6.407933e-23; + + // Set the initial conditions (celltype = EPI) + extra_data->initial_ss_epi[0] = -8.917755e+01; + extra_data->initial_ss_epi[1] = 1.288116e-02; + extra_data->initial_ss_epi[2] = 5.767956e-05; + extra_data->initial_ss_epi[3] = 1.284260e+01; + extra_data->initial_ss_epi[4] = 1.284291e+01; + extra_data->initial_ss_epi[5] = 1.429114e+02; + extra_data->initial_ss_epi[6] = 1.429113e+02; + extra_data->initial_ss_epi[7] = 1.812268e+00; + extra_data->initial_ss_epi[8] = 1.810520e+00; + extra_data->initial_ss_epi[9] = 6.631866e-05; + extra_data->initial_ss_epi[10] = 7.370422e-04; + extra_data->initial_ss_epi[11] = 8.366816e-01; + extra_data->initial_ss_epi[12] = 8.366012e-01; + extra_data->initial_ss_epi[13] = 6.840260e-01; + extra_data->initial_ss_epi[14] = 8.363958e-01; + extra_data->initial_ss_epi[15] = 1.505860e-04; + extra_data->initial_ss_epi[16] = 5.412669e-01; + extra_data->initial_ss_epi[17] = 3.043382e-01; + extra_data->initial_ss_epi[18] = 9.248184e-04; + extra_data->initial_ss_epi[19] = 9.996371e-01; + extra_data->initial_ss_epi[20] = 9.996342e-01; + extra_data->initial_ss_epi[21] = 4.712023e-04; + extra_data->initial_ss_epi[22] = 9.996371e-01; + extra_data->initial_ss_epi[23] = 9.996366e-01; + extra_data->initial_ss_epi[24] = 4.333129e-43; + extra_data->initial_ss_epi[25] = 1.000000e+00; + extra_data->initial_ss_epi[26] = 9.485160e-01; + extra_data->initial_ss_epi[27] = 1.000000e+00; + extra_data->initial_ss_epi[28] = 9.999339e-01; + extra_data->initial_ss_epi[29] = 9.999822e-01; + extra_data->initial_ss_epi[30] = 1.000000e+00; + extra_data->initial_ss_epi[31] = 1.000000e+00; + extra_data->initial_ss_epi[32] = 3.086885e-04; + extra_data->initial_ss_epi[33] = 5.303737e-04; + extra_data->initial_ss_epi[34] = 6.775197e-04; + extra_data->initial_ss_epi[35] = 8.264829e-04; + extra_data->initial_ss_epi[36] = 9.982135e-01; + extra_data->initial_ss_epi[37] = 9.433146e-06; + extra_data->initial_ss_epi[38] = 2.730221e-04; + extra_data->initial_ss_epi[39] = 2.308784e-01; + extra_data->initial_ss_epi[40] = 1.690386e-04; + extra_data->initial_ss_epi[41] = -1.103286e-23; + extra_data->initial_ss_epi[42] = -6.177055e-22; + + // Set the initial conditions (celltype = MCELL) + extra_data->initial_ss_mid[0] = -8.924177e+01; + extra_data->initial_ss_mid[1] = 1.922391e-02; + extra_data->initial_ss_mid[2] = 6.585066e-05; + extra_data->initial_ss_mid[3] = 1.503347e+01; + extra_data->initial_ss_mid[4] = 1.503401e+01; + extra_data->initial_ss_mid[5] = 1.434407e+02; + extra_data->initial_ss_mid[6] = 1.434406e+02; + extra_data->initial_ss_mid[7] = 1.959747e+00; + extra_data->initial_ss_mid[8] = 1.963459e+00; + extra_data->initial_ss_mid[9] = 8.177438e-05; + extra_data->initial_ss_mid[10] = 7.269124e-04; + extra_data->initial_ss_mid[11] = 8.379059e-01; + extra_data->initial_ss_mid[12] = 8.377164e-01; + extra_data->initial_ss_mid[13] = 6.860578e-01; + extra_data->initial_ss_mid[14] = 8.372100e-01; + extra_data->initial_ss_mid[15] = 1.487602e-04; + extra_data->initial_ss_mid[16] = 5.350003e-01; + extra_data->initial_ss_mid[17] = 2.851164e-01; + extra_data->initial_ss_mid[18] = 9.208259e-04; + extra_data->initial_ss_mid[19] = 9.996411e-01; + extra_data->initial_ss_mid[20] = 5.673539e-01; + extra_data->initial_ss_mid[21] = 4.691672e-04; + extra_data->initial_ss_mid[22] = 9.996412e-01; + extra_data->initial_ss_mid[23] = 6.265825e-01; + extra_data->initial_ss_mid[24] = -4.922960e-40; + extra_data->initial_ss_mid[25] = 1.000000e+00; + extra_data->initial_ss_mid[26] = 9.200354e-01; + extra_data->initial_ss_mid[27] = 1.000000e+00; + extra_data->initial_ss_mid[28] = 9.997888e-01; + extra_data->initial_ss_mid[29] = 9.999665e-01; + extra_data->initial_ss_mid[30] = 1.000000e+00; + extra_data->initial_ss_mid[31] = 1.000000e+00; + extra_data->initial_ss_mid[32] = 5.161178e-04; + extra_data->initial_ss_mid[33] = 1.189422e-03; + extra_data->initial_ss_mid[34] = 6.917041e-04; + extra_data->initial_ss_mid[35] = 8.225453e-04; + extra_data->initial_ss_mid[36] = 9.979358e-01; + extra_data->initial_ss_mid[37] = 1.835276e-05; + extra_data->initial_ss_mid[38] = 5.316232e-04; + extra_data->initial_ss_mid[39] = 2.650323e-01; + extra_data->initial_ss_mid[40] = 1.678628e-04; + extra_data->initial_ss_mid[41] = 2.091039e-25; + extra_data->initial_ss_mid[42] = 2.438403e-23; + + extra_data->transmurality = MALLOC_ARRAY_OF_TYPE(real, num_cells); + extra_data->sf_IKs = 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 eca7df6c..d8936a38 100644 --- a/src/extra_data_library/helper_functions.h +++ b/src/extra_data_library/helper_functions.h @@ -156,7 +156,7 @@ struct extra_data_for_torord_land_twave { real *basetoapex; }; -struct extra_data_for_torord_general { +struct extra_data_for_torord_gksgkrtjca_twave { real INa_Multiplier; real INaL_Multiplier; real INaCa_Multiplier; @@ -182,7 +182,6 @@ struct extra_data_for_torord_general { real *initial_ss_mid; real *transmurality; real *sf_IKs; - real *basetoapex; }; struct extra_data_for_fibrosis * set_common_schemia_data(struct config *config, uint32_t num_cells); @@ -191,14 +190,7 @@ struct extra_data_for_torord * set_common_torord_data (struct config *config, ui 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_torord_land_twave * set_common_torord_Land_twave_data (struct config *config, uint32_t num_cells); +struct extra_data_for_torord_gksgkrtjca_twave * set_common_torord_gksgkrtjca_twave_data (struct config *config, uint32_t num_cells); struct extra_data_for_trovato * set_common_trovato_data (struct config *config, uint32_t num_cells); -struct extra_data_for_torord_general * set_common_torord_general (struct config *config, uint32_t num_cells, uint32_t model_id); - -void configure_torord_extra_data_with_user_input (struct extra_data_for_torord_general *extra_data, struct config *config, uint32_t model_id); -void configure_torord_extra_data_initial_conditions (struct extra_data_for_torord_general *extra_data, uint32_t model_id); -void configure_torord_fkatp_initial_conditions (struct extra_data_for_torord_general *extra_data); -void configure_torord_dynCl_initial_conditions (struct extra_data_for_torord_general *extra_data); -void configure_torord_Land_initial_conditions (struct extra_data_for_torord_general *extra_data); - #endif // MONOALG3D_C_EXTRA_DATA_HELPER_FUNCTIONS_H diff --git a/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.c b/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.c new file mode 100644 index 00000000..ba2eafdd --- /dev/null +++ b/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.c @@ -0,0 +1,537 @@ +#include "ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.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_fkatp_2019 CPU model with GKs GKr and tjca adjustments\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; + real *sf_Iks = NULL; + if(solver->ode_extra_data) { + struct extra_data_for_torord_gksgkrtjca_twave *extra_data = (struct extra_data_for_torord_gksgkrtjca_twave*)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; + sf_Iks = extra_data->sf_IKs; + + 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]; + + // Steady-state after 200 beats (endocardium cell) + sv[0] = -8.890585e+01; + sv[1] = 1.107642e-02; + sv[2] = 6.504164e-05; + sv[3] = 1.210818e+01; + sv[4] = 1.210851e+01; + sv[5] = 1.426206e+02; + sv[6] = 1.426205e+02; + sv[7] = 1.530373e+00; + sv[8] = 1.528032e+00; + sv[9] = 7.455488e-05; + sv[10] = 7.814592e-04; + sv[11] = 8.313839e-01; + sv[12] = 8.311938e-01; + sv[13] = 6.752873e-01; + sv[14] = 8.308255e-01; + sv[15] = 1.585610e-04; + sv[16] = 5.294475e-01; + sv[17] = 2.896996e-01; + sv[18] = 9.419166e-04; + sv[19] = 9.996194e-01; + sv[20] = 5.938602e-01; + sv[21] = 4.799180e-04; + sv[22] = 9.996194e-01; + sv[23] = 6.543754e-01; + sv[24] = -2.898677e-33; + sv[25] = 1.000000e+00; + sv[26] = 9.389659e-01; + sv[27] = 1.000000e+00; + sv[28] = 9.999003e-01; + sv[29] = 9.999773e-01; + sv[30] = 1.000000e+00; + sv[31] = 1.000000e+00; + sv[32] = 4.920606e-04; + sv[33] = 8.337021e-04; + sv[34] = 6.962775e-04; + sv[35] = 8.425453e-04; + sv[36] = 9.980807e-01; + sv[37] = 1.289824e-05; + sv[38] = 3.675442e-04; + sv[39] = 2.471690e-01; + sv[40] = 1.742987e-04; + sv[41] = 5.421027e-24; + sv[42] = 6.407933e-23; + } + } +} + +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; + real *sf_Iks = NULL; + if (ode_solver->ode_extra_data) { + struct extra_data_for_torord_gksgkrtjca_twave *extra_data = (struct extra_data_for_torord_gksgkrtjca_twave*)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; + sf_Iks = extra_data->sf_IKs; + 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; + } + } + + 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], sf_Iks[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, 1.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], sf_Iks[i], extra_par); + } + else { + solve_model_ode_cpu(dt, sv + (sv_id * NEQ), stim_currents[i], 0.0, 1.0, extra_par); + } + } + } + } +} + +void solve_model_ode_cpu(real dt, real *sv, real stim_current, real transmurality, real sf_Iks, 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, sf_Iks, 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); // CaMKt + SOLVE_EQUATION_EULER_CPU(2); // cass + SOLVE_EQUATION_EULER_CPU(3); // nai + SOLVE_EQUATION_EULER_CPU(4); // nass + SOLVE_EQUATION_EULER_CPU(5); // ki + SOLVE_EQUATION_EULER_CPU(6); // kss + SOLVE_EQUATION_EULER_CPU(7); // cansr + SOLVE_EQUATION_EULER_CPU(8); // cajsr + SOLVE_EQUATION_EULER_CPU(9); // cai + SOLVE_EQUATION_RUSH_LARSEN_CPU(10); // m + SOLVE_EQUATION_RUSH_LARSEN_CPU(11); // h + SOLVE_EQUATION_RUSH_LARSEN_CPU(12); // j + SOLVE_EQUATION_RUSH_LARSEN_CPU(13); // hp + SOLVE_EQUATION_RUSH_LARSEN_CPU(14); // jp + SOLVE_EQUATION_RUSH_LARSEN_CPU(15); // mL + SOLVE_EQUATION_RUSH_LARSEN_CPU(16); // hL + SOLVE_EQUATION_RUSH_LARSEN_CPU(17); // hLp + SOLVE_EQUATION_RUSH_LARSEN_CPU(18); // a + SOLVE_EQUATION_RUSH_LARSEN_CPU(19); // iF + SOLVE_EQUATION_RUSH_LARSEN_CPU(20); // iS + SOLVE_EQUATION_RUSH_LARSEN_CPU(21); // ap + SOLVE_EQUATION_RUSH_LARSEN_CPU(22); // iFp + SOLVE_EQUATION_RUSH_LARSEN_CPU(23); // iSp + SOLVE_EQUATION_RUSH_LARSEN_CPU(24); // d + SOLVE_EQUATION_RUSH_LARSEN_CPU(25); // ff + SOLVE_EQUATION_RUSH_LARSEN_CPU(26); // fs + SOLVE_EQUATION_RUSH_LARSEN_CPU(27); // fcaf + SOLVE_EQUATION_RUSH_LARSEN_CPU(28); // fcas + SOLVE_EQUATION_RUSH_LARSEN_CPU(29); // jca + SOLVE_EQUATION_RUSH_LARSEN_CPU(30); // ffp + SOLVE_EQUATION_RUSH_LARSEN_CPU(31); // fcafp + SOLVE_EQUATION_EULER_CPU(32); // nca + SOLVE_EQUATION_EULER_CPU(33); // nca_i + SOLVE_EQUATION_EULER_CPU(34); // ikr_c0 + SOLVE_EQUATION_EULER_CPU(35); // ikr_c1 + SOLVE_EQUATION_EULER_CPU(36); // ikr_c2 + SOLVE_EQUATION_EULER_CPU(37); // ikr_i + SOLVE_EQUATION_EULER_CPU(38); // ikr_o + SOLVE_EQUATION_RUSH_LARSEN_CPU(39); // xs1 + SOLVE_EQUATION_RUSH_LARSEN_CPU(40); // xs2 + SOLVE_EQUATION_RUSH_LARSEN_CPU(41); // Jrel_np + SOLVE_EQUATION_RUSH_LARSEN_CPU(42); // Jrel_p +} + +void solve_forward_euler_cpu_adpt(real *sv, real stim_curr, real transmurality, real sf_Iks, 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, sf_Iks, 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,sf_Iks, 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 sf_Iks, 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 + real v = sv[0]; + real CaMKt = sv[1]; + real cass = sv[2]; + real nai = sv[3]; + real nass = sv[4]; + real ki = sv[5]; + real kss = sv[6]; + real cansr = sv[7]; + real cajsr = sv[8]; + real cai = sv[9]; + real m = sv[10]; + real h = sv[11]; + real j = sv[12]; + real hp = sv[13]; + real jp = sv[14]; + real mL = sv[15]; + real hL = sv[16]; + real hLp = sv[17]; + real a = sv[18]; + real iF = sv[19]; + real iS = sv[20]; + real ap = sv[21]; + real iFp = sv[22]; + real iSp = sv[23]; + real d = sv[24]; + real ff = sv[25]; + real fs = sv[26]; + real fcaf = sv[27]; + real fcas = sv[28]; + real jca = sv[29]; + real ffp = sv[30]; + real fcafp = sv[31]; + real nca = sv[32]; + real nca_i = sv[33]; + real ikr_c0 = sv[34]; + real ikr_c1 = sv[35]; + real ikr_c2 = sv[36]; + real ikr_i = sv[37]; + real ikr_o = sv[38]; + real xs1 = sv[39]; + real xs2 = sv[40]; + real Jrel_np = sv[41]; + real Jrel_p = sv[42]; + + #include "ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.common.c" +} + +void RHS_RL_cpu(real *a_, real *b_, const real *sv, real *rDY_, real stim_current, real dt, real transmurality, real sf_Iks, 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 + real v = sv[0]; + real CaMKt = sv[1]; + real cass = sv[2]; + real nai = sv[3]; + real nass = sv[4]; + real ki = sv[5]; + real kss = sv[6]; + real cansr = sv[7]; + real cajsr = sv[8]; + real cai = sv[9]; + real m = sv[10]; + real h = sv[11]; + real j = sv[12]; + real hp = sv[13]; + real jp = sv[14]; + real mL = sv[15]; + real hL = sv[16]; + real hLp = sv[17]; + real a = sv[18]; + real iF = sv[19]; + real iS = sv[20]; + real ap = sv[21]; + real iFp = sv[22]; + real iSp = sv[23]; + real d = sv[24]; + real ff = sv[25]; + real fs = sv[26]; + real fcaf = sv[27]; + real fcas = sv[28]; + real jca = sv[29]; + real ffp = sv[30]; + real fcafp = sv[31]; + real nca = sv[32]; + real nca_i = sv[33]; + real ikr_c0 = sv[34]; + real ikr_c1 = sv[35]; + real ikr_c2 = sv[36]; + real ikr_i = sv[37]; + real ikr_o = sv[38]; + real xs1 = sv[39]; + real xs2 = sv[40]; + real Jrel_np = sv[41]; + real Jrel_p = sv[42]; + + #include "ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments_RL.common.c" +} diff --git a/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.common.c b/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.common.c new file mode 100644 index 00000000..f20b04be --- /dev/null +++ b/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.common.c @@ -0,0 +1,678 @@ +// 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; + +// 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 cli = 24; // Intracellular Cl [mM] +real clo = 150; // Extracellular Cl [mM] + +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)); + +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))); + +// convenient shorthand calculations +real vffrt=v*F*F/(R*T); +real vfrt=v*F/(R*T); +real frt = F/(R*T); + +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 = 60; // 75 original. Adjustments to enable GKs to play a bigger role in APD modulation. T wave personalisation. +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 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 +// 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; + +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 +real c0 = ikr_c0; +real c1 = ikr_c1; +real c2 = ikr_c2; +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 * 0.6; // 1st element compensates for change to ko (sqrt(5/5.4)* 0.0362) Adjustment for T wave personalisation. +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 * 5.0 * IKs_Multiplier; +real GKs= 0.0011 * 5.0 * sf_Iks; // Adjustment for T wave personalisation +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); + +// 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 * INaK_Multiplier; +if (celltype==EPI) + Pnak=Pnak*0.9; +else if (celltype==2) + Pnak=Pnak*0.7; + +real INaK=Pnak*(zna*JnakNa+zk*JnakK); + +// Minor/background currents +// calculate IKb +real xkb=1.0/(1.0+exp(-(v-10.8968)/(23.9871))); +real GKb=0.0189*IKb_Multiplier; +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 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-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 +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+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; + +// 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)); +real dcai=Bcai*(-(ICaL_i + IpCa+ICab-2.0*INaCa_i)*Acap/(2.0*F*vmyo)-Jup*vnsr/vmyo+Jdiff*vss/vmyo); + +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); + +// Right-hand side +rDY_[0] = dv; +rDY_[1] = dCaMKt; +rDY_[2] = dcass; +rDY_[3] = dnai; +rDY_[4] = dnass; +rDY_[5] = dki; +rDY_[6] = dkss; +rDY_[7] = dcansr; +rDY_[8] = dcajsr; +rDY_[9] = dcai; +rDY_[10] = dm; +rDY_[11] = dh; +rDY_[12] = dj; +rDY_[13] = dhp; +rDY_[14] = djp; +rDY_[15] = dmL; +rDY_[16] = dhL; +rDY_[17] = dhLp; +rDY_[18] = da; +rDY_[19] = diF; +rDY_[20] = diS; +rDY_[21] = dap; +rDY_[22] = diFp; +rDY_[23] = diSp; +rDY_[24] = dd; +rDY_[25] = dff; +rDY_[26] = dfs; +rDY_[27] = dfcaf; +rDY_[28] = dfcas; +rDY_[29] = djca; +rDY_[30] = dffp; +rDY_[31] = dfcafp; +rDY_[32] = dnca; +rDY_[33] = dnca_i; +rDY_[34] = dc0; +rDY_[35] = dc1; +rDY_[36] = dc2; +rDY_[37] = di; +rDY_[38] = delta_o; +rDY_[39] = dxs1; +rDY_[40] = dxs2; +rDY_[41] = dJrelnp; +rDY_[42] = dJrelp; diff --git a/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.cu b/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.cu new file mode 100644 index 00000000..09d7b899 --- /dev/null +++ b/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.cu @@ -0,0 +1,812 @@ +#include "ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.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; + + // Default initial conditions (endocardium cell) + if (threadID < num_volumes) { + for (int i = 0; i < NEQ; i++) { + // Steady-state after 200 beats (endocardium cell) + *((real * )((char *) sv + pitch * 0) + threadID) = -8.890585e+01; + *((real * )((char *) sv + pitch * 1) + threadID) = 1.107642e-02; + *((real * )((char *) sv + pitch * 2) + threadID) = 6.504164e-05; + *((real * )((char *) sv + pitch * 3) + threadID) = 1.210818e+01; + *((real * )((char *) sv + pitch * 4) + threadID) = 1.210851e+01; + *((real * )((char *) sv + pitch * 5) + threadID) = 1.426206e+02; + *((real * )((char *) sv + pitch * 6) + threadID) = 1.426205e+02; + *((real * )((char *) sv + pitch * 7) + threadID) = 1.530373e+00; + *((real * )((char *) sv + pitch * 8) + threadID) = 1.528032e+00; + *((real * )((char *) sv + pitch * 9) + threadID) = 7.455488e-05; + *((real * )((char *) sv + pitch * 10) + threadID) = 7.814592e-04; + *((real * )((char *) sv + pitch * 11) + threadID) = 8.313839e-01; + *((real * )((char *) sv + pitch * 12) + threadID) = 8.311938e-01; + *((real * )((char *) sv + pitch * 13) + threadID) = 6.752873e-01; + *((real * )((char *) sv + pitch * 14) + threadID) = 8.308255e-01; + *((real * )((char *) sv + pitch * 15) + threadID) = 1.585610e-04; + *((real * )((char *) sv + pitch * 16) + threadID) = 5.294475e-01; + *((real * )((char *) sv + pitch * 17) + threadID) = 2.896996e-01; + *((real * )((char *) sv + pitch * 18) + threadID) = 9.419166e-04; + *((real * )((char *) sv + pitch * 19) + threadID) = 9.996194e-01; + *((real * )((char *) sv + pitch * 20) + threadID) = 5.938602e-01; + *((real * )((char *) sv + pitch * 21) + threadID) = 4.799180e-04; + *((real * )((char *) sv + pitch * 22) + threadID) = 9.996194e-01; + *((real * )((char *) sv + pitch * 23) + threadID) = 6.543754e-01; + *((real * )((char *) sv + pitch * 24) + threadID) = -2.898677e-33; + *((real * )((char *) sv + pitch * 25) + threadID) = 1.000000e+00; + *((real * )((char *) sv + pitch * 26) + threadID) = 9.389659e-01; + *((real * )((char *) sv + pitch * 27) + threadID) = 1.000000e+00; + *((real * )((char *) sv + pitch * 28) + threadID) = 9.999003e-01; + *((real * )((char *) sv + pitch * 29) + threadID) = 9.999773e-01; + *((real * )((char *) sv + pitch * 30) + threadID) = 1.000000e+00; + *((real * )((char *) sv + pitch * 31) + threadID) = 1.000000e+00; + *((real * )((char *) sv + pitch * 32) + threadID) = 4.920606e-04; + *((real * )((char *) sv + pitch * 33) + threadID) = 8.337021e-04; + *((real * )((char *) sv + pitch * 34) + threadID) = 6.962775e-04; + *((real * )((char *) sv + pitch * 35) + threadID) = 8.425453e-04; + *((real * )((char *) sv + pitch * 36) + threadID) = 9.980807e-01; + *((real * )((char *) sv + pitch * 37) + threadID) = 1.289824e-05; + *((real * )((char *) sv + pitch * 38) + threadID) = 3.675442e-04; + *((real * )((char *) sv + pitch * 39) + threadID) = 2.471690e-01; + *((real * )((char *) sv + pitch * 40) + threadID) = 1.742987e-04; + *((real * )((char *) sv + pitch * 41) + threadID) = 5.421027e-24; + *((real * )((char *) sv + pitch * 42) + threadID) = 6.407933e-23; + } + + if(use_adpt_dt) { + *((real *)((char *)sv + pitch * 43) + threadID) = min_dt; // dt + *((real *)((char *)sv + pitch * 44) + threadID) = 0.0; // time_new + *((real *)((char *)sv + pitch * 45) + 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, real *sf_Iks) { + 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 * 43) + threadID) = min_dt; // dt + *((real *)((char *)sv + pitch * 44) + threadID) = 0.0; // time_new + *((real *)((char *)sv + pitch * 45) + 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 *sf_Iks = NULL; + real *initial_conditions_endo_device = NULL; + real *initial_conditions_epi_device = NULL; + real *initial_conditions_mid_device = NULL; + real *transmurality_device = NULL; + real *sf_Iks_device = NULL; + + if(solver->ode_extra_data) { + struct extra_data_for_torord_gksgkrtjca_twave *extra_data = (struct extra_data_for_torord_gksgkrtjca_twave*)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; + sf_Iks = extra_data->sf_IKs; + 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)); + check_cuda_error(cudaMalloc((void **)&sf_Iks_device, sizeof(real)*num_volumes)); + check_cuda_error(cudaMemcpy(sf_Iks_device, sf_Iks, 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, sf_Iks_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; + bool has_extra_params = (ode_solver->ode_extra_data) ? true : false; + + // 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)); + } + + uint32_t num_volumes = ode_solver->original_num_cells; + real *transmurality = NULL; + real *transmurality_device = NULL; + real *sf_Iks = NULL; + real *sf_Iks_device = NULL; + int num_extra_parameters = 20; + real extra_par[num_extra_parameters]; + real *extra_par_device = NULL; + + // Get the extra data array if exists + // Transmurality and sf_IKs mapping defined on 'extra_data' function + if(ode_solver->ode_extra_data) { + struct extra_data_for_torord_gksgkrtjca_twave *extra_data = (struct extra_data_for_torord_gksgkrtjca_twave*)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; + sf_Iks = extra_data->sf_IKs; + 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 **)&sf_Iks_device, sizeof(real)*num_volumes)); + check_cuda_error(cudaMemcpy(sf_Iks_device, sf_Iks, 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)); + } + // No [extra_data] section, we consider all cells ENDO! + 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)); + } + + // Call the kernel function to solve the cellular model on the GPU + solve_endo_mid_epi_gpu<<>>(current_t, dt, sv, stims_currents_device, cells_to_solve_device, transmurality_device, sf_Iks_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, has_extra_params); + + 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 (sf_Iks_device) check_cuda_error(cudaFree(sf_Iks_device)); + if (extra_par_device) check_cuda_error(cudaFree(extra_par_device)); +} + +__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 *sf_Iks, 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, bool has_extra_params) { + 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) { + + // [!] Transmurality and sf_IKs are given by the columns of ALG mesh file + if (has_extra_params) { + RHS_RL_gpu(a, b, sv, rDY, stim_currents[threadID], transmurality[threadID], sf_Iks[threadID], extra_params, sv_id, dt, pitch, false); + } + // Default: all cells are ENDO and sf_IKs equals to 1 + else { + RHS_RL_gpu(a, b, sv, rDY, stim_currents[threadID], 0.0, 1.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); // CaMKt + SOLVE_EQUATION_EULER_GPU(2); // cass + SOLVE_EQUATION_EULER_GPU(3); // nai + SOLVE_EQUATION_EULER_GPU(4); // nass + SOLVE_EQUATION_EULER_GPU(5); // ki + SOLVE_EQUATION_EULER_GPU(6); // kss + SOLVE_EQUATION_EULER_GPU(7); // cansr + SOLVE_EQUATION_EULER_GPU(8); // cajsr + SOLVE_EQUATION_EULER_GPU(9); // cai + SOLVE_EQUATION_RUSH_LARSEN_GPU(10); // m + SOLVE_EQUATION_RUSH_LARSEN_GPU(11); // h + SOLVE_EQUATION_RUSH_LARSEN_GPU(12); // j + SOLVE_EQUATION_RUSH_LARSEN_GPU(13); // hp + SOLVE_EQUATION_RUSH_LARSEN_GPU(14); // jp + SOLVE_EQUATION_RUSH_LARSEN_GPU(15); // mL + SOLVE_EQUATION_RUSH_LARSEN_GPU(16); // hL + SOLVE_EQUATION_RUSH_LARSEN_GPU(17); // hLp + SOLVE_EQUATION_RUSH_LARSEN_GPU(18); // a + SOLVE_EQUATION_RUSH_LARSEN_GPU(19); // iF + SOLVE_EQUATION_RUSH_LARSEN_GPU(20); // iS + SOLVE_EQUATION_RUSH_LARSEN_GPU(21); // ap + SOLVE_EQUATION_RUSH_LARSEN_GPU(22); // iFp + SOLVE_EQUATION_RUSH_LARSEN_GPU(23); // iSp + SOLVE_EQUATION_RUSH_LARSEN_GPU(24); // d + SOLVE_EQUATION_RUSH_LARSEN_GPU(25); // ff + SOLVE_EQUATION_RUSH_LARSEN_GPU(26); // fs + SOLVE_EQUATION_RUSH_LARSEN_GPU(27); // fcaf + SOLVE_EQUATION_RUSH_LARSEN_GPU(28); // fcas + SOLVE_EQUATION_RUSH_LARSEN_GPU(29); // jca + SOLVE_EQUATION_RUSH_LARSEN_GPU(30); // ffp + SOLVE_EQUATION_RUSH_LARSEN_GPU(31); // fcafp + SOLVE_EQUATION_EULER_GPU(32); // nca + SOLVE_EQUATION_EULER_GPU(33); // nca_i + SOLVE_EQUATION_EULER_GPU(34); // ikr_c0 + SOLVE_EQUATION_EULER_GPU(35); // ikr_c1 + SOLVE_EQUATION_EULER_GPU(36); // ikr_c2 + SOLVE_EQUATION_EULER_GPU(37); // ikr_i + SOLVE_EQUATION_EULER_GPU(38); // ikr_o + SOLVE_EQUATION_RUSH_LARSEN_GPU(39); // xs1 + SOLVE_EQUATION_RUSH_LARSEN_GPU(40); // xs2 + SOLVE_EQUATION_RUSH_LARSEN_GPU(41); // Jrel_np + SOLVE_EQUATION_RUSH_LARSEN_GPU(42); // Jrel_p + } + } else { + solve_forward_euler_gpu_adpt(sv, stim_currents[threadID], transmurality[threadID], sf_Iks[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 transmurality, real sf_Iks, 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, transmurality, sf_Iks, 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, transmurality, sf_Iks, 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 transmurality, real sf_Iks, 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 = transmurality; + + // Get the stimulus current from the current cell + real calc_I_stim = stim_current; + + // State variables + real v; + real CaMKt; + real cass; + real nai; + real nass; + real ki; + real kss; + real cansr; + real cajsr; + real cai; + real m; + real h; + real j; + real hp; + real jp; + real mL; + real hL; + real hLp; + real a; + real iF; + real iS; + real ap; + real iFp; + real iSp; + real d; + real ff; + real fs; + real fcaf; + real fcas; + real jca; + real ffp; + real fcafp; + real nca; + real nca_i; + real ikr_c0; + real ikr_c1; + real ikr_c2; + real ikr_i; + real ikr_o; + real xs1; + real xs2; + real Jrel_np; + real Jrel_p; + + if (use_adpt_dt) { + v = sv[0]; + CaMKt = sv[1]; + cass = sv[2]; + nai = sv[3]; + nass = sv[4]; + ki = sv[5]; + kss = sv[6]; + cansr = sv[7]; + cajsr = sv[8]; + cai = sv[9]; + m = sv[10]; + h = sv[11]; + j = sv[12]; + hp = sv[13]; + jp = sv[14]; + mL = sv[15]; + hL = sv[16]; + hLp = sv[17]; + a = sv[18]; + iF = sv[19]; + iS = sv[20]; + ap = sv[21]; + iFp = sv[22]; + iSp = sv[23]; + d = sv[24]; + ff = sv[25]; + fs = sv[26]; + fcaf = sv[27]; + fcas = sv[28]; + jca = sv[29]; + ffp = sv[30]; + fcafp = sv[31]; + nca = sv[32]; + nca_i = sv[33]; + ikr_c0 = sv[34]; + ikr_c1 = sv[35]; + ikr_c2 = sv[36]; + ikr_i = sv[37]; + ikr_o = sv[38]; + xs1 = sv[39]; + xs2 = sv[40]; + Jrel_np = sv[41]; + Jrel_p = sv[42]; + } else { + v = *((real *)((char *)sv + pitch * 0) + threadID_); + CaMKt = *((real *)((char *)sv + pitch * 1) + threadID_); + cass = *((real *)((char *)sv + pitch * 2) + threadID_); + nai = *((real *)((char *)sv + pitch * 3) + threadID_); + nass = *((real *)((char *)sv + pitch * 4) + threadID_); + ki = *((real *)((char *)sv + pitch * 5) + threadID_); + kss = *((real *)((char *)sv + pitch * 6) + threadID_); + cansr = *((real *)((char *)sv + pitch * 7) + threadID_); + cajsr = *((real *)((char *)sv + pitch * 8) + threadID_); + cai = *((real *)((char *)sv + pitch * 9) + threadID_); + m = *((real *)((char *)sv + pitch * 10) + threadID_); + h = *((real *)((char *)sv + pitch * 11) + threadID_); + j = *((real *)((char *)sv + pitch * 12) + threadID_); + hp = *((real *)((char *)sv + pitch * 13) + threadID_); + jp = *((real *)((char *)sv + pitch * 14) + threadID_); + mL = *((real *)((char *)sv + pitch * 15) + threadID_); + hL = *((real *)((char *)sv + pitch * 16) + threadID_); + hLp = *((real *)((char *)sv + pitch * 17) + threadID_); + a = *((real *)((char *)sv + pitch * 18) + threadID_); + iF = *((real *)((char *)sv + pitch * 19) + threadID_); + iS = *((real *)((char *)sv + pitch * 20) + threadID_); + ap = *((real *)((char *)sv + pitch * 21) + threadID_); + iFp = *((real *)((char *)sv + pitch * 22) + threadID_); + iSp = *((real *)((char *)sv + pitch * 23) + threadID_); + d = *((real *)((char *)sv + pitch * 24) + threadID_); + ff = *((real *)((char *)sv + pitch * 25) + threadID_); + fs = *((real *)((char *)sv + pitch * 26) + threadID_); + fcaf = *((real *)((char *)sv + pitch * 27) + threadID_); + fcas = *((real *)((char *)sv + pitch * 28) + threadID_); + jca = *((real *)((char *)sv + pitch * 29) + threadID_); + ffp = *((real *)((char *)sv + pitch * 30) + threadID_); + fcafp = *((real *)((char *)sv + pitch * 31) + threadID_); + nca = *((real *)((char *)sv + pitch * 32) + threadID_); + nca_i = *((real *)((char *)sv + pitch * 33) + threadID_); + ikr_c0 = *((real *)((char *)sv + pitch * 34) + threadID_); + ikr_c1 = *((real *)((char *)sv + pitch * 35) + threadID_); + ikr_c2 = *((real *)((char *)sv + pitch * 36) + threadID_); + ikr_i = *((real *)((char *)sv + pitch * 37) + threadID_); + ikr_o = *((real *)((char *)sv + pitch * 38) + threadID_); + xs1 = *((real *)((char *)sv + pitch * 39) + threadID_); + xs2 = *((real *)((char *)sv + pitch * 40) + threadID_); + Jrel_np = *((real *)((char *)sv + pitch * 41) + threadID_); + Jrel_p = *((real *)((char *)sv + pitch * 42) + threadID_); + } + + #include "ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.common.c" +} + + +inline __device__ void RHS_RL_gpu(real *a_, real *b_, real *sv, real *rDY_, real stim_current, real transmurality, real sf_Iks, 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 = transmurality; + + // Get the stimulus current from the current cell + real calc_I_stim = stim_current; + + // State variables + real v; + real CaMKt; + real cass; + real nai; + real nass; + real ki; + real kss; + real cansr; + real cajsr; + real cai; + real m; + real h; + real j; + real hp; + real jp; + real mL; + real hL; + real hLp; + real a; + real iF; + real iS; + real ap; + real iFp; + real iSp; + real d; + real ff; + real fs; + real fcaf; + real fcas; + real jca; + real ffp; + real fcafp; + real nca; + real nca_i; + real ikr_c0; + real ikr_c1; + real ikr_c2; + real ikr_i; + real ikr_o; + real xs1; + real xs2; + real Jrel_np; + real Jrel_p; + + if (use_adpt_dt) { + v = sv[0]; + CaMKt = sv[1]; + cass = sv[2]; + nai = sv[3]; + nass = sv[4]; + ki = sv[5]; + kss = sv[6]; + cansr = sv[7]; + cajsr = sv[8]; + cai = sv[9]; + m = sv[10]; + h = sv[11]; + j = sv[12]; + hp = sv[13]; + jp = sv[14]; + mL = sv[15]; + hL = sv[16]; + hLp = sv[17]; + a = sv[18]; + iF = sv[19]; + iS = sv[20]; + ap = sv[21]; + iFp = sv[22]; + iSp = sv[23]; + d = sv[24]; + ff = sv[25]; + fs = sv[26]; + fcaf = sv[27]; + fcas = sv[28]; + jca = sv[29]; + ffp = sv[30]; + fcafp = sv[31]; + nca = sv[32]; + nca_i = sv[33]; + ikr_c0 = sv[34]; + ikr_c1 = sv[35]; + ikr_c2 = sv[36]; + ikr_i = sv[37]; + ikr_o = sv[38]; + xs1 = sv[39]; + xs2 = sv[40]; + Jrel_np = sv[41]; + Jrel_p = sv[42]; + } else { + v = *((real *)((char *)sv + pitch * 0) + threadID_); + CaMKt = *((real *)((char *)sv + pitch * 1) + threadID_); + cass = *((real *)((char *)sv + pitch * 2) + threadID_); + nai = *((real *)((char *)sv + pitch * 3) + threadID_); + nass = *((real *)((char *)sv + pitch * 4) + threadID_); + ki = *((real *)((char *)sv + pitch * 5) + threadID_); + kss = *((real *)((char *)sv + pitch * 6) + threadID_); + cansr = *((real *)((char *)sv + pitch * 7) + threadID_); + cajsr = *((real *)((char *)sv + pitch * 8) + threadID_); + cai = *((real *)((char *)sv + pitch * 9) + threadID_); + m = *((real *)((char *)sv + pitch * 10) + threadID_); + h = *((real *)((char *)sv + pitch * 11) + threadID_); + j = *((real *)((char *)sv + pitch * 12) + threadID_); + hp = *((real *)((char *)sv + pitch * 13) + threadID_); + jp = *((real *)((char *)sv + pitch * 14) + threadID_); + mL = *((real *)((char *)sv + pitch * 15) + threadID_); + hL = *((real *)((char *)sv + pitch * 16) + threadID_); + hLp = *((real *)((char *)sv + pitch * 17) + threadID_); + a = *((real *)((char *)sv + pitch * 18) + threadID_); + iF = *((real *)((char *)sv + pitch * 19) + threadID_); + iS = *((real *)((char *)sv + pitch * 20) + threadID_); + ap = *((real *)((char *)sv + pitch * 21) + threadID_); + iFp = *((real *)((char *)sv + pitch * 22) + threadID_); + iSp = *((real *)((char *)sv + pitch * 23) + threadID_); + d = *((real *)((char *)sv + pitch * 24) + threadID_); + ff = *((real *)((char *)sv + pitch * 25) + threadID_); + fs = *((real *)((char *)sv + pitch * 26) + threadID_); + fcaf = *((real *)((char *)sv + pitch * 27) + threadID_); + fcas = *((real *)((char *)sv + pitch * 28) + threadID_); + jca = *((real *)((char *)sv + pitch * 29) + threadID_); + ffp = *((real *)((char *)sv + pitch * 30) + threadID_); + fcafp = *((real *)((char *)sv + pitch * 31) + threadID_); + nca = *((real *)((char *)sv + pitch * 32) + threadID_); + nca_i = *((real *)((char *)sv + pitch * 33) + threadID_); + ikr_c0 = *((real *)((char *)sv + pitch * 34) + threadID_); + ikr_c1 = *((real *)((char *)sv + pitch * 35) + threadID_); + ikr_c2 = *((real *)((char *)sv + pitch * 36) + threadID_); + ikr_i = *((real *)((char *)sv + pitch * 37) + threadID_); + ikr_o = *((real *)((char *)sv + pitch * 38) + threadID_); + xs1 = *((real *)((char *)sv + pitch * 39) + threadID_); + xs2 = *((real *)((char *)sv + pitch * 40) + threadID_); + Jrel_np = *((real *)((char *)sv + pitch * 41) + threadID_); + Jrel_p = *((real *)((char *)sv + pitch * 42) + threadID_); + } + + #include "ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments_RL.common.c" +} diff --git a/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.h b/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.h new file mode 100644 index 00000000..33b0a174 --- /dev/null +++ b/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.h @@ -0,0 +1,94 @@ +#ifndef MONOALG3D_MODEL_TORORD_FKATP_MIXED_ENDO_MID_EPI_H +#define MONOALG3D_MODEL_TORORD_FKATP_MIXED_ENDO_MID_EPI_H + +// TOMEK, Jakub et al. Development, calibration, and validation of a novel human ventricular myocyte model in health, disease, and drug block. +// Elife, v. 8, p. e48890, 2019. + +#include "../model_common.h" +#include "../../extra_data_library/helper_functions.h" + +#define NEQ 43 +#define INITIAL_V (-8.890585e+01) + +#define ENDO 0.0 +#define MID 1.0 +#define EPI 2.0 + +// CPU macros +#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_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_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_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, real *sf_Iks); + +__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 *sf_Iks, 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, bool has_extra_params); + +inline __device__ void RHS_gpu(real *sv, real *rDY_, real stim_current, real transmurality, real sf_Iks, 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 sf_Iks, 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 sf_Iks, 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 sf_Iks, 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 sf_Iks, real const *extra_params); +void solve_forward_euler_cpu_adpt(real *sv, real stim_curr, real transmurality, real sf_Iks, 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 sf_Iks, real const *extra_params); + +#endif //MONOALG3D_MODEL_TORORD_FKATP_MIXED_ENDO_MID_EPI_H + diff --git a/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments_RL.common.c b/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments_RL.common.c new file mode 100644 index 00000000..1efd0ea0 --- /dev/null +++ b/src/models_library/ToROrd/ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments_RL.common.c @@ -0,0 +1,734 @@ +// 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; + +// 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 cli = 24; // Intracellular Cl [mM] +real clo = 150; // Extracellular Cl [mM] + +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)); + +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))); + +// convenient shorthand calculations +real vffrt=v*F*F/(R*T); +real vfrt=v*F/(R*T); +real frt = F/(R*T); + +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 = 60; // 75 original. Adjustments to enable GKs to play a bigger role in APD modulation. T wave personalisation. +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 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 +// 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; + +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 +real c0 = ikr_c0; +real c1 = ikr_c1; +real c2 = ikr_c2; +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 * 0.6; // 1st element compensates for change to ko (sqrt(5/5.4)* 0.0362) Adjustment for T wave personalisation. +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 * 5.0 * IKs_Multiplier; +real GKs= 0.0011 * 5.0 * sf_Iks; // Adjustment for T wave personalisation +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); + +// 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 * INaK_Multiplier; +if (celltype==EPI) + Pnak=Pnak*0.9; +else if (celltype==2) + Pnak=Pnak*0.7; + +real INaK=Pnak*(zna*JnakNa+zk*JnakK); + +// Minor/background currents +// calculate IKb +real xkb=1.0/(1.0+exp(-(v-10.8968)/(23.9871))); +real GKb=0.0189*IKb_Multiplier; +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 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-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 +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+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; + +// 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)); +real dcai=Bcai*(-(ICaL_i + IpCa+ICab-2.0*INaCa_i)*Acap/(2.0*F*vmyo)-Jup*vnsr/vmyo+Jdiff*vss/vmyo); + +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); + +// Compute 'a' coefficients for the Hodkin-Huxley variables +a_[10] = -1.0 / tm; +a_[11] = -1.0 / tauh; +a_[12] = -1.0 / tauj; +a_[13] = -1.0 / tauh; +a_[14] = -1.0 / taujp; +a_[15] = -1.0 / tmL; +a_[16] = -1.0 / thL; +a_[17] = -1.0 / thLp; +a_[18] = -1.0 / ta; +a_[19] = -1.0 / tiF; +a_[20] = -1.0 / tiS; +a_[21] = -1.0 / ta; +a_[22] = -1.0 / tiFp; +a_[23] = -1.0 / tiSp; +a_[24] = -1.0 / td; +a_[25] = -1.0 / tff; +a_[26] = -1.0 / tfs; +a_[27] = -1.0 / tfcaf; +a_[28] = -1.0 / tfcas; +a_[29] = -1.0 / tjca; +a_[30] = -1.0 / tffp; +a_[31] = -1.0 / tfcafp; +a_[39] = -1.0 / txs1; +a_[40] = -1.0 / txs2; +a_[41] = -1.0 / tau_rel; +a_[42] = -1.0 / tau_relp; + +// Compute 'b' coefficients for the Hodkin-Huxley variables +b_[10] = mss / tm; +b_[11] = hss / tauh; +b_[12] = jss / tauj; +b_[13] = hssp / tauh; +b_[14] = jss / taujp; +b_[15] = mLss / tmL; +b_[16] = hLss / thL; +b_[17] = hLssp / thLp; +b_[18] = ass / ta; +b_[19] = iss / tiF; +b_[20] = iss / tiS; +b_[21] = assp / ta; +b_[22] = iss / tiFp; +b_[23] = iss / tiSp; +b_[24] = dss / td; +b_[25] = fss / tff; +b_[26] = fss / tfs; +b_[27] = fcass / tfcaf; +b_[28] = fcass / tfcas; +b_[29] = jcass / tjca; +b_[30] = fss / tffp; +b_[31] = fcass / tfcafp; +b_[39] = xs1ss / txs1; +b_[40] = xs2ss / txs2; +b_[41] = Jrel_inf / tau_rel; +b_[42] = Jrel_infp / tau_relp; + +// Right-hand side +rDY_[0] = dv; +rDY_[1] = dCaMKt; +rDY_[2] = dcass; +rDY_[3] = dnai; +rDY_[4] = dnass; +rDY_[5] = dki; +rDY_[6] = dkss; +rDY_[7] = dcansr; +rDY_[8] = dcajsr; +rDY_[9] = dcai; +rDY_[10] = dm; +rDY_[11] = dh; +rDY_[12] = dj; +rDY_[13] = dhp; +rDY_[14] = djp; +rDY_[15] = dmL; +rDY_[16] = dhL; +rDY_[17] = dhLp; +rDY_[18] = da; +rDY_[19] = diF; +rDY_[20] = diS; +rDY_[21] = dap; +rDY_[22] = diFp; +rDY_[23] = diSp; +rDY_[24] = dd; +rDY_[25] = dff; +rDY_[26] = dfs; +rDY_[27] = dfcaf; +rDY_[28] = dfcas; +rDY_[29] = djca; +rDY_[30] = dffp; +rDY_[31] = dfcafp; +rDY_[32] = dnca; +rDY_[33] = dnca_i; +rDY_[34] = dc0; +rDY_[35] = dc1; +rDY_[36] = dc2; +rDY_[37] = di; +rDY_[38] = delta_o; +rDY_[39] = dxs1; +rDY_[40] = dxs2; +rDY_[41] = dJrelnp; +rDY_[42] = dJrelp; diff --git a/src/models_library/ToROrd/build.sh b/src/models_library/ToROrd/build.sh index 0da7b92e..6b66e048 100644 --- a/src/models_library/ToROrd/build.sh +++ b/src/models_library/ToROrd/build.sh @@ -25,3 +25,11 @@ 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" + +############## ToRORd fkatp GKsGKrtjca adjusted Mixed ENDO_MID_EPI ############################## +MODEL_FILE_CPU="ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.c" +MODEL_FILE_GPU="ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.cu" +COMMON_HEADERS="ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.h" + +COMPILE_MODEL_LIB "ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments" "$MODEL_FILE_CPU" "$MODEL_FILE_GPU" "$COMMON_HEADERS" +