diff --git a/src/3dparty/ini_parser/ini.h b/src/3dparty/ini_parser/ini.h
index 5ce9ba06..27059a83 100644
--- a/src/3dparty/ini_parser/ini.h
+++ b/src/3dparty/ini_parser/ini.h
@@ -100,7 +100,7 @@ int ini_parse_string(const char* string, ini_handler handler, void* user);
 
 /* Maximum line length for any line in INI file. */
 #ifndef INI_MAX_LINE
-#define INI_MAX_LINE 200
+#define INI_MAX_LINE 2048
 #endif
 
 #ifdef __cplusplus
diff --git a/src/config/ecg_config.h b/src/config/ecg_config.h
index ba936312..bbe12178 100644
--- a/src/config/ecg_config.h
+++ b/src/config/ecg_config.h
@@ -11,7 +11,7 @@
 #include "../monodomain/monodomain_solver.h"
 #include "config_common.h"
 
-#define CALC_ECG(name) void name(struct time_info *time_info, struct config *config, struct grid *the_grid)
+#define CALC_ECG(name) void name(struct time_info *time_info, struct config *config, struct grid *the_grid, char *output_dir)
 typedef CALC_ECG(calc_ecg_fn);
 
 #define INIT_CALC_ECG(name) void name(struct config *config, struct ode_solver *the_ode_solver, struct grid *the_grid, char *output_dir)
diff --git a/src/domains_library/mesh_info_data.h b/src/domains_library/mesh_info_data.h
index 06c74704..7af3d886 100644
--- a/src/domains_library/mesh_info_data.h
+++ b/src/domains_library/mesh_info_data.h
@@ -12,17 +12,12 @@ struct default_fibrotic_mesh_info {
 
 // TODO: Move this struct and its macros to "custom_mesh_info_data.h"!
 struct dti_mesh_info {
-	enum dti_transmurality_labels {
-		DTI_FAST_ENDO,
-		DTI_ENDO,
-		DTI_MID,
-		DTI_EPI
-	} dti_transmurality_labels;
-	float transmurality;
-    float apicobasal;
-	float base_apex_heterogeneity;
-    float sf_IKs;
-    int fast_endo;
+		float dti_transmurality_labels;
+		float transmurality;
+    	float apicobasal;
+		float base_apex_heterogeneity;
+    	float sf_IKs;
+    	int fast_endo;
 };
 
 #define FIBROTIC_INFO(grid_cell) (struct default_fibrotic_mesh_info *)(grid_cell)->mesh_extra_info
diff --git a/src/ecg_library/ecg.c b/src/ecg_library/ecg.c
index 095ad2ea..5a5fd221 100644
--- a/src/ecg_library/ecg.c
+++ b/src/ecg_library/ecg.c
@@ -482,12 +482,12 @@ CALC_ECG(pseudo_bidomain) {
     GET_PARAMETER_BOOLEAN_VALUE_OR_USE_DEFAULT(gpu, config, "use_gpu");
     if(gpu) {
 #ifdef COMPILE_CUDA
-        pseudo_bidomain_gpu(time_info, config, the_grid);
+        pseudo_bidomain_gpu(time_info, config, the_grid, output_dir);
 #else
-        pseudo_bidomain_cpu(time_info, config, the_grid);
+        pseudo_bidomain_cpu(time_info, config, the_grid, output_dir);
 #endif
     } else {
-        pseudo_bidomain_cpu(time_info, config, the_grid);
+        pseudo_bidomain_cpu(time_info, config, the_grid, output_dir);
     }
 }
 
@@ -504,3 +504,316 @@ END_CALC_ECG(end_pseudo_bidomain) {
         end_pseudo_bidomain_cpu(config);
     }
 }
+
+INIT_CALC_ECG(init_pseudo_bidomain_with_diffusive_current_cpu) {
+    config->persistent_data = CALLOC_ONE_TYPE(struct pseudo_bidomain_persistent_data);
+
+    uint32_t nlen_output_dir = strlen(output_dir);
+    char *filename = (char*)malloc(sizeof(char)*nlen_output_dir+10);
+    sprintf(filename, "%s/ecg.txt", output_dir);
+    
+    char *dir = get_dir_from_path(filename);
+    create_dir(dir);
+    free(dir);
+
+    PSEUDO_BIDOMAIN_DATA->output_file = fopen(filename, "w");
+
+    if(PSEUDO_BIDOMAIN_DATA->output_file == NULL) {
+        log_error_and_exit("init_pseudo_bidomain - Unable to open file %s!\n", filename);
+    }
+
+    real_cpu sigma_b = 1.0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_b, config, "sigma_b");
+    uint32_t diff_curr_rate = 100;
+    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(uint32_t, diff_curr_rate, config, "diff_curr_rate");
+    PSEUDO_BIDOMAIN_DATA->diff_curr_rate = diff_curr_rate;
+    real_cpu diff_curr_max_time = 100.0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, diff_curr_max_time, config, "diff_curr_max_time");
+    PSEUDO_BIDOMAIN_DATA->diff_curr_max_time = diff_curr_max_time;
+
+    if(sigma_b == 0.0) {
+        log_error_and_exit("init_pseudo_bidomain - sigma_b can't be 0!\n");
+    }
+
+    PSEUDO_BIDOMAIN_DATA->scale_factor = 1.0 / (4.0 * M_PI * sigma_b);
+
+    get_leads(config);
+    PSEUDO_BIDOMAIN_DATA->n_leads = arrlen(PSEUDO_BIDOMAIN_DATA->leads);
+
+    uint32_t n_active = the_grid->num_active_cells;
+    struct cell_node **ac = the_grid->active_cells;
+
+    PSEUDO_BIDOMAIN_DATA->distances = MALLOC_ARRAY_OF_TYPE(real, PSEUDO_BIDOMAIN_DATA->n_leads * n_active);
+
+    PSEUDO_BIDOMAIN_DATA->beta_im = MALLOC_ARRAY_OF_TYPE(real, n_active);
+
+    // calc the distances from each volume to each electrode (r)
+    for(uint32_t i = 0; i < PSEUDO_BIDOMAIN_DATA->n_leads; i++) {
+
+        struct point_3d lead = PSEUDO_BIDOMAIN_DATA->leads[i];
+
+        OMP(parallel for)
+        for(int j = 0; j < n_active; j++) {
+            uint32_t index = i * n_active + j;
+            struct point_3d center = ac[j]->center;
+            PSEUDO_BIDOMAIN_DATA->distances[index] = EUCLIDIAN_DISTANCE(lead, center);
+        }
+    }
+
+    assembly_divergent(the_grid);
+
+    free(filename);
+}
+
+CALC_ECG(pseudo_bidomain_with_diffusive_current_cpu) {
+    // use the equation described in https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3378475/#FD7
+    uint32_t n_active = the_grid->num_active_cells;
+    struct cell_node **ac = the_grid->active_cells;
+
+    OMP(parallel for)
+    for(uint32_t i = 0; i < n_active; i++) {
+        struct element *cell_elements = ac[i]->elements;
+        size_t max_el = arrlen(cell_elements);
+
+        struct point_3d d = ac[i]->discretization;
+        real_cpu volume = d.x * d.y * d.z;
+
+        PSEUDO_BIDOMAIN_DATA->beta_im[i] = 0;
+
+        for(size_t el = 0; el < max_el; el++) {
+            PSEUDO_BIDOMAIN_DATA->beta_im[i] += cell_elements[el].value_ecg * cell_elements[el].cell->v / volume;
+        }
+    }
+
+    // DIFFUSIVE CURRENT PLOT
+    if (time_info->iteration % PSEUDO_BIDOMAIN_DATA->diff_curr_rate == 0 && time_info->current_t <= PSEUDO_BIDOMAIN_DATA->diff_curr_max_time) {
+        char *filename_diffusive_current[200];
+        sprintf(filename_diffusive_current, "%s/diffusive_current_t_%g_ms.txt", output_dir, time_info->current_t);
+        FILE *output_diffusive_current = fopen(filename_diffusive_current, "w");
+        for(uint32_t i = 0; i < n_active; i++) {
+            fprintf(output_diffusive_current, "%g\n", PSEUDO_BIDOMAIN_DATA->beta_im[i]);
+        }
+        fclose(output_diffusive_current);
+    }
+
+    fprintf(PSEUDO_BIDOMAIN_DATA->output_file, "%lf ", time_info->current_t);
+
+    for(uint32_t i = 0; i < PSEUDO_BIDOMAIN_DATA->n_leads; i++) {
+        real_cpu local_sum = 0.0;
+
+        OMP(parallel for reduction(+:local_sum))
+        for(uint32_t j = 0; j < n_active; j++) {
+            struct point_3d d = ac[j]->discretization;
+            real_cpu volume = d.x * d.y * d.z;
+
+            uint32_t index = i * n_active + j;
+            local_sum += ((PSEUDO_BIDOMAIN_DATA->beta_im[j] / PSEUDO_BIDOMAIN_DATA->distances[index])) * volume;
+        }
+
+        fprintf(PSEUDO_BIDOMAIN_DATA->output_file, "%lf ", -PSEUDO_BIDOMAIN_DATA->scale_factor * local_sum);
+    }
+
+    fprintf(PSEUDO_BIDOMAIN_DATA->output_file, "\n");
+}
+    
+END_CALC_ECG(end_pseudo_bidomain_with_diffusive_current_cpu) {
+    fclose(PSEUDO_BIDOMAIN_DATA->output_file);
+    free(PSEUDO_BIDOMAIN_DATA->beta_im);
+    arrfree(PSEUDO_BIDOMAIN_DATA->leads);
+}
+
+#ifdef COMPILE_CUDA
+
+INIT_CALC_ECG(init_pseudo_bidomain_with_diffusive_current_gpu) {
+
+    init_pseudo_bidomain_with_diffusive_current_cpu(config, NULL, the_grid, output_dir);
+
+    if(!the_ode_solver->gpu) {
+        log_warn("The current implementation of pseudo_bidomain_gpu only works when the odes are also being solved using the GPU! Falling back to CPU version\n");
+        shput(config->config_data, "use_gpu", "false");
+        return;
+    }
+
+    //This is allocated when using the CPU code, but we do not need it in the gpu version
+    free(PSEUDO_BIDOMAIN_DATA->beta_im);
+
+    check_cublas_error(cusparseCreate(&(PSEUDO_BIDOMAIN_DATA->cusparseHandle)));
+    check_cublas_error(cublasCreate(&(PSEUDO_BIDOMAIN_DATA->cublasHandle)));
+
+    int_array I = NULL, J = NULL;
+    f32_array val = NULL;
+
+    grid_to_csr_for_ecg(the_grid, &val, &I, &J, false, true);
+
+    int nz = arrlen(val);
+    uint32_t N = the_grid->num_active_cells;
+
+    // DIFFUSIVE CURRENT
+    // Allocate memory for the CPU beta_im
+    PSEUDO_BIDOMAIN_DATA->beta_im_cpu = MALLOC_ARRAY_OF_TYPE(real, N);
+
+    // We need to do that because val is originally a float value. Maybe we should make it real.
+    real *new_val = NULL;
+    arrsetlen(new_val, nz);
+    for(int i = 0; i < nz; i++) {
+        new_val[i] = (real)val[i];
+    }
+
+    PSEUDO_BIDOMAIN_DATA->nz = nz;
+
+    check_cuda_error(cudaMalloc((void **)&(PSEUDO_BIDOMAIN_DATA->d_col), nz * sizeof(int)));
+    check_cuda_error(cudaMalloc((void **)&(PSEUDO_BIDOMAIN_DATA->d_row), (N + 1) * sizeof(int)));
+    check_cuda_error(cudaMalloc((void **)&(PSEUDO_BIDOMAIN_DATA->d_val), nz * sizeof(real)));
+    check_cuda_error(cudaMalloc((void **)&(PSEUDO_BIDOMAIN_DATA->beta_im), N * sizeof(real)));
+    check_cuda_error(cudaMalloc((void **)&(PSEUDO_BIDOMAIN_DATA->d_distances), PSEUDO_BIDOMAIN_DATA->n_leads * N * sizeof(real)));
+    check_cuda_error(cudaMalloc((void **)&(PSEUDO_BIDOMAIN_DATA->d_volumes), N * sizeof(real)));
+    check_cuda_error(cudaMalloc((void **)&(PSEUDO_BIDOMAIN_DATA->tmp_data), N * sizeof(real)));
+
+#if CUBLAS_VER_MAJOR >= 11
+    check_cuda_error(cusparseCreateCsr(&(PSEUDO_BIDOMAIN_DATA->matA), N, N, nz, PSEUDO_BIDOMAIN_DATA->d_row, PSEUDO_BIDOMAIN_DATA->d_col,
+                                       PSEUDO_BIDOMAIN_DATA->d_val, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUBLAS_SIZE));
+    check_cuda_error(cusparseCreateDnVec(&(PSEUDO_BIDOMAIN_DATA->vec_beta_im), N, PSEUDO_BIDOMAIN_DATA->beta_im, CUBLAS_SIZE));
+    check_cuda_error(cusparseCreateDnVec(&(PSEUDO_BIDOMAIN_DATA->vec_vm), N, the_ode_solver->sv, CUBLAS_SIZE));
+#else
+    check_cuda_error((cudaError_t)cusparseCreateMatDescr(&(PSEUDO_BIDOMAIN_DATA->descr)));
+    cusparseSetMatType(PSEUDO_BIDOMAIN_DATA->descr, CUSPARSE_MATRIX_TYPE_GENERAL);
+    cusparseSetMatIndexBase(PSEUDO_BIDOMAIN_DATA->descr, CUSPARSE_INDEX_BASE_ZERO);
+    PSEUDO_BIDOMAIN_DATA->local_sv = the_ode_solver->sv;
+#endif
+
+    check_cuda_error(cudaMemcpy(PSEUDO_BIDOMAIN_DATA->d_col, J, nz * sizeof(int), cudaMemcpyHostToDevice));        // JA
+    check_cuda_error(cudaMemcpy(PSEUDO_BIDOMAIN_DATA->d_row, I, (N + 1) * sizeof(int), cudaMemcpyHostToDevice));   // IA
+    check_cuda_error(cudaMemcpy(PSEUDO_BIDOMAIN_DATA->d_val, new_val, nz * sizeof(real), cudaMemcpyHostToDevice)); // A
+    check_cuda_error(cudaMemcpy(PSEUDO_BIDOMAIN_DATA->d_distances, PSEUDO_BIDOMAIN_DATA->distances, PSEUDO_BIDOMAIN_DATA->n_leads * N * sizeof(real),
+                                cudaMemcpyHostToDevice));
+
+    PSEUDO_BIDOMAIN_DATA->volumes = MALLOC_ARRAY_OF_TYPE(real, N);
+
+    struct cell_node **ac = the_grid->active_cells;
+    OMP(parallel for)
+    for(int i = 0; i < N; i += 1) {
+        struct point_3d d = ac[i]->discretization;
+        PSEUDO_BIDOMAIN_DATA->volumes[i] = d.x * d.y * d.z;
+    }
+
+    check_cuda_error(cudaMemcpy(PSEUDO_BIDOMAIN_DATA->d_volumes, PSEUDO_BIDOMAIN_DATA->volumes, N * sizeof(real), cudaMemcpyHostToDevice));
+
+#if CUSPARSE_VER_MAJOR >= 11
+    real alpha = 1.0;
+    real beta = 0.0;
+
+    check_cuda_error(cusparseSpMV_bufferSize(PSEUDO_BIDOMAIN_DATA->cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, PSEUDO_BIDOMAIN_DATA->matA,
+                                             PSEUDO_BIDOMAIN_DATA->vec_vm, &beta, PSEUDO_BIDOMAIN_DATA->vec_beta_im, CUBLAS_SIZE, CUSPARSE_ALG,
+                                             &(PSEUDO_BIDOMAIN_DATA->bufferSize)));
+
+    check_cuda_error(cudaMalloc(&(PSEUDO_BIDOMAIN_DATA->buffer), PSEUDO_BIDOMAIN_DATA->bufferSize));
+#endif
+
+    arrfree(I);
+    arrfree(J);
+    arrfree(val);
+    arrfree(new_val);
+
+    free(PSEUDO_BIDOMAIN_DATA->volumes);
+    free(PSEUDO_BIDOMAIN_DATA->distances);
+}
+
+CALC_ECG(pseudo_bidomain_with_diffusive_current_gpu) {
+
+    uint32_t n_active = the_grid->num_active_cells;
+
+    // VM is correct
+    real alpha = 1.0;
+    real beta = 0.0;
+#if CUSPARSE_VER_MAJOR >= 11
+    check_cublas_error(cusparseSpMV(PSEUDO_BIDOMAIN_DATA->cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, PSEUDO_BIDOMAIN_DATA->matA,
+                                    PSEUDO_BIDOMAIN_DATA->vec_vm, &beta, PSEUDO_BIDOMAIN_DATA->vec_beta_im, CUBLAS_SIZE, CUSPARSE_ALG,
+                                    PSEUDO_BIDOMAIN_DATA->buffer));
+#else
+
+#ifdef CELL_MODEL_REAL_DOUBLE
+    cusparseDcsrmv(PSEUDO_BIDOMAIN_DATA->cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, n_active, n_active, PSEUDO_BIDOMAIN_DATA->nz, &alpha,
+                   PSEUDO_BIDOMAIN_DATA->descr, PSEUDO_BIDOMAIN_DATA->d_val, PSEUDO_BIDOMAIN_DATA->d_row, PSEUDO_BIDOMAIN_DATA->d_col,
+                   PSEUDO_BIDOMAIN_DATA->local_sv, &beta, PSEUDO_BIDOMAIN_DATA->beta_im);
+#else
+    cusparseScsrmv(PSEUDO_BIDOMAIN_DATA->cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, n_active, n_active, PSEUDO_BIDOMAIN_DATA->nz, &alpha,
+                   PSEUDO_BIDOMAIN_DATA->descr, PSEUDO_BIDOMAIN_DATA->d_val, PSEUDO_BIDOMAIN_DATA->d_row, PSEUDO_BIDOMAIN_DATA->d_col,
+                   PSEUDO_BIDOMAIN_DATA->local_sv, &beta, PSEUDO_BIDOMAIN_DATA->beta_im);
+#endif
+
+#endif
+
+    fprintf(PSEUDO_BIDOMAIN_DATA->output_file, "%lf ", time_info->current_t);
+
+    gpu_vec_div_vec(PSEUDO_BIDOMAIN_DATA->beta_im, PSEUDO_BIDOMAIN_DATA->d_volumes, PSEUDO_BIDOMAIN_DATA->beta_im, n_active);
+
+    // DIFFUSIVE CURRENT
+    if (time_info->iteration % PSEUDO_BIDOMAIN_DATA->diff_curr_rate == 0 && time_info->current_t <= PSEUDO_BIDOMAIN_DATA->diff_curr_max_time) {
+        char *filename_diffusive_current[200];
+        sprintf(filename_diffusive_current, "%s/diffusive_current_t_%g_ms.txt", output_dir, time_info->current_t);
+        FILE *output_diffusive_current = fopen(filename_diffusive_current, "w");
+        check_cuda_error(cudaMemcpy(PSEUDO_BIDOMAIN_DATA->beta_im_cpu, PSEUDO_BIDOMAIN_DATA->beta_im, n_active * sizeof(real), cudaMemcpyDeviceToHost));
+        for(uint32_t i = 0; i < n_active; i++) {
+            fprintf(output_diffusive_current, "%g\n", PSEUDO_BIDOMAIN_DATA->beta_im_cpu[i]);
+        }
+        fclose(output_diffusive_current);
+    }
+
+    for(int i = 0; i < PSEUDO_BIDOMAIN_DATA->n_leads; i++) {
+
+        // beta_im / distance
+        gpu_vec_div_vec(PSEUDO_BIDOMAIN_DATA->beta_im, PSEUDO_BIDOMAIN_DATA->d_distances + n_active * i, PSEUDO_BIDOMAIN_DATA->tmp_data, n_active);
+        real local_sum;
+
+#ifdef CELL_MODEL_REAL_DOUBLE
+        check_cublas_error(
+            cublasDdot(PSEUDO_BIDOMAIN_DATA->cublasHandle, n_active, PSEUDO_BIDOMAIN_DATA->tmp_data, 1, PSEUDO_BIDOMAIN_DATA->d_volumes, 1, &local_sum));
+#else
+        check_cublas_error(
+            cublasSdot(PSEUDO_BIDOMAIN_DATA->cublasHandle, n_active, PSEUDO_BIDOMAIN_DATA->tmp_data, 1, PSEUDO_BIDOMAIN_DATA->d_volumes, 1, &local_sum));
+#endif
+
+        fprintf(PSEUDO_BIDOMAIN_DATA->output_file, "%lf ", -PSEUDO_BIDOMAIN_DATA->scale_factor * local_sum);
+    }
+
+    fprintf(PSEUDO_BIDOMAIN_DATA->output_file, "\n");
+
+}
+
+END_CALC_ECG(end_pseudo_bidomain_with_diffusive_current_gpu) {
+
+    struct pseudo_bidomain_persistent_data *persistent_data = (struct pseudo_bidomain_persistent_data *)config->persistent_data;
+
+    if(!persistent_data)
+        return;
+
+    check_cuda_error((cudaError_t)cusparseDestroy(persistent_data->cusparseHandle));
+    check_cuda_error((cudaError_t)cublasDestroy(persistent_data->cublasHandle));
+
+#if CUBLAS_VER_MAJOR >= 11
+    if(persistent_data->matA) {
+        check_cuda_error(cusparseDestroySpMat(persistent_data->matA));
+    }
+    if(persistent_data->vec_beta_im) {
+        check_cuda_error(cusparseDestroyDnVec(persistent_data->vec_beta_im));
+    }
+    if(persistent_data->vec_vm) {
+        check_cuda_error(cusparseDestroyDnVec(persistent_data->vec_vm));
+    }
+#else
+    check_cuda_error((cudaError_t)cusparseDestroyMatDescr(persistent_data->descr));
+#endif
+    check_cuda_error(cudaFree(persistent_data->d_col));
+    check_cuda_error(cudaFree(persistent_data->d_row));
+    check_cuda_error(cudaFree(persistent_data->d_val));
+    check_cuda_error(cudaFree(persistent_data->beta_im));
+    check_cuda_error(cudaFree(persistent_data->tmp_data));
+    check_cuda_error(cudaFree(persistent_data->d_distances));
+    check_cuda_error(cudaFree(persistent_data->d_volumes));
+
+    free(persistent_data->beta_im_cpu);
+    fclose(persistent_data->output_file);
+
+    free(persistent_data);
+}
+#endif
\ No newline at end of file
diff --git a/src/ecg_library/ecg.h b/src/ecg_library/ecg.h
index 11e2ae76..3d6b5bc3 100644
--- a/src/ecg_library/ecg.h
+++ b/src/ecg_library/ecg.h
@@ -13,6 +13,8 @@ struct pseudo_bidomain_persistent_data {
     FILE *output_file;
     real_cpu scale_factor;
     uint32_t n_leads;
+    uint32_t diff_curr_rate;
+    real_cpu diff_curr_max_time;
 
 #ifdef COMPILE_CUDA
     cusparseHandle_t cusparseHandle;
@@ -23,6 +25,7 @@ struct pseudo_bidomain_persistent_data {
     real *volumes;
     real *tmp_data;
     real *d_val;
+    real *beta_im_cpu;
 
     size_t bufferSize;
     void *buffer;
diff --git a/src/extra_data_library/helper_functions.c b/src/extra_data_library/helper_functions.c
index f2e42ae8..e2140f0f 100644
--- a/src/extra_data_library/helper_functions.c
+++ b/src/extra_data_library/helper_functions.c
@@ -1031,6 +1031,8 @@ struct extra_data_for_torord_gksgkrtjca_twave * set_common_torord_gksgkrtjca_twa
     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");
+	real cycle_length = 1000;
+	GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, cycle_length, config, "cycle_length");
 
     extra_data->INa_Multiplier = INa_Multiplier;
     extra_data->INaL_Multiplier = INaL_Multiplier;
@@ -1057,140 +1059,280 @@ struct extra_data_for_torord_gksgkrtjca_twave * set_common_torord_gksgkrtjca_twa
     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;
+    // (DTI032) Basic cycle length => 810ms
+    if (cycle_length == 810) {
+        printf("[ToRORd] Setting hacked initial conditions for DTI032 CL 810!!!");
+        extra_data->initial_ss_endo[0 ]  = -88.676; // Vm
+        extra_data->initial_ss_endo[1 ]  = 12.084; // Nai
+        extra_data->initial_ss_endo[2 ]  = 12.085; // Nass
+        extra_data->initial_ss_endo[3 ]  = 141.28; // ki
+        extra_data->initial_ss_endo[4 ]  = 141.28; // kss
+        extra_data->initial_ss_endo[5 ]  = 7.456e-05; //cai
+        extra_data->initial_ss_endo[6 ]  = 6.504164e-05;  // cass
+        extra_data->initial_ss_endo[7 ]  = 1.5145; //cansr
+        extra_data->initial_ss_endo[8 ]  = 1.5; // cajsr
+        extra_data->initial_ss_endo[9] = 0.00082113; // m
+        extra_data->initial_ss_endo[10] = 0.66775; // hp
+        extra_data->initial_ss_endo[11] = 0.82677; // h
+        extra_data->initial_ss_endo[12] = 0.82599; // j
+        extra_data->initial_ss_endo[13] = 0.82257; // jp
+        extra_data->initial_ss_endo[14] = 0.00016564; // mL
+        extra_data->initial_ss_endo[15] = 0.49351; // hL
+        extra_data->initial_ss_endo[16] = 0.24641; // hLp
+        extra_data->initial_ss_endo[17] = 0.00095665; // a
+        extra_data->initial_ss_endo[18] = 0.9996; // iF
+        extra_data->initial_ss_endo[19] = 0.46261; // iS
+        extra_data->initial_ss_endo[20] = 0.00048743; // ap
+        extra_data->initial_ss_endo[21] = 0.9996; // iFp
+        extra_data->initial_ss_endo[22] = 0.51774; // iSp
+        extra_data->initial_ss_endo[23] = -5.22e-27; // d
+        extra_data->initial_ss_endo[24] = 1.000000e+00; // ff
+        extra_data->initial_ss_endo[25] = 9.0939e-01; // fs
+        extra_data->initial_ss_endo[26] = 1.000000e+00; // fcaf
+        extra_data->initial_ss_endo[27] = 9.9902e-01; // fcas
+        extra_data->initial_ss_endo[28] = 9.9995e-01; // jca
+        extra_data->initial_ss_endo[29] = 5.032e-04; // nca
+        extra_data->initial_ss_endo[30] = 8.3459e-04; // nca_i
+        extra_data->initial_ss_endo[31] = 1.000000e+00; // ffp
+        extra_data->initial_ss_endo[32] = 1.000000e+00; // fcafp
+        extra_data->initial_ss_endo[33] = 0.30978; // xs1
+        extra_data->initial_ss_endo[34] = 0.00017898; // xs2
+        extra_data->initial_ss_endo[35] = 1.1046e-22; // Jrel_np
+        extra_data->initial_ss_endo[36]  = 0.014729; // CaMKt
+        extra_data->initial_ss_endo[37] = 0.99455; // ikr_c0
+        extra_data->initial_ss_endo[38] = 0.00085423; // ikr_c1
+        extra_data->initial_ss_endo[39] = 0.00091282; // ikr_c2
+        extra_data->initial_ss_endo[40] = 0.0035557; // ikr_o
+        extra_data->initial_ss_endo[41] = 0.00012677; // ikr_i
+        extra_data->initial_ss_endo[42] = -3.6032e-20; // Jrel_p
+    }
+    else if (cycle_length == 909) {
+        printf("[ToRORd] Setting hacked initial conditions for DTI024 CL 909!!!");
+        extra_data->initial_ss_endo[0 ]  = -88.606; // Vm
+        extra_data->initial_ss_endo[1 ]  = 11.951; // Nai
+        extra_data->initial_ss_endo[2 ]  = 11.951; // Nass
+        extra_data->initial_ss_endo[3 ]  = 141.37; // ki
+        extra_data->initial_ss_endo[4 ]  = 141.37; // kss
+        extra_data->initial_ss_endo[5 ]  = 7.3303e-05; //cai
+        extra_data->initial_ss_endo[6 ]  = 6.4258e-05;  // cass
+        extra_data->initial_ss_endo[7 ]  = 1.4869; //cansr
+        extra_data->initial_ss_endo[8 ]  = 1.4799; // cajsr
+        extra_data->initial_ss_endo[9 ] = 0.00082298; // m
+        extra_data->initial_ss_endo[10] = 0.66742; // hp
+        extra_data->initial_ss_endo[11] = 0.82657; // h
+        extra_data->initial_ss_endo[12] = 0.8262; // j
+        extra_data->initial_ss_endo[13] = 0.82484; // jp
+        extra_data->initial_ss_endo[14] = 0.00016597; // mL
+        extra_data->initial_ss_endo[15] = 0.50798; // hL
+        extra_data->initial_ss_endo[16] = 0.26085; // hLp
+        extra_data->initial_ss_endo[17] = 0.00095732; // a
+        extra_data->initial_ss_endo[18] = 0.9996; // iF
+        extra_data->initial_ss_endo[19] = 0.51827; // iS
+        extra_data->initial_ss_endo[20] = 0.00048777; // ap
+        extra_data->initial_ss_endo[21] = 0.9996; // iFp
+        extra_data->initial_ss_endo[22] = 0.57465; // iSp
+        extra_data->initial_ss_endo[23] = -3.7135e-25; // d
+        extra_data->initial_ss_endo[24] = 1.000000e+00; // ff
+        extra_data->initial_ss_endo[25] = 9.2181e-01; // fs
+        extra_data->initial_ss_endo[26] = 1.000000e+00; // fcaf
+        extra_data->initial_ss_endo[27] = 9.996e-01; // fcas
+        extra_data->initial_ss_endo[28] = 0.99999; // jca
+        extra_data->initial_ss_endo[29] = 4.6963e-04; // nca
+        extra_data->initial_ss_endo[30] = 7.8127e-04; // nca_i
+        extra_data->initial_ss_endo[31] = 1.000000e+00; // ffp
+        extra_data->initial_ss_endo[32] = 1.000000e+00; // fcafp
+        extra_data->initial_ss_endo[33] = 0.2845; // xs1
+        extra_data->initial_ss_endo[34] = 0.0001791; // xs2
+        extra_data->initial_ss_endo[35] = 8.7097e-22; // Jrel_np
+        extra_data->initial_ss_endo[36]  = 0.011988; // CaMKt  
+        extra_data->initial_ss_endo[37] = 0.99695; // ikr_c0
+        extra_data->initial_ss_endo[38] = 0.00085638; // ikr_c1
+        extra_data->initial_ss_endo[39] = 0.00077257; // ikr_c2
+        extra_data->initial_ss_endo[40] = 0.00013754; // ikr_o
+        extra_data->initial_ss_endo[41] = 4.9031e-05; // ikr_i
+        extra_data->initial_ss_endo[42] = -1.724e-20; // Jrel_p
+    } else if (cycle_length == 1250) {
+	    printf("[ToRORd] Setting hacked initial conditions for DTI004 CL 1250!!!");
+        extra_data->initial_ss_endo[0 ]  = -88.606; // Vm
+        extra_data->initial_ss_endo[1 ]  = 11.523; // Nai
+        extra_data->initial_ss_endo[2 ]  = 11.523; // Nass
+        extra_data->initial_ss_endo[3 ]  = 141.63; // ki
+        extra_data->initial_ss_endo[4 ]  = 141.63; // kss
+        extra_data->initial_ss_endo[5 ]  = 7.1347e-05; //cai
+        extra_data->initial_ss_endo[6 ]  = 6.2589e-05;  // cass
+        extra_data->initial_ss_endo[7 ]  = 1.4164; //cansr
+        extra_data->initial_ss_endo[8 ]  = 1.4177; // cajsr
+        extra_data->initial_ss_endo[9 ] = 0.00083359; // m
+        extra_data->initial_ss_endo[10] = 0.66549; // hp
+        extra_data->initial_ss_endo[11] = 0.82538; // h
+        extra_data->initial_ss_endo[12] = 0.82534; // j
+        extra_data->initial_ss_endo[13] = 0.82526; // jp
+        extra_data->initial_ss_endo[14] = 0.00016786; // mL
+        extra_data->initial_ss_endo[15] = 0.52741; // hL
+        extra_data->initial_ss_endo[16] = 0.2928; // hLp
+        extra_data->initial_ss_endo[17] = 0.00096116; // a
+        extra_data->initial_ss_endo[18] = 0.9996; // iF
+        extra_data->initial_ss_endo[19] = 0.67268; // iS
+        extra_data->initial_ss_endo[20] = 0.00048973; // ap
+        extra_data->initial_ss_endo[21] = 0.9996; // iFp
+        extra_data->initial_ss_endo[22] = 0.72701; // iSp
+        extra_data->initial_ss_endo[23] = -7.7171e-27; // d
+        extra_data->initial_ss_endo[24] = 1.000000e+00; // ff
+        extra_data->initial_ss_endo[25] = 9.5029e-01; // fs
+        extra_data->initial_ss_endo[26] = 1.000000e+00; // fcaf
+        extra_data->initial_ss_endo[27] = 9.998e-01; // fcas
+        extra_data->initial_ss_endo[28] = 1.0; // jca
+        extra_data->initial_ss_endo[29] = 4.238e-04; // nca
+        extra_data->initial_ss_endo[30] = 7.0341e-04; // nca_i
+        extra_data->initial_ss_endo[31] = 1.000000e+00; // ffp
+        extra_data->initial_ss_endo[32] = 1.000000e+00; // fcafp
+        extra_data->initial_ss_endo[33] = 0.21689; // xs1
+        extra_data->initial_ss_endo[34] = 0.00018021; // xs2
+        extra_data->initial_ss_endo[35] = -5.1079e-25; // Jrel_np
+        extra_data->initial_ss_endo[36]  = 0.0068498; // CaMKt  
+        extra_data->initial_ss_endo[37] = 0.99827; // ikr_c0
+        extra_data->initial_ss_endo[38] = 0.00086086; // ikr_c1
+        extra_data->initial_ss_endo[39] = 0.00069758; // ikr_c2
+        extra_data->initial_ss_endo[40] = 0.00016919; // ikr_o
+        extra_data->initial_ss_endo[41] = 6.0064e-06; // ikr_i
+        extra_data->initial_ss_endo[42] = -8.5262e-21; // Jrel_p
+    }else {
+        printf("[ToRORd] Setting default initial conditions!!!");
+        // Set the default initial conditions from Matlab (celltype = ENDO)
+        extra_data->initial_ss_endo[0]  = -88.6369922306458;
+        extra_data->initial_ss_endo[1]  = 11.8973412949238;
+        extra_data->initial_ss_endo[2]  = 11.897661047085;
+        extra_data->initial_ss_endo[3]  = 141.234464714982;
+        extra_data->initial_ss_endo[4]  = 141.234423402713;
+        extra_data->initial_ss_endo[5]  = 7.26747296460659e-05;
+        extra_data->initial_ss_endo[6]  = 6.33786975780735e-05;
+        extra_data->initial_ss_endo[7]  = 1.5326530637197;
+        extra_data->initial_ss_endo[8]  = 1.53394579180493;
+        extra_data->initial_ss_endo[9]  = 0.000828007761976018;
+        extra_data->initial_ss_endo[10] = 0.666527193684116;
+        extra_data->initial_ss_endo[11] = 0.826020806005678;
+        extra_data->initial_ss_endo[12] = 0.826055985895856;
+        extra_data->initial_ss_endo[13] = 0.825850881115628;
+        extra_data->initial_ss_endo[14] = 0.000166868626513013;
+        extra_data->initial_ss_endo[15] = 0.522830604669169;
+        extra_data->initial_ss_endo[16] = 0.285969584294187;
+        extra_data->initial_ss_endo[17] = 0.000959137028030184;
+        extra_data->initial_ss_endo[18] = 0.999601150012565;
+        extra_data->initial_ss_endo[19] = 0.5934016398361;
+        extra_data->initial_ss_endo[20] = 0.000488696137242056;
+        extra_data->initial_ss_endo[21] = 0.999601147267179;
+        extra_data->initial_ss_endo[22] = 0.654668660159696;
+        extra_data->initial_ss_endo[23] = 9.50007519781516e-32;
+        extra_data->initial_ss_endo[24] = 0.999999992317577;
+        extra_data->initial_ss_endo[25] = 0.939258048397962;
+        extra_data->initial_ss_endo[26] = 0.999999992317557;
+        extra_data->initial_ss_endo[27] = 0.999898379647465;
+        extra_data->initial_ss_endo[28] = 0.99997825156004;
+        extra_data->initial_ss_endo[29] = 0.000444816183420527;
+        extra_data->initial_ss_endo[30] = 0.000755072490632667;
+        extra_data->initial_ss_endo[31] = 0.999999992318446;
+        extra_data->initial_ss_endo[32] = 0.999999992318445;
+        extra_data->initial_ss_endo[33] = 0.24240468344952;
+        extra_data->initial_ss_endo[34] = 0.000179537726989804;
+        extra_data->initial_ss_endo[35] = -6.88308558109975e-25;
+        extra_data->initial_ss_endo[36] = 0.0111749845355653;
+        extra_data->initial_ss_endo[37] = 0.998036620213316;
+        extra_data->initial_ss_endo[38] = 0.000858801779013532;
+        extra_data->initial_ss_endo[39] = 0.000709744678350176;
+        extra_data->initial_ss_endo[40] = 0.000381261722195702;
+        extra_data->initial_ss_endo[41] = 1.35711566929992e-05;
+        extra_data->initial_ss_endo[42] = 2.30252452954649e-23;
+
+        // Set the default initial conditions from Matlab (celltype = EPI)
+        extra_data->initial_ss_epi[0] = -89.0462806262884;
+        extra_data->initial_ss_epi[1] = 12.7218980311997;
+        extra_data->initial_ss_epi[2] = 12.7222039977392;
+        extra_data->initial_ss_epi[3] = 142.248960281735;
+        extra_data->initial_ss_epi[4] = 142.248911688304;
+        extra_data->initial_ss_epi[5] = 6.54105789316085e-05;
+        extra_data->initial_ss_epi[6] = 5.68443136844764e-05;
+        extra_data->initial_ss_epi[7] = 1.80911728399381;
+        extra_data->initial_ss_epi[8] = 1.80970235621251;
+        extra_data->initial_ss_epi[9] = 0.000758182108180449;
+        extra_data->initial_ss_epi[10] = 0.679839847935577;
+        extra_data->initial_ss_epi[11] = 0.834150231581688;
+        extra_data->initial_ss_epi[12] = 0.834188252920967;
+        extra_data->initial_ss_epi[13] = 0.834081731522592;
+        extra_data->initial_ss_epi[14] = 0.000154387698861246;
+        extra_data->initial_ss_epi[15] = 0.538295069820379;
+        extra_data->initial_ss_epi[16] = 0.302769394159465;
+        extra_data->initial_ss_epi[17] = 0.000933035060391086;
+        extra_data->initial_ss_epi[18] = 0.999628705730844;
+        extra_data->initial_ss_epi[19] = 0.999626204093615;
+        extra_data->initial_ss_epi[20] = 0.00047539066209218;
+        extra_data->initial_ss_epi[21] = 0.999628705544664;
+        extra_data->initial_ss_epi[22] = 0.999628513430851;
+        extra_data->initial_ss_epi[23] = 1.74213411952898e-37;
+        extra_data->initial_ss_epi[24] = 0.999999993122906;
+        extra_data->initial_ss_epi[25] = 0.947952168523141;
+        extra_data->initial_ss_epi[26] = 0.999999993122889;
+        extra_data->initial_ss_epi[27] = 0.999932686646139;
+        extra_data->initial_ss_epi[28] = 0.999982915381882;
+        extra_data->initial_ss_epi[29] = 0.000291544679470133;
+        extra_data->initial_ss_epi[30] = 0.000502604507932921;
+        extra_data->initial_ss_epi[31] = 0.999999993124187;
+        extra_data->initial_ss_epi[32] = 0.999999993123756;
+        extra_data->initial_ss_epi[33] = 0.22881550094027;
+        extra_data->initial_ss_epi[34] = 0.000171497784228012;
+        extra_data->initial_ss_epi[35] = -1.13118992668881e-26;
+        extra_data->initial_ss_epi[36] = 0.0129505221481656;
+        extra_data->initial_ss_epi[37] = 0.998194356754674;
+        extra_data->initial_ss_epi[38] = 0.000834232097912889;
+        extra_data->initial_ss_epi[39] = 0.000683865770895308;
+        extra_data->initial_ss_epi[40] = 0.00027787850109644;
+        extra_data->initial_ss_epi[41] = 9.66775862738005e-06;
+        extra_data->initial_ss_epi[42] = 8.16930403133409e-24;
+
+        // Set the default initial conditions from Matlab (celltype = MCELL)
+        extra_data->initial_ss_mid[0] = -89.5379994049964;
+        extra_data->initial_ss_mid[1] = 14.9292004720038;
+        extra_data->initial_ss_mid[2] = 14.9296673679334;
+        extra_data->initial_ss_mid[3] = 144.84471868881;
+        extra_data->initial_ss_mid[4] = 144.844658476157;
+        extra_data->initial_ss_mid[5] = 7.50228807455408e-05;
+        extra_data->initial_ss_mid[6] = 6.10763598140135e-05;
+        extra_data->initial_ss_mid[7] = 1.79043480744558;
+        extra_data->initial_ss_mid[8] = 1.79484249993962;
+        extra_data->initial_ss_mid[9] = 0.000681936485046493;
+        extra_data->initial_ss_mid[10] = 0.695380653101535;
+        extra_data->initial_ss_mid[11] = 0.843488797335149;
+        extra_data->initial_ss_mid[12] = 0.843520761455969;
+        extra_data->initial_ss_mid[13] = 0.843226224403045;
+        extra_data->initial_ss_mid[14] = 0.000140621109700401;
+        extra_data->initial_ss_mid[15] = 0.545314876174586;
+        extra_data->initial_ss_mid[16] = 0.292496735833565;
+        extra_data->initial_ss_mid[17] = 0.000902612655601118;
+        extra_data->initial_ss_mid[18] = 0.999659345906191;
+        extra_data->initial_ss_mid[19] = 0.56311967936689;
+        extra_data->initial_ss_mid[20] = 0.000459883274920751;
+        extra_data->initial_ss_mid[21] = 0.999659343029625;
+        extra_data->initial_ss_mid[22] = 0.623696443871387;
+        extra_data->initial_ss_mid[23] = -1.31418873360667e-33;
+        extra_data->initial_ss_mid[24] = 0.999999993979673;
+        extra_data->initial_ss_mid[25] = 0.920408593154793;
+        extra_data->initial_ss_mid[26] = 0.999999993979652;
+        extra_data->initial_ss_mid[27] = 0.999761950174748;
+        extra_data->initial_ss_mid[28] = 0.999962530196306;
+        extra_data->initial_ss_mid[29] = 0.0003853594696671;
+        extra_data->initial_ss_mid[30] = 0.000853529194511867;
+        extra_data->initial_ss_mid[31] = 0.999999993978835;
+        extra_data->initial_ss_mid[32] = 0.999999993980401;
+        extra_data->initial_ss_mid[33] = 0.266415111925392;
+        extra_data->initial_ss_mid[34] = 0.000162310655612839;
+        extra_data->initial_ss_mid[35] = 1.20976169203982e-24;
+        extra_data->initial_ss_mid[36] = 0.0178243652102213;
+        extra_data->initial_ss_mid[37] = 0.997971986641796;
+        extra_data->initial_ss_mid[38] = 0.000805399061926759;
+        extra_data->initial_ss_mid[39] = 0.000678179976274546;
+        extra_data->initial_ss_mid[40] = 0.00052653630893167;
+        extra_data->initial_ss_mid[41] = 1.78956481798154e-05;
+        extra_data->initial_ss_mid[42] = 7.0591623795627e-23;
+    }
 
     extra_data->transmurality = MALLOC_ARRAY_OF_TYPE(real, num_cells);
     extra_data->sf_IKs = MALLOC_ARRAY_OF_TYPE(real, num_cells);
@@ -1293,3 +1435,4 @@ struct extra_data_for_trovato * set_common_trovato_data (struct config *config,
     return extra_data;
 }
 
+
diff --git a/src/extra_data_library/helper_functions.h b/src/extra_data_library/helper_functions.h
index d7fe01d5..4c00b31c 100644
--- a/src/extra_data_library/helper_functions.h
+++ b/src/extra_data_library/helper_functions.h
@@ -39,6 +39,8 @@ struct extra_data_for_tt3 {
     real *transmurality;
 };
 
+
+
 struct extra_data_for_torord {
     real INa_Multiplier; 
     real ICaL_Multiplier;
@@ -191,4 +193,5 @@ struct extra_data_for_torord_land_twave * set_common_torord_Land_twave_data (str
 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);
 
+
 #endif // MONOALG3D_C_EXTRA_DATA_HELPER_FUNCTIONS_H
diff --git a/src/main_clip_mesh.c b/src/main_clip_mesh.c
index 77ce2143..41d49096 100644
--- a/src/main_clip_mesh.c
+++ b/src/main_clip_mesh.c
@@ -87,22 +87,17 @@ static void expand_scar(const char *input, const char *output, int n_rows, char*
     struct config *domain_config;
 
     domain_config = alloc_and_init_config_data();
-    char *discretization = "300";
-    //char *discretization = "500";
+    char *discretization = "500";
 
-    shput_dup_value(domain_config->config_data, "start_discretization", strdup(discretization));
-    shput_dup_value(domain_config->config_data, "maximum_discretization", strdup(discretization));
+    shput_dup_value(domain_config->config_data, "original_discretization", strdup(discretization));
+    shput_dup_value(domain_config->config_data, "desired_discretization", strdup(discretization));
     shput_dup_value(domain_config->config_data, "mesh_file", strdup(input));
 
-    domain_config->main_function_name = strdup("initialize_grid_with_dti_mesh");
-    shput_dup_value(domain_config->config_data, "name", "Oxford DTI004 with Transmurality and Fiber orientation");
-
-    //shput(domain_config->config_data, "side_length_x", strdup(discretization));
-    //shput(domain_config->config_data, "side_length_y", strdup(discretization));
-    //shput(domain_config->config_data, "side_length_z", strdup(discretization));
+    domain_config->main_function_name = strdup("initialize_grid_with_dti_mesh_with_dense_and_sparse_regions_twave");
+    shput_dup_value(domain_config->config_data, "name", "Oxford DTI004 with Transmurality, Fiber orientation and Dense/Sparse regions");
+    
     shput(domain_config->config_data, "num_volumes", strdup(n_volumes));
-    //shput(domain_config->config_data, "num_extra_fields", "5");
-
+    
     init_config_functions(domain_config, "./shared_libs/libdefault_domains.so", "domain");
 
     int success = ((set_spatial_domain_fn*)domain_config->main_function)(domain_config, grid);
@@ -117,15 +112,16 @@ static void expand_scar(const char *input, const char *output, int n_rows, char*
     real_cpu dx, dy, dz;
     real_cpu min_x, max_x, min_y, max_y, min_z, max_z;
     real *f, *s, *n;
-    real_cpu transmurality, base_apex_heterogeneity, apicobasal;
-    uint32_t transmurality_labels;
+    real_cpu transmurality, transmurality_labels, base_apex;
+    real_cpu dense_sparse_tag, healthy, sf_IKs;
 
-    min_x=76803.6;
-    min_y=51176.3;
-    min_z=8740.9;
-    max_x=96803.6;
-    max_y=71176.3;
-    max_z=28740.9;
+    // Box bounds
+    min_x=67108;
+    min_y=52945;
+    min_z=14308;
+    max_x=77108;
+    max_y=62945;
+    max_z=24308;
 
     FILE *out = fopen(output, "w");
     FOR_EACH_CELL(grid) {
@@ -141,22 +137,20 @@ static void expand_scar(const char *input, const char *output, int n_rows, char*
 
             extra_data = (struct dti_mesh_info *)cell->mesh_extra_info;
             transmurality = extra_data->transmurality;
-            base_apex_heterogeneity = extra_data->base_apex_heterogeneity;
-            apicobasal = extra_data->apicobasal;
             transmurality_labels = extra_data->dti_transmurality_labels;
-            
-            // Change the FAST_ENDO tag to ENDO
-            if (transmurality_labels == 3) {
-                transmurality_labels = 0;
-            }
+            base_apex = extra_data->base_apex_heterogeneity;
+            dense_sparse_tag = extra_data->fast_endo;
+            healthy = 1.0;
+            sf_IKs = extra_data->sf_IKs;
 
             ignore_cell = cell->center.x < min_x || cell->center.x > max_x || cell->center.y < min_y || cell->center.y > max_y || cell->center.z < min_z || cell->center.z > max_z;
 
             if(!ignore_cell) {
-                fprintf(out, "%g,%g,%g,%g,%g,%g,%g,%g,%g,%u,%g,%g,%g,%g,%g,%g,%g,%g,%g\n", cell->center.x, cell->center.y, cell->center.z, dx, dy, dz, \
-                                                    transmurality, base_apex_heterogeneity, apicobasal, \
-                                                    transmurality_labels, \
-                                                    f[0], f[1], f[2], s[0], s[1], s[2], n[0], n[1], n[2]);
+                fprintf(out, "%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g\n", cell->center.x, cell->center.y, cell->center.z, dx, dy, dz, \
+                                                    transmurality, transmurality_labels, base_apex, \
+                                                    dense_sparse_tag, healthy, \
+                                                    f[0], f[1], f[2], s[0], s[1], s[2], n[0], n[1], n[2], \
+                                                    sf_IKs);
             }
 
         }
@@ -183,7 +177,7 @@ int main(int argc, char **argv) {
     }
 
     if(!output) {
-        output = "expanded_mesh.alg";
+        output = "clipped_mesh.alg";
     }
 
     expand_scar(input, output, options->n_rows, options->n_volumes);
diff --git a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.c b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.c
index bf32adbd..44dbe570 100644
--- a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.c
+++ b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.c
@@ -117,6 +117,7 @@ SET_ODE_INITIAL_CONDITIONS_CPU(set_model_initial_conditions_cpu) {
             sv[46] = 9.993734e-01;
             sv[47] = 0.000000e+00;
             sv[48] = 0.000000e+00;
+	        sv[49] = 0.000000e+00;
 
             // Default initial conditions for MID cell (from original Matlab script)
             //sv[0] = -8.953800e+01;
@@ -265,7 +266,7 @@ SOLVE_MODEL_ODES(solve_model_odes_cpu) {
     else {
         // Default: initialize all current modifiers
         for (uint32_t i = 0; i < num_extra_parameters; i++) {
-            if (i == 9)
+            if (i == 10)
                 extra_par[i] = 0.0;
             else
                 extra_par[i] = 1.0;
@@ -368,6 +369,7 @@ void solve_model_ode_cpu(real dt, real *sv, real stim_current, real transmuralit
     SOLVE_EQUATION_EULER_CPU(46);       // TmBlocked
     SOLVE_EQUATION_EULER_CPU(47);       // ZETAS
     SOLVE_EQUATION_EULER_CPU(48);       // ZETAW 
+    SOLVE_EQUATION_CONSTANT_CPU(49);	// Ta
 }
 
 void solve_forward_euler_cpu_adpt(real *sv, real stim_curr, real transmurality, real final_time, int sv_id, struct ode_solver *solver, real const *extra_params) {
@@ -583,6 +585,7 @@ void RHS_cpu(const real *sv, real *rDY_, real stim_current, real dt, real transm
     real TmBlocked = sv[46];
     real ZETAS = sv[47];
     real ZETAW = sv[48];
+    real TA = sv[49];
 
     #include "ToRORd_Land_mixed_endo_mid_epi.common.c"
 }
@@ -678,5 +681,7 @@ void RHS_RL_cpu(real *a_, real *b_, const real *sv, real *rDY_, real stim_curren
     real ZETAS = sv[47];
     real ZETAW = sv[48];
 
+    real TA = sv[49];
+
     #include "ToRORd_Land_mixed_endo_mid_epi_RL.common.c"
 }
diff --git a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.common.c b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.common.c
index 065a8c08..e371e395 100644
--- a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.common.c
+++ b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.common.c
@@ -828,3 +828,4 @@ rDY_[45] = dCa_TRPN;
 rDY_[46] = dTmBlocked;
 rDY_[47] = dZETAS;
 rDY_[48] = dZETAW;
+rDY_[49] = Ta;
diff --git a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.cu b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.cu
index 5b0f64c7..d048e054 100644
--- a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.cu
+++ b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.cu
@@ -57,12 +57,13 @@ __global__ void kernel_set_model_initial_conditions(real *sv, int num_volumes, s
             *((real * )((char *) sv + pitch * 46) + threadID) = 9.993734e-01;
             *((real * )((char *) sv + pitch * 47) + threadID) = 0.000000e+00;
             *((real * )((char *) sv + pitch * 48) + threadID) = 0.000000e+00;
+	        *((real * )((char *) sv + pitch * 49) + threadID) = 0.000000e+00;
         }
             
         if(use_adpt_dt) {
-            *((real *)((char *)sv + pitch * 49) + threadID) = min_dt; // dt
-            *((real *)((char *)sv + pitch * 50) + threadID) = 0.0;    // time_new
-            *((real *)((char *)sv + pitch * 51) + threadID) = 0.0;    // previous dt
+            *((real *)((char *)sv + pitch * 50) + threadID) = min_dt; // dt
+            *((real *)((char *)sv + pitch * 51) + threadID) = 0.0;    // time_new
+            *((real *)((char *)sv + pitch * 52) + threadID) = 0.0;    // previous dt
         }
     }
 }
@@ -230,7 +231,7 @@ extern "C" SOLVE_MODEL_ODES(solve_model_odes_gpu) {
     else {
         // Default: initialize all current modifiers
         for (uint32_t i = 0; i < num_extra_parameters; i++) {
-            if (i == 9)
+            if (i == 10)
                 extra_par[i] = 0.0;
             else 
                 extra_par[i] = 1.0;
@@ -334,6 +335,7 @@ __global__ void solve_gpu(real cur_time, real dt, real *sv, real *stim_currents,
                 SOLVE_EQUATION_EULER_GPU(46);       // TmBlocked
                 SOLVE_EQUATION_EULER_GPU(47);       // ZETAS
                 SOLVE_EQUATION_EULER_GPU(48);       // ZETAW
+		        SOLVE_EQUATION_CONSTANT_GPU(49);    // Ta
             }
         } else {
             solve_forward_euler_gpu_adpt(sv, stim_currents[threadID], 0.0, extra_params, cur_time + max_dt, sv_id, pitch, abstol,  reltol,  dt,  max_dt);
@@ -416,6 +418,7 @@ __global__ void solve_endo_mid_epi_gpu(real cur_time, real dt, real *sv, real *s
                 SOLVE_EQUATION_EULER_GPU(46);       // TmBlocked
                 SOLVE_EQUATION_EULER_GPU(47);       // ZETAS
                 SOLVE_EQUATION_EULER_GPU(48);       // ZETAW
+		        SOLVE_EQUATION_CONSTANT_GPU(49);    // Ta
             }
         } else {
             solve_forward_euler_gpu_adpt(sv, stim_currents[threadID], transmurality[threadID], extra_params, cur_time + max_dt, sv_id, pitch, abstol,  reltol,  dt,  max_dt);
@@ -645,6 +648,7 @@ inline __device__ void RHS_gpu(real *sv, real *rDY_, real stim_current, \
     real ZETAS;
     real ZETAW;
 
+    real TA;
     if (use_adpt_dt) {
         v = sv[0];
         nai = sv[1];
@@ -701,6 +705,7 @@ inline __device__ void RHS_gpu(real *sv, real *rDY_, real stim_current, \
         TmBlocked = sv[46];
         ZETAS = sv[47];
         ZETAW = sv[48];
+        TA = sv[49];
     } else {
         v = *((real *)((char *)sv + pitch * 0) + threadID_);
         nai = *((real *)((char *)sv + pitch * 1) + threadID_);
@@ -758,6 +763,8 @@ inline __device__ void RHS_gpu(real *sv, real *rDY_, real stim_current, \
         TmBlocked = (*((real *)((char *)sv + pitch * 46) + threadID_));
         ZETAS = (*((real *)((char *)sv + pitch * 47) + threadID_));
         ZETAW = (*((real *)((char *)sv + pitch * 48) + threadID_));
+
+	TA = (*((real *)((char *)sv + pitch * 49) + threadID_));
     }
 
     #include "ToRORd_Land_mixed_endo_mid_epi.common.c"
@@ -855,6 +862,8 @@ inline __device__ void RHS_RL_gpu(real *a_, real *b_, real *sv, real *rDY_, real
     real ZETAS;
     real ZETAW;
 
+    real TA;
+
     if (use_adpt_dt) {
         v = sv[0];
         nai = sv[1];
@@ -911,6 +920,7 @@ inline __device__ void RHS_RL_gpu(real *a_, real *b_, real *sv, real *rDY_, real
         TmBlocked = sv[46];
         ZETAS = sv[47];
         ZETAW = sv[48];
+	    TA = sv[49];
     } else {
         v = *((real *)((char *)sv + pitch * 0) + threadID_);
         nai = *((real *)((char *)sv + pitch * 1) + threadID_);
@@ -968,6 +978,7 @@ inline __device__ void RHS_RL_gpu(real *a_, real *b_, real *sv, real *rDY_, real
         TmBlocked = (*((real *)((char *)sv + pitch * 46) + threadID_));
         ZETAS = (*((real *)((char *)sv + pitch * 47) + threadID_));
         ZETAW = (*((real *)((char *)sv + pitch * 48) + threadID_));
+	    TA = (*((real *)((char *)sv + pitch * 49) + threadID_));
     }
 
     #include "ToRORd_Land_mixed_endo_mid_epi_RL.common.c"
diff --git a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.h b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.h
index a2a874e3..258780fa 100644
--- a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.h
+++ b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.h
@@ -10,7 +10,7 @@
 #include "../../extra_data_library/helper_functions.h"
 #include <stdlib.h>
 
-#define NEQ 49
+#define NEQ 50
 #define INITIAL_V (-8.863699e+01)
 
 #define ENDO 0.0
@@ -18,7 +18,7 @@
 #define EPI  2.0
 
 // CPU macros
-#define SOLVE_EQUATION_CONSTANT_CPU(id) sv[id] = rY[id]
+#define SOLVE_EQUATION_CONSTANT_CPU(id) sv[id] = rDY[id]
 
 #define SOLVE_EQUATION_EULER_CPU(id) sv[id] = dt * rDY[id] + rY[id]
 
@@ -49,7 +49,7 @@
                                                     greatestError = (auxError > greatestError) ? auxError : greatestError
 
 // GPU macros
-#define SOLVE_EQUATION_CONSTANT_GPU(id) *((real *)((char *)sv + pitch * id) + sv_id) = *((real *)((char *)sv + pitch * id) + sv_id)
+#define SOLVE_EQUATION_CONSTANT_GPU(id) *((real *)((char *)sv + pitch * id) + sv_id) = rDY[id]
 
 #define SOLVE_EQUATION_EULER_GPU(id) *((real *)((char *)sv + pitch * id) + sv_id) = *((real *)((char *)sv + pitch * id) + sv_id) + dt * rDY[id]
 
diff --git a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi_RL.common.c b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi_RL.common.c
index 66cb582a..331c446f 100644
--- a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi_RL.common.c
+++ b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi_RL.common.c
@@ -884,3 +884,4 @@ rDY_[45] = dCa_TRPN;
 rDY_[46] = dTmBlocked;
 rDY_[47] = dZETAS;
 rDY_[48] = dZETAW;
+rDY_[49] = Ta;
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
index ed62d253..a1591e30 100644
--- 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
@@ -68,50 +68,50 @@ SET_ODE_INITIAL_CONDITIONS_CPU(set_model_initial_conditions_cpu) {
         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;
+            // Default initial condition from Matlab (endocardium cell)
+            sv[0]  = -88.6369922306458;
+            sv[1]  = 11.8973412949238;
+            sv[2]  = 11.897661047085;
+            sv[3]  = 141.234464714982;
+            sv[4]  = 141.234423402713;
+            sv[5]  = 7.26747296460659e-05;
+            sv[6]  = 6.33786975780735e-05;
+            sv[7]  = 1.5326530637197;
+            sv[8]  = 1.53394579180493;
+            sv[9]  = 0.000828007761976018;
+            sv[10] = 0.666527193684116;
+            sv[11] = 0.826020806005678;
+            sv[12] = 0.826055985895856;
+            sv[13] = 0.825850881115628;
+            sv[14] = 0.000166868626513013;
+            sv[15] = 0.522830604669169;
+            sv[16] = 0.285969584294187;
+            sv[17] = 0.000959137028030184;
+            sv[18] = 0.999601150012565;
+            sv[19] = 0.5934016398361;
+            sv[20] = 0.000488696137242056;
+            sv[21] = 0.999601147267179;
+            sv[22] = 0.654668660159696;
+            sv[23] = 9.50007519781516e-32;
+            sv[24] = 0.999999992317577;
+            sv[25] = 0.939258048397962;
+            sv[26] = 0.999999992317557;
+            sv[27] = 0.999898379647465;
+            sv[28] = 0.99997825156004;
+            sv[29] = 0.000444816183420527;
+            sv[30] = 0.000755072490632667;
+            sv[31] = 0.999999992318446;
+            sv[32] = 0.999999992318445;
+            sv[33] = 0.24240468344952;
+            sv[34] = 0.000179537726989804;
+            sv[35] = -6.88308558109975e-25;
+            sv[36] = 0.0111749845355653;
+            sv[37] = 0.998036620213316;
+            sv[38] = 0.000858801779013532;
+            sv[39] = 0.000709744678350176;
+            sv[40] = 0.000381261722195702;
+            sv[41] = 1.35711566929992e-05;
+            sv[42] = 2.30252452954649e-23;
         }
     }        
 }
@@ -160,7 +160,7 @@ SOLVE_MODEL_ODES(solve_model_odes_cpu) {
     else {
         // Default: initialize all current modifiers
         for (uint32_t i = 0; i < num_extra_parameters; i++) {
-            if (i == 9)
+            if (i == 10)
                 extra_par[i] = 0.0;
             else 
                 extra_par[i] = 1.0;
@@ -212,47 +212,47 @@ void solve_model_ode_cpu(real dt, real *sv, real stim_current, real transmuralit
     //  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(1);        // nai
+    SOLVE_EQUATION_EULER_CPU(2);        // nass
+    SOLVE_EQUATION_EULER_CPU(3);        // ki   
+    SOLVE_EQUATION_EULER_CPU(4);        // kss
+    SOLVE_EQUATION_EULER_CPU(5);        // cai
+    SOLVE_EQUATION_EULER_CPU(6);        // cass 
     SOLVE_EQUATION_EULER_CPU(7);        // cansr
     SOLVE_EQUATION_EULER_CPU(8);        // cajsr
-    SOLVE_EQUATION_EULER_CPU(9);        // cai
-    SOLVE_EQUATION_RUSH_LARSEN_CPU(10); // m
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(9);  // m
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(10); // hp
     SOLVE_EQUATION_RUSH_LARSEN_CPU(11); // h
     SOLVE_EQUATION_RUSH_LARSEN_CPU(12); // j
-    SOLVE_EQUATION_RUSH_LARSEN_CPU(13); // 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(13); // jp
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(14); // mL
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(15); // hL
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(16); // hLp
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(17); // a
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(18); // iF
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(19); // iS
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(20); // ap
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(21); // iFp
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(22); // iSp
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(23); // d
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(24); // ff
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(25); // fs
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(26); // fcaf
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(27); // fcas
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(28); // jca
+    SOLVE_EQUATION_EULER_CPU(29);       // nca
+    SOLVE_EQUATION_EULER_CPU(30);       // nca_i
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(31); // ffp
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(32); // fcafp
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(33); // xs1
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(34); // xs2
+    SOLVE_EQUATION_RUSH_LARSEN_CPU(35); // Jrel_np
+    SOLVE_EQUATION_EULER_CPU(36);        // CaMKt 
+    SOLVE_EQUATION_EULER_CPU(37);       // ikr_c0
+    SOLVE_EQUATION_EULER_CPU(38);       // ikr_c1
+    SOLVE_EQUATION_EULER_CPU(39);       // ikr_c2
+    SOLVE_EQUATION_EULER_CPU(40);       // ikr_o
+    SOLVE_EQUATION_EULER_CPU(41);       // ikr_i
     SOLVE_EQUATION_RUSH_LARSEN_CPU(42); // Jrel_p
 }
 
@@ -410,49 +410,49 @@ void RHS_cpu(const real *sv, real *rDY_, real stim_current, real dt, real transm
     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];
+    real v       = sv[0];
+    real nai     = sv[1];
+    real nass    = sv[2];
+    real ki      = sv[3];
+    real kss     = sv[4];
+    real cai     = sv[5];
+    real cass    = sv[6];
+    real cansr   = sv[7];
+    real cajsr   = sv[8];
+    real m       = sv[9];
+    real hp      = sv[10];
+    real h       = sv[11];
+    real j       = sv[12];
+    real jp      = sv[13];
+    real mL      = sv[14];
+    real hL      = sv[15];
+    real hLp     = sv[16];
+    real a       = sv[17];
+    real iF      = sv[18];
+    real iS      = sv[19];
+    real ap      = sv[20];
+    real iFp     = sv[21];
+    real iSp     = sv[22];
+    real d       = sv[23];
+    real ff      = sv[24];
+    real fs      = sv[25];
+    real fcaf    = sv[26];
+    real fcas    = sv[27];
+    real jca     = sv[28];
+    real nca     = sv[29];
+    real nca_i   = sv[30];
+    real ffp     = sv[31];
+    real fcafp   = sv[32];
+    real xs1     = sv[33];
+    real xs2     = sv[34];
+    real Jrel_np = sv[35];
+    real CaMKt   = sv[36];
+    real ikr_c0  = sv[37];
+    real ikr_c1  = sv[38];
+    real ikr_c2  = sv[39];
+    real ikr_o   = sv[40];
+    real ikr_i   = sv[41];
+    real Jrel_p  = sv[42];
 
     #include "ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.common.c"
 }
@@ -488,49 +488,49 @@ void RHS_RL_cpu(real *a_, real *b_, const real *sv, real *rDY_, real stim_curren
     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];
+    real v       = sv[0];
+    real nai     = sv[1];
+    real nass    = sv[2];
+    real ki      = sv[3];
+    real kss     = sv[4];
+    real cai     = sv[5];
+    real cass    = sv[6];
+    real cansr   = sv[7];
+    real cajsr   = sv[8];
+    real m       = sv[9];
+    real hp      = sv[10];
+    real h       = sv[11];
+    real j       = sv[12];
+    real jp      = sv[13];
+    real mL      = sv[14];
+    real hL      = sv[15];
+    real hLp     = sv[16];
+    real a       = sv[17];
+    real iF      = sv[18];
+    real iS      = sv[19];
+    real ap      = sv[20];
+    real iFp     = sv[21];
+    real iSp     = sv[22];
+    real d       = sv[23];
+    real ff      = sv[24];
+    real fs      = sv[25];
+    real fcaf    = sv[26];
+    real fcas    = sv[27];
+    real jca     = sv[28];
+    real nca     = sv[29];
+    real nca_i   = sv[30];
+    real ffp     = sv[31];
+    real fcafp   = sv[32];
+    real xs1     = sv[33];
+    real xs2     = sv[34];
+    real Jrel_np = sv[35];
+    real CaMKt   = sv[36];
+    real ikr_c0  = sv[37];
+    real ikr_c1  = sv[38];
+    real ikr_c2  = sv[39];
+    real ikr_o   = sv[40];
+    real ikr_i   = sv[41];
+    real Jrel_p  = sv[42];
 
     #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
index f20b04be..83e314e6 100644
--- 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
@@ -225,7 +225,7 @@ real PCa=8.3757e-05 * ICaL_Multiplier;
 if (celltype==EPI)
     PCa=PCa*1.2;
 else if (celltype==MID)
-    PCa=PCa*2;
+    PCa=PCa*1.8;
 
 real PCap=1.1*PCa;
 real PCaNa=0.00125*PCa;
@@ -300,7 +300,7 @@ real dc2 = c1 * alpha1 + o*beta2 + I*betaItoC2 - c2 * (beta1 + alpha2 + alphac2T
 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. 
+real GKr = 0.0321 * sqrt(ko/5) * IKr_Multiplier * 0.5; // 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)
@@ -593,7 +593,8 @@ real I_katp = (fkatp * gkatp * akik * bkik * (v - EK));
 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
+//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
+real dv=-(INa+INaL+Ito+ICaL+ICaNa+ICaK+IKr+IKs+IK1+INaCa_i+INaCa_ss+INaK+INab+IKb+IpCa+ICab+I_ClCa+I_Clbk+Istim);    // Euler
 
 // calculate diffusion fluxes
 real JdiffNa=(nass-nai)/2.0;
@@ -634,45 +635,45 @@ 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_[1]  = dnai;
+rDY_[2]  = dnass;
+rDY_[3]  = dki;
+rDY_[4]  = dkss;
+rDY_[5]  = dcai;
+rDY_[6]  = dcass;
 rDY_[7]  = dcansr;
 rDY_[8]  = dcajsr;
-rDY_[9]  = dcai;
-rDY_[10] = dm;
+rDY_[9]  = dm;
+rDY_[10] = dhp;
 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_[13] = djp;
+rDY_[14] = dmL;
+rDY_[15] = dhL;
+rDY_[16] = dhLp;
+rDY_[17] = da;
+rDY_[18] = diF;
+rDY_[19] = diS;
+rDY_[20] = dap;
+rDY_[21] = diFp;
+rDY_[22] = diSp;
+rDY_[23] = dd;
+rDY_[24] = dff;
+rDY_[25] = dfs;
+rDY_[26] = dfcaf;
+rDY_[27] = dfcas;
+rDY_[28] = djca;
+rDY_[29] = dnca;
+rDY_[30] = dnca_i;
+rDY_[31] = dffp;
+rDY_[32] = dfcafp;
+rDY_[33] = dxs1;
+rDY_[34] = dxs2;
+rDY_[35] = dJrelnp;
+rDY_[36] = dCaMKt;
+rDY_[37] = dc0;
+rDY_[38] = dc1;
+rDY_[39] = dc2;
+rDY_[40] = delta_o;
+rDY_[41] = di;
 rDY_[42] = dJrelp;
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
index 41bb3213..79301d0d 100644
--- 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
@@ -5,53 +5,52 @@
 __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;
+            // Default initial condition from Matlab (endocardium cell)
+            *((real * )((char *) sv + pitch * 0) + threadID) = -88.6369922306458;
+            *((real * )((char *) sv + pitch * 1) + threadID) = 11.8973412949238;
+            *((real * )((char *) sv + pitch * 2) + threadID) = 11.897661047085;
+            *((real * )((char *) sv + pitch * 3) + threadID) = 141.234464714982;
+            *((real * )((char *) sv + pitch * 4) + threadID) = 141.234423402713;
+            *((real * )((char *) sv + pitch * 5) + threadID) = 7.26747296460659e-05;
+            *((real * )((char *) sv + pitch * 6) + threadID) = 6.33786975780735e-05;
+            *((real * )((char *) sv + pitch * 7) + threadID) = 1.5326530637197;
+            *((real * )((char *) sv + pitch * 8) + threadID) = 1.53394579180493;
+            *((real * )((char *) sv + pitch * 9) + threadID) = 0.000828007761976018;
+            *((real * )((char *) sv + pitch * 10) + threadID) = 0.666527193684116;
+            *((real * )((char *) sv + pitch * 11) + threadID) = 0.826020806005678;
+            *((real * )((char *) sv + pitch * 12) + threadID) = 0.826055985895856;
+            *((real * )((char *) sv + pitch * 13) + threadID) = 0.825850881115628;
+            *((real * )((char *) sv + pitch * 14) + threadID) = 0.000166868626513013;
+            *((real * )((char *) sv + pitch * 15) + threadID) = 0.522830604669169;
+            *((real * )((char *) sv + pitch * 16) + threadID) = 0.285969584294187;
+            *((real * )((char *) sv + pitch * 17) + threadID) = 0.000959137028030184;
+            *((real * )((char *) sv + pitch * 18) + threadID) = 0.999601150012565;
+            *((real * )((char *) sv + pitch * 19) + threadID) = 0.5934016398361;
+            *((real * )((char *) sv + pitch * 20) + threadID) = 0.000488696137242056;
+            *((real * )((char *) sv + pitch * 21) + threadID) = 0.999601147267179;
+            *((real * )((char *) sv + pitch * 22) + threadID) = 0.654668660159696;
+            *((real * )((char *) sv + pitch * 23) + threadID) = 9.50007519781516e-32;
+            *((real * )((char *) sv + pitch * 24) + threadID) = 0.999999992317577;
+            *((real * )((char *) sv + pitch * 25) + threadID) = 0.939258048397962;
+            *((real * )((char *) sv + pitch * 26) + threadID) = 0.999999992317557;
+            *((real * )((char *) sv + pitch * 27) + threadID) = 0.999898379647465;
+            *((real * )((char *) sv + pitch * 28) + threadID) = 0.99997825156004;
+            *((real * )((char *) sv + pitch * 29) + threadID) = 0.000444816183420527;
+            *((real * )((char *) sv + pitch * 30) + threadID) = 0.000755072490632667;
+            *((real * )((char *) sv + pitch * 31) + threadID) = 0.999999992318446;
+            *((real * )((char *) sv + pitch * 32) + threadID) = 0.999999992318445;
+            *((real * )((char *) sv + pitch * 33) + threadID) = 0.24240468344952;
+            *((real * )((char *) sv + pitch * 34) + threadID) = 0.000179537726989804;
+            *((real * )((char *) sv + pitch * 35) + threadID) = -6.88308558109975e-25;
+            *((real * )((char *) sv + pitch * 36) + threadID) = 0.0111749845355653;
+            *((real * )((char *) sv + pitch * 37) + threadID) = 0.998036620213316;
+            *((real * )((char *) sv + pitch * 38) + threadID) = 0.000858801779013532;
+            *((real * )((char *) sv + pitch * 39) + threadID) = 0.000709744678350176;
+            *((real * )((char *) sv + pitch * 40) + threadID) = 0.000381261722195702;
+            *((real * )((char *) sv + pitch * 41) + threadID) = 1.35711566929992e-05;
+            *((real * )((char *) sv + pitch * 42) + threadID) = 2.30252452954649e-23;
         }
             
         if(use_adpt_dt) {
@@ -240,7 +239,7 @@ extern "C" SOLVE_MODEL_ODES(solve_model_odes_gpu) {
 
         // Default: initialize all current modifiers
         for (uint32_t i = 0; i < num_extra_parameters; i++) {
-            if (i == 9)
+            if (i == 10)
                 extra_par[i] = 0.0;
             else 
                 extra_par[i] = 1.0;
@@ -294,47 +293,47 @@ __global__ void solve_endo_mid_epi_gpu(real cur_time, real dt, real *sv, real *s
                 //  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(1);        // nai
+                SOLVE_EQUATION_EULER_GPU(2);        // nass
+                SOLVE_EQUATION_EULER_GPU(3);        // ki   
+                SOLVE_EQUATION_EULER_GPU(4);        // kss
+                SOLVE_EQUATION_EULER_GPU(5);        // cai
+                SOLVE_EQUATION_EULER_GPU(6);        // cass 
                 SOLVE_EQUATION_EULER_GPU(7);        // cansr
                 SOLVE_EQUATION_EULER_GPU(8);        // cajsr
-                SOLVE_EQUATION_EULER_GPU(9);        // cai
-                SOLVE_EQUATION_RUSH_LARSEN_GPU(10); // m
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(9);  // m
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(10); // hp
                 SOLVE_EQUATION_RUSH_LARSEN_GPU(11); // h
                 SOLVE_EQUATION_RUSH_LARSEN_GPU(12); // j
-                SOLVE_EQUATION_RUSH_LARSEN_GPU(13); // 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(13); // jp
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(14); // mL
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(15); // hL
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(16); // hLp
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(17); // a
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(18); // iF
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(19); // iS
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(20); // ap
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(21); // iFp
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(22); // iSp
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(23); // d
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(24); // ff
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(25); // fs
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(26); // fcaf
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(27); // fcas
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(28); // jca
+                SOLVE_EQUATION_EULER_GPU(29);       // nca
+                SOLVE_EQUATION_EULER_GPU(30);       // nca_i
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(31); // ffp
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(32); // fcafp
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(33); // xs1
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(34); // xs2
+                SOLVE_EQUATION_RUSH_LARSEN_GPU(35); // Jrel_np
+                SOLVE_EQUATION_EULER_GPU(36);        // CaMKt 
+                SOLVE_EQUATION_EULER_GPU(37);       // ikr_c0
+                SOLVE_EQUATION_EULER_GPU(38);       // ikr_c1
+                SOLVE_EQUATION_EULER_GPU(39);       // ikr_c2
+                SOLVE_EQUATION_EULER_GPU(40);       // ikr_o
+                SOLVE_EQUATION_EULER_GPU(41);       // ikr_i
                 SOLVE_EQUATION_RUSH_LARSEN_GPU(42); // Jrel_p
             }
         } else {
@@ -549,93 +548,93 @@ inline __device__ void RHS_gpu(real *sv, real *rDY_, real stim_current, real tra
     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];
+        v       = sv[0];
+        nai     = sv[1];
+        nass    = sv[2];
+        ki      = sv[3];
+        kss     = sv[4];
+        cai     = sv[5];
+        cass    = sv[6];
+        cansr   = sv[7];
+        cajsr   = sv[8];
+        m       = sv[9];
+        hp      = sv[10];
+        h       = sv[11];
+        j       = sv[12];
+        jp      = sv[13];
+        mL      = sv[14];
+        hL      = sv[15];
+        hLp     = sv[16];
+        a       = sv[17];
+        iF      = sv[18];
+        iS      = sv[19];
+        ap      = sv[20];
+        iFp     = sv[21];
+        iSp     = sv[22];
+        d       = sv[23];
+        ff      = sv[24];
+        fs      = sv[25];
+        fcaf    = sv[26];
+        fcas    = sv[27];
+        jca     = sv[28];
+        nca     = sv[29];
+        nca_i   = sv[30];
+        ffp     = sv[31];
+        fcafp   = sv[32];
+        xs1     = sv[33];
+        xs2     = sv[34];
+        Jrel_np = sv[35];
+        CaMKt   = sv[36];
+        ikr_c0  = sv[37];
+        ikr_c1  = sv[38];
+        ikr_c2  = sv[39];
+        ikr_o   = sv[40];
+        ikr_i   = 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_);
+        v       = *((real *)((char *)sv + pitch * 0) + threadID_);
+        nai     = *((real *)((char *)sv + pitch * 1) + threadID_);
+        nass    = *((real *)((char *)sv + pitch * 2) + threadID_);
+        ki      = *((real *)((char *)sv + pitch * 3) + threadID_);
+        kss     = *((real *)((char *)sv + pitch * 4) + threadID_);
+        cai     = *((real *)((char *)sv + pitch * 5) + threadID_);
+        cass    = *((real *)((char *)sv + pitch * 6) + threadID_);
+        cansr   = *((real *)((char *)sv + pitch * 7) + threadID_);
+        cajsr   = *((real *)((char *)sv + pitch * 8) + threadID_);
+        m       = *((real *)((char *)sv + pitch * 9) + threadID_);
+        hp      = *((real *)((char *)sv + pitch * 10) + threadID_);
+        h       = *((real *)((char *)sv + pitch * 11) + threadID_);
+        j       = *((real *)((char *)sv + pitch * 12) + threadID_);
+        jp      = *((real *)((char *)sv + pitch * 13) + threadID_);
+        mL      = *((real *)((char *)sv + pitch * 14) + threadID_);
+        hL      = *((real *)((char *)sv + pitch * 15) + threadID_);
+        hLp     = *((real *)((char *)sv + pitch * 16) + threadID_);
+        a       = *((real *)((char *)sv + pitch * 17) + threadID_);
+        iF      = *((real *)((char *)sv + pitch * 18) + threadID_);
+        iS      = *((real *)((char *)sv + pitch * 19) + threadID_);
+        ap      = *((real *)((char *)sv + pitch * 20) + threadID_);
+        iFp     = *((real *)((char *)sv + pitch * 21) + threadID_);
+        iSp     = *((real *)((char *)sv + pitch * 22) + threadID_);
+        d       = *((real *)((char *)sv + pitch * 23) + threadID_);
+        ff      = *((real *)((char *)sv + pitch * 24) + threadID_);
+        fs      = *((real *)((char *)sv + pitch * 25) + threadID_);
+        fcaf    = *((real *)((char *)sv + pitch * 26) + threadID_);
+        fcas    = *((real *)((char *)sv + pitch * 27) + threadID_);
+        jca     = *((real *)((char *)sv + pitch * 28) + threadID_);
+        nca     = *((real *)((char *)sv + pitch * 29) + threadID_);
+        nca_i   = *((real *)((char *)sv + pitch * 30) + threadID_);
+        ffp     = *((real *)((char *)sv + pitch * 31) + threadID_);
+        fcafp   = *((real *)((char *)sv + pitch * 32) + threadID_);
+        xs1     = *((real *)((char *)sv + pitch * 33) + threadID_);
+        xs2     = *((real *)((char *)sv + pitch * 34) + threadID_);
+        Jrel_np = *((real *)((char *)sv + pitch * 35) + threadID_);
+        CaMKt   = *((real *)((char *)sv + pitch * 36) + threadID_);
+        ikr_c0  = *((real *)((char *)sv + pitch * 37) + threadID_);
+        ikr_c1  = *((real *)((char *)sv + pitch * 38) + threadID_);
+        ikr_c2  = *((real *)((char *)sv + pitch * 39) + threadID_);
+        ikr_o   = *((real *)((char *)sv + pitch * 40) + threadID_);
+        ikr_i   = *((real *)((char *)sv + pitch * 41) + threadID_);
+        Jrel_p  = *((real *)((char *)sv + pitch * 42) + threadID_);
     }
 
     #include "ToRORd_fkatp_mixed_endo_mid_epi_GKsGKrtjca_adjustments.common.c"
@@ -718,93 +717,93 @@ inline __device__ void RHS_RL_gpu(real *a_, real *b_, real *sv, real *rDY_, real
     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];
+        v       = sv[0];
+        nai     = sv[1];
+        nass    = sv[2];
+        ki      = sv[3];
+        kss     = sv[4];
+        cai     = sv[5];
+        cass    = sv[6];
+        cansr   = sv[7];
+        cajsr   = sv[8];
+        m       = sv[9];
+        hp      = sv[10];
+        h       = sv[11];
+        j       = sv[12];
+        jp      = sv[13];
+        mL      = sv[14];
+        hL      = sv[15];
+        hLp     = sv[16];
+        a       = sv[17];
+        iF      = sv[18];
+        iS      = sv[19];
+        ap      = sv[20];
+        iFp     = sv[21];
+        iSp     = sv[22];
+        d       = sv[23];
+        ff      = sv[24];
+        fs      = sv[25];
+        fcaf    = sv[26];
+        fcas    = sv[27];
+        jca     = sv[28];
+        nca     = sv[29];
+        nca_i   = sv[30];
+        ffp     = sv[31];
+        fcafp   = sv[32];
+        xs1     = sv[33];
+        xs2     = sv[34];
+        Jrel_np = sv[35];
+        CaMKt   = sv[36];
+        ikr_c0  = sv[37];
+        ikr_c1  = sv[38];
+        ikr_c2  = sv[39];
+        ikr_o   = sv[40];
+        ikr_i   = 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_);
+        v       = *((real *)((char *)sv + pitch * 0) + threadID_);
+        nai     = *((real *)((char *)sv + pitch * 1) + threadID_);
+        nass    = *((real *)((char *)sv + pitch * 2) + threadID_);
+        ki      = *((real *)((char *)sv + pitch * 3) + threadID_);
+        kss     = *((real *)((char *)sv + pitch * 4) + threadID_);
+        cai     = *((real *)((char *)sv + pitch * 5) + threadID_);
+        cass    = *((real *)((char *)sv + pitch * 6) + threadID_);
+        cansr   = *((real *)((char *)sv + pitch * 7) + threadID_);
+        cajsr   = *((real *)((char *)sv + pitch * 8) + threadID_);
+        m       = *((real *)((char *)sv + pitch * 9) + threadID_);
+        hp      = *((real *)((char *)sv + pitch * 10) + threadID_);
+        h       = *((real *)((char *)sv + pitch * 11) + threadID_);
+        j       = *((real *)((char *)sv + pitch * 12) + threadID_);
+        jp      = *((real *)((char *)sv + pitch * 13) + threadID_);
+        mL      = *((real *)((char *)sv + pitch * 14) + threadID_);
+        hL      = *((real *)((char *)sv + pitch * 15) + threadID_);
+        hLp     = *((real *)((char *)sv + pitch * 16) + threadID_);
+        a       = *((real *)((char *)sv + pitch * 17) + threadID_);
+        iF      = *((real *)((char *)sv + pitch * 18) + threadID_);
+        iS      = *((real *)((char *)sv + pitch * 19) + threadID_);
+        ap      = *((real *)((char *)sv + pitch * 20) + threadID_);
+        iFp     = *((real *)((char *)sv + pitch * 21) + threadID_);
+        iSp     = *((real *)((char *)sv + pitch * 22) + threadID_);
+        d       = *((real *)((char *)sv + pitch * 23) + threadID_);
+        ff      = *((real *)((char *)sv + pitch * 24) + threadID_);
+        fs      = *((real *)((char *)sv + pitch * 25) + threadID_);
+        fcaf    = *((real *)((char *)sv + pitch * 26) + threadID_);
+        fcas    = *((real *)((char *)sv + pitch * 27) + threadID_);
+        jca     = *((real *)((char *)sv + pitch * 28) + threadID_);
+        nca     = *((real *)((char *)sv + pitch * 29) + threadID_);
+        nca_i   = *((real *)((char *)sv + pitch * 30) + threadID_);
+        ffp     = *((real *)((char *)sv + pitch * 31) + threadID_);
+        fcafp   = *((real *)((char *)sv + pitch * 32) + threadID_);
+        xs1     = *((real *)((char *)sv + pitch * 33) + threadID_);
+        xs2     = *((real *)((char *)sv + pitch * 34) + threadID_);
+        Jrel_np = *((real *)((char *)sv + pitch * 35) + threadID_);
+        CaMKt   = *((real *)((char *)sv + pitch * 36) + threadID_);
+        ikr_c0  = *((real *)((char *)sv + pitch * 37) + threadID_);
+        ikr_c1  = *((real *)((char *)sv + pitch * 38) + threadID_);
+        ikr_c2  = *((real *)((char *)sv + pitch * 39) + threadID_);
+        ikr_o   = *((real *)((char *)sv + pitch * 40) + threadID_);
+        ikr_i   = *((real *)((char *)sv + pitch * 41) + threadID_);
+        Jrel_p  = *((real *)((char *)sv + pitch * 42) + threadID_);
     }
 
     #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
index 33b0a174..53cdb76d 100644
--- 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
@@ -8,7 +8,7 @@
 #include "../../extra_data_library/helper_functions.h"
 
 #define NEQ 43
-#define INITIAL_V (-8.890585e+01)
+#define INITIAL_V (-88.6369922306458)
 
 #define ENDO 0.0
 #define MID  1.0
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
index 1efd0ea0..4dd686d1 100644
--- 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
@@ -225,7 +225,7 @@ real PCa=8.3757e-05 * ICaL_Multiplier;
 if (celltype==EPI)
     PCa=PCa*1.2;
 else if (celltype==MID)
-    PCa=PCa*2;
+    PCa=PCa*1.8;
 
 real PCap=1.1*PCa;
 real PCaNa=0.00125*PCa;
@@ -300,7 +300,7 @@ real dc2 = c1 * alpha1 + o*beta2 + I*betaItoC2 - c2 * (beta1 + alpha2 + alphac2T
 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. 
+real GKr = 0.0321 * sqrt(ko/5) * IKr_Multiplier * 0.5; // 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)
@@ -594,6 +594,7 @@ 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
+//real dv=-(INa+INaL+Ito+ICaL+ICaNa+ICaK+IKr+IKs+IK1+INaCa_i+INaCa_ss+INaK+INab+IKb+IpCa+ICab+I_ClCa+I_Clbk+Istim);    // Euler
 
 // calculate diffusion fluxes
 real JdiffNa=(nass-nai)/2.0;
@@ -633,102 +634,102 @@ 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_[9 ] = -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_[10] = -1.0 / tauh;
+a_[13] = -1.0 / taujp;
+a_[14] = -1.0 / tmL;
+a_[15] = -1.0 / thL;
+a_[16] = -1.0 / thLp;
+a_[17] = -1.0 / ta;
+a_[18] = -1.0 / tiF;
+a_[19] = -1.0 / tiS;
+a_[20] = -1.0 / ta;
+a_[21] = -1.0 / tiFp;
+a_[22] = -1.0 / tiSp;
+a_[23] = -1.0 / td;
+a_[24] = -1.0 / tff;
+a_[25] = -1.0 / tfs;
+a_[26] = -1.0 / tfcaf;
+a_[27] = -1.0 / tfcas;
+a_[28] = -1.0 / tjca;
+a_[31] = -1.0 / tffp;
+a_[32] = -1.0 / tfcafp;
+a_[33] = -1.0 / txs1;
+a_[34] = -1.0 / txs2;
+a_[35] = -1.0 / tau_rel;
 a_[42] = -1.0 / tau_relp;    
 
 // Compute 'b' coefficients for the Hodkin-Huxley variables
-b_[10] = mss / tm;
+b_[9 ] = 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_[10] = hssp / tauh;
+b_[13] = jss / taujp;
+b_[14] = mLss / tmL;
+b_[15] = hLss / thL;
+b_[16] = hLssp / thLp;
+b_[17] = ass / ta;
+b_[18] = iss / tiF;
+b_[19] = iss / tiS;
+b_[20] = assp / ta;
+b_[21] = iss / tiFp;
+b_[22] = iss / tiSp;
+b_[23] = dss / td;
+b_[24] = fss / tff;
+b_[25] = fss / tfs;
+b_[26] = fcass / tfcaf;
+b_[27] = fcass / tfcas;
+b_[28] = jcass / tjca;
+b_[31] = fss / tffp;
+b_[32] = fcass / tfcafp;
+b_[33] = xs1ss / txs1;
+b_[34] = xs2ss / txs2;
+b_[35] = Jrel_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_[1]  = dnai;
+rDY_[2]  = dnass;
+rDY_[3]  = dki;
+rDY_[4]  = dkss;
+rDY_[5]  = dcai;
+rDY_[6]  = dcass;
 rDY_[7]  = dcansr;
 rDY_[8]  = dcajsr;
-rDY_[9]  = dcai;
-rDY_[10] = dm;
+rDY_[9]  = dm;
+rDY_[10] = dhp;
 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_[13] = djp;
+rDY_[14] = dmL;
+rDY_[15] = dhL;
+rDY_[16] = dhLp;
+rDY_[17] = da;
+rDY_[18] = diF;
+rDY_[19] = diS;
+rDY_[20] = dap;
+rDY_[21] = diFp;
+rDY_[22] = diSp;
+rDY_[23] = dd;
+rDY_[24] = dff;
+rDY_[25] = dfs;
+rDY_[26] = dfcaf;
+rDY_[27] = dfcas;
+rDY_[28] = djca;
+rDY_[29] = dnca;
+rDY_[30] = dnca_i;
+rDY_[31] = dffp;
+rDY_[32] = dfcafp;
+rDY_[33] = dxs1;
+rDY_[34] = dxs2;
+rDY_[35] = dJrelnp;
+rDY_[36] = dCaMKt;
+rDY_[37] = dc0;
+rDY_[38] = dc1;
+rDY_[39] = dc2;
+rDY_[40] = delta_o;
+rDY_[41] = di;
 rDY_[42] = dJrelp;
diff --git a/src/models_library/minimal_model/build.sh b/src/models_library/minimal_model/build.sh
new file mode 100644
index 00000000..f02ddf1b
--- /dev/null
+++ b/src/models_library/minimal_model/build.sh
@@ -0,0 +1,6 @@
+############## minimal_model ##############################
+MODEL_FILE_CPU="minimal_model.c"
+MODEL_FILE_GPU="minimal_model.cu"
+COMMON_HEADERS="minimal_model.h"
+
+COMPILE_MODEL_LIB "minimal_model" "$MODEL_FILE_CPU" "$MODEL_FILE_GPU" "$COMMON_HEADERS"
diff --git a/src/models_library/minimal_model/minimal_model.c b/src/models_library/minimal_model/minimal_model.c
new file mode 100644
index 00000000..02a696c2
--- /dev/null
+++ b/src/models_library/minimal_model/minimal_model.c
@@ -0,0 +1,185 @@
+#include <assert.h>
+#include <stdlib.h>
+#include <math.h>
+#include "minimal_model.h"
+#include <stdio.h>
+
+GET_CELL_MODEL_DATA(init_cell_model_data) {
+
+   assert(cell_model);
+
+    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 Minimal Model Mixed CPU model\n");
+
+    uint32_t num_cells = solver->original_num_cells;
+    solver->sv = (real*)malloc(NEQ*num_cells*sizeof(real));
+
+    OMP(parallel for)
+        for(uint32_t i = 0; i < num_cells; i++) {
+            real *sv = &solver->sv[i * NEQ];
+            sv[0] = 0.0f; //u
+            sv[1] = 1.0f; //v
+            sv[2] = 1.0f; //w
+            sv[3] = 0.0f; //s
+        }
+}
+
+SOLVE_MODEL_ODES(solve_model_odes_cpu) {
+
+    uint32_t sv_id;
+
+    uint32_t num_cells_to_solve = ode_solver->num_cells_to_solve;
+    size_t num_cells_to_solve_size = num_cells_to_solve*sizeof(uint32_t);
+    uint32_t * cells_to_solve = ode_solver->cells_to_solve;
+    real *sv = ode_solver->sv;
+    real dt = ode_solver->min_dt;
+    uint32_t num_steps = ode_solver->num_steps;
+
+    int i;
+
+    // Mapping array
+    uint32_t *mapping = NULL;
+    
+    // Extra data function will tag the cells in the grid
+    if (ode_solver->ode_extra_data) {
+        mapping = (uint32_t*)ode_solver->ode_extra_data; 
+    
+        OMP(parallel for private(sv_id))
+        for (i = 0; i < num_cells_to_solve; i++) {
+            if(cells_to_solve)
+                sv_id = cells_to_solve[i];
+            else
+                sv_id = i;
+
+            for (int j = 0; j < num_steps; ++j) {
+                solve_model_ode_cpu(dt, sv + (sv_id * NEQ), stim_currents[i], mapping[i]);
+            }
+        }
+    }
+    // Default: All cells are ENDO type
+    else {
+        OMP(parallel for private(sv_id))
+        for (i = 0; i < num_cells_to_solve; i++) {
+            if(cells_to_solve)
+                sv_id = cells_to_solve[i];
+            else
+                sv_id = i;
+
+            for (int j = 0; j < num_steps; ++j) {
+                solve_model_ode_cpu(dt, sv + (sv_id * NEQ), stim_currents[i], 0);
+            }
+        }
+    }
+}
+
+
+void solve_model_ode_cpu(real dt, real *sv, real stim_current, int type_cell)  {
+
+    assert(sv);
+
+    real rY[NEQ], rDY[NEQ];
+
+    for(int i = 0; i < NEQ; i++)
+        rY[i] = sv[i];
+
+    RHS_cpu(rY, rDY, stim_current, dt, type_cell);
+
+    // Transmembrane potencial is solved using Explicit Euler
+    sv[0] = dt*rDY[0] + sv[0];
+
+    // The other gate variables are solved using Rush-Larsen
+    for(int i = 1; i < NEQ; i++) {
+        sv[i] = rDY[i];
+    }
+}
+
+void RHS_cpu(const real *sv, real *rDY_, real stim_current, real dt, int type_cell) {
+    const real u   = sv[0];
+    const real v   = sv[1];
+    const real w   = sv[2];
+    const real s   = sv[3];
+
+    const real u_o = 0.0;
+    const real theta_v = 0.3;
+    const real theta_w = 0.13;
+    const real tau_vplus = 1.4506;
+    const real tau_s1 = 2.7342;
+    const real k_s = 2.0994;
+    const real u_s = 0.9087;
+
+    real u_u, theta_vminus, theta_o, tau_v1minus, tau_v2minus, tau_w1minus, tau_w2minus;
+    real k_wminus, u_wminus, tau_wplus, tau_fi, tau_o1, tau_o2, tau_so1, tau_so2;
+    real k_so, u_so, tau_s2, tau_si, tau_winf, w_infstar;
+
+    if (type_cell == 0) {        // ENDO
+        u_u = 1.56; theta_vminus = 0.2; theta_o = 0.006; tau_v1minus = 75.0; tau_v2minus = 10.0;
+        tau_w1minus = 6.0; tau_w2minus = 140.0; k_wminus = 200.0; u_wminus = 0.016;
+        tau_wplus = 280.0; tau_fi = 1.5*0.1; tau_o1 = 470.0; tau_o2 = 6.0; tau_so1 = 40.0;
+        tau_so2 = 1.2; k_so = 2.0; u_so = 0.65; tau_s2 = 2.0; tau_si = 2.9013;
+        tau_winf = 0.0273; w_infstar = 0.78;
+    } else if (type_cell == 1) { // MYO
+        u_u = 1.61; theta_vminus = 0.1; theta_o = 0.005; tau_v1minus = 80.0; tau_v2minus = 1.4506;
+        tau_w1minus = 70.0; tau_w2minus = 8.0; k_wminus = 200.0; u_wminus = 0.016;
+        tau_wplus = 280.0; tau_fi = 1.5*0.078; tau_o1 = 410.0; tau_o2 = 7.0; tau_so1 = 91.0;
+        tau_so2 = 0.8; k_so = 2.1; u_so = 0.6; tau_s2 = 4.0; tau_si = 3.3849;
+        tau_winf = 0.01; w_infstar = 0.5;
+    } else {                    // EPI
+        u_u = 1.55; theta_vminus = 0.006; theta_o = 0.006; tau_v1minus = 60.0; tau_v2minus = 1150.0;
+        tau_w1minus = 60.0; tau_w2minus = 15.0; k_wminus = 65.0; u_wminus = 0.03;
+        tau_wplus = 200.0; tau_fi = 1.5*0.11; tau_o1 = 400.0; tau_o2 = 6.0; tau_so1 = 30.0181;
+        tau_so2 = 0.9957; k_so = 2.0458; u_so = 0.65; tau_s2 = 16.0; tau_si = 1.8875;
+        tau_winf = 0.07; w_infstar = 0.94;
+    }
+
+    real H = (u - theta_v > 0) ? 1.0 : 0.0;
+    real h_o = (u - theta_o > 0) ? 1.0 : 0.0;
+    real h_w = (u - theta_w > 0) ? 1.0 : 0.0;
+    real h_v_minus = (u - theta_vminus > 0) ? 1.0 : 0.0;
+
+    real tau_o = (1.0 - h_o) * tau_o1 + h_o * tau_o2;
+    real tau_so = tau_so1 + (tau_so2 - tau_so1) * (1.0 + tanh(k_so * (u - u_so))) * 0.5;
+    real tau_vminus = (1.0 - h_v_minus) * tau_v1minus + h_v_minus * tau_v2minus;
+
+    real J_fi = -v * H * (u - theta_v) * (u_u - u) / tau_fi;
+    real J_so = (u - u_o) * (1.0 - h_w) / tau_o + h_w / tau_so;
+    real J_si = -h_w * w * s / tau_si;
+
+    rDY_[0] = -(J_fi + J_so + J_si) + stim_current;
+
+    real v_inf = (u < theta_vminus) ? 1.0 : 0.0;
+    real tau_v_rl = (tau_vplus * tau_vminus) / (tau_vplus - tau_vplus * H + tau_vminus * H);
+    real v_inf_rl = (tau_vplus * v_inf * (1 - H)) / (tau_vplus - tau_vplus * H + tau_vminus * H);
+
+    if (tau_v_rl > 1e-10) {
+        rDY_[1] = v_inf_rl - (v_inf_rl - v) * exp(-dt / tau_v_rl);
+    } else {
+        rDY_[1] = dt * ((1.0 - H) * (v_inf - v) / tau_vminus - H * v / tau_vplus) + v;
+    }
+
+    real w_inf = (1.0 - h_o) * (1.0 - (u / tau_winf)) + h_o * w_infstar;
+    real tau_wminus = tau_w1minus + (tau_w2minus - tau_w1minus) * (1.0 + tanh(k_wminus * (u - u_wminus))) * 0.5;
+    real tau_w_rl = (tau_wplus * tau_wminus) / (tau_wplus - tau_wplus * h_w + tau_wminus * h_w);
+    real w_inf_rl = (tau_wplus * w_inf * (1 - h_w)) / (tau_wplus - tau_wplus * h_w + tau_wminus * h_w);
+
+    if (tau_w_rl > 1e-10) {
+        rDY_[2] = w_inf_rl - (w_inf_rl - w) * exp(-dt / tau_w_rl);
+    } else {
+        rDY_[2] = dt * ((1.0 - h_w) * (w_inf - w) / tau_wminus - h_w * w / tau_wplus) + w;
+    }
+
+    real tau_s = (1.0 - h_w) * tau_s1 + h_w * tau_s2;
+    real s_inf_rl = (1.0 + tanh(k_s * (u - u_s))) / 2;
+
+    if (tau_s > 1e-10) {
+        rDY_[3] = s_inf_rl - (s_inf_rl - s) * exp(-dt / tau_s);
+    } else {
+        rDY_[3] = dt * (((1.0 + tanh(k_s * (u - u_s))) * 0.5 - s) / tau_s) + s;
+    }
+}
diff --git a/src/models_library/minimal_model/minimal_model.cu b/src/models_library/minimal_model/minimal_model.cu
new file mode 100644
index 00000000..ae38e99c
--- /dev/null
+++ b/src/models_library/minimal_model/minimal_model.cu
@@ -0,0 +1,209 @@
+#include "../../gpu_utils/gpu_utils.h"
+#include "../../monodomain/constants.h"
+#include <stddef.h>
+#include <assert.h>
+#include <stdlib.h>
+#include <math.h>
+#include "minimal_model.h"
+#include <stdio.h>
+
+extern "C" SET_ODE_INITIAL_CONDITIONS_GPU(set_model_initial_conditions_gpu) {
+
+    log_info("Using Minimal Model Mixed GPU\n");
+
+    uint32_t num_volumes = solver->original_num_cells;
+
+    // execution configuration
+    const int GRID  = (num_volumes + BLOCK_SIZE - 1)/BLOCK_SIZE;
+
+    size_t size = num_volumes*sizeof(real);
+
+    check_cuda_error(cudaMallocPitch((void **) &(solver->sv), &pitch_h, size, (size_t )NEQ));
+    check_cuda_error(cudaMemcpyToSymbol(pitch, &pitch_h, sizeof(size_t)));
+
+    kernel_set_model_inital_conditions <<<GRID, BLOCK_SIZE>>>(solver->sv, NULL, num_volumes);
+
+    check_cuda_error( cudaPeekAtLastError() );
+    cudaDeviceSynchronize();
+    return pitch_h;
+}
+
+extern "C" SOLVE_MODEL_ODES(solve_model_odes_gpu) {
+
+    size_t num_cells_to_solve = ode_solver->num_cells_to_solve;
+    uint32_t * cells_to_solve = ode_solver->cells_to_solve;
+    real *sv = ode_solver->sv;
+    real dt = ode_solver->min_dt;
+    uint32_t num_steps = ode_solver->num_steps;
+
+    // execution configuration
+    const int GRID  = ((int)num_cells_to_solve + BLOCK_SIZE - 1)/BLOCK_SIZE;
+
+    size_t stim_currents_size = sizeof(real)*num_cells_to_solve;
+    size_t cells_to_solve_size = sizeof(uint32_t)*num_cells_to_solve;
+
+    real *stims_currents_device;
+    check_cuda_error(cudaMalloc((void **) &stims_currents_device, stim_currents_size));
+    check_cuda_error(cudaMemcpy(stims_currents_device, stim_currents, stim_currents_size, cudaMemcpyHostToDevice));
+
+    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 *mapping = NULL;
+    uint32_t *mapping_device = NULL;
+    // Extra data function will tag the cells in the grid
+    if(ode_solver->ode_extra_data)
+    {
+        mapping = (uint32_t*)ode_solver->ode_extra_data;
+        check_cuda_error(cudaMalloc((void **)&mapping_device, ode_solver->extra_data_size));
+        check_cuda_error(cudaMemcpy(mapping_device, mapping, ode_solver->extra_data_size, cudaMemcpyHostToDevice));
+    }
+    // Default: All cells in the grid are ENDO type (tag=0.0)
+    else
+    {
+        check_cuda_error(cudaMalloc((void **)&mapping_device, ode_solver->extra_data_size));
+        check_cuda_error(cudaMemset(mapping_device, 0, cells_to_solve_size));
+    }
+
+    solve_gpu<<<GRID, BLOCK_SIZE>>>(dt, sv, stims_currents_device, cells_to_solve_device, num_cells_to_solve, num_steps, mapping_device);
+
+    check_cuda_error( cudaPeekAtLastError() );
+    check_cuda_error(cudaFree(stims_currents_device));
+
+    if(cells_to_solve_device) check_cuda_error(cudaFree(cells_to_solve_device));
+    if(mapping_device) check_cuda_error(cudaFree(mapping_device));
+}
+
+__global__ void kernel_set_model_inital_conditions(real *sv, real*IC, int num_volumes)
+{
+    // Thread ID
+    int threadID = blockDim.x * blockIdx.x + threadIdx.x;
+
+    if(threadID < num_volumes) {
+
+        *((real *) ((char *) sv + pitch * 0) + threadID) = 0.0f;        // u
+        *((real *) ((char *) sv + pitch * 1) + threadID) = 1.0f;        // v
+        *((real *) ((char *) sv + pitch * 2) + threadID) = 1.0f;        // w
+        *((real *) ((char *) sv + pitch * 3) + threadID) = 0.0f;        // s
+    }
+}
+
+// Solving the model for each cell in the tissue matrix ni x nj
+__global__ void solve_gpu(real dt, real *sv, real* stim_currents,
+                          uint32_t *cells_to_solve, uint32_t num_cells_to_solve,
+                          int num_steps, uint32_t *mapping)
+{
+    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;
+
+        real rDY[NEQ];
+
+        for (int n = 0; n < num_steps; ++n) {
+
+            RHS_gpu(sv, rDY, stim_currents[threadID], sv_id, dt, mapping[sv_id]);
+
+            // Transmembrane potencial is solved using Explicit Euler
+            *((real*)((char*)sv) + sv_id) = dt*rDY[0] + *((real*)((char*)sv) + sv_id);
+
+            // The other gate variables are solved using Rush-Larsen
+            for(int i = 1; i < NEQ; i++) {
+                *((real*)((char*)sv + pitch * i) + sv_id) = rDY[i];
+            }
+        }
+    }
+}
+
+inline __device__ void RHS_gpu(real *sv_, real *rDY_, real stim_current, int threadID_, real dt, int type_cell) {
+
+    const real u   = *((real*)((char*)sv_ + pitch * 0) + threadID_);
+    const real v   = *((real*)((char*)sv_ + pitch * 1) + threadID_);
+    const real w   = *((real*)((char*)sv_ + pitch * 2) + threadID_);
+    const real s   = *((real*)((char*)sv_ + pitch * 3) + threadID_);
+
+    const real u_o = 0.0;
+    const real theta_v = 0.3;
+    const real theta_w = 0.13;
+    const real tau_vplus = 1.4506;
+    const real tau_s1 = 2.7342;
+    const real k_s = 2.0994;
+    const real u_s = 0.9087;
+
+    real u_u, theta_vminus, theta_o, tau_v1minus, tau_v2minus, tau_w1minus, tau_w2minus;
+    real k_wminus, u_wminus, tau_wplus, tau_fi, tau_o1, tau_o2, tau_so1, tau_so2;
+    real k_so, u_so, tau_s2, tau_si, tau_winf, w_infstar;
+
+    if (type_cell == 0) {        // ENDO
+        u_u = 1.56; theta_vminus = 0.2; theta_o = 0.006; tau_v1minus = 75.0; tau_v2minus = 10.0;
+        tau_w1minus = 6.0; tau_w2minus = 140.0; k_wminus = 200.0; u_wminus = 0.016;
+        tau_wplus = 280.0; tau_fi = 1.5*0.1; tau_o1 = 470.0; tau_o2 = 6.0; tau_so1 = 40.0;
+        tau_so2 = 1.2; k_so = 2.0; u_so = 0.65; tau_s2 = 2.0; tau_si = 2.9013;
+        tau_winf = 0.0273; w_infstar = 0.78;
+    } else if (type_cell == 1) { // MYO
+        u_u = 1.61; theta_vminus = 0.1; theta_o = 0.005; tau_v1minus = 80.0; tau_v2minus = 1.4506;
+        tau_w1minus = 70.0; tau_w2minus = 8.0; k_wminus = 200.0; u_wminus = 0.016;
+        tau_wplus = 280.0; tau_fi = 1.5*0.078; tau_o1 = 410.0; tau_o2 = 7.0; tau_so1 = 91.0;
+        tau_so2 = 0.8; k_so = 2.1; u_so = 0.6; tau_s2 = 4.0; tau_si = 3.3849;
+        tau_winf = 0.01; w_infstar = 0.5;
+    } else {                    // EPI
+        u_u = 1.55; theta_vminus = 0.006; theta_o = 0.006; tau_v1minus = 60.0; tau_v2minus = 1150.0;
+        tau_w1minus = 60.0; tau_w2minus = 15.0; k_wminus = 65.0; u_wminus = 0.03;
+        tau_wplus = 200.0; tau_fi = 1.5*0.11; tau_o1 = 400.0; tau_o2 = 6.0; tau_so1 = 30.0181;
+        tau_so2 = 0.9957; k_so = 2.0458; u_so = 0.65; tau_s2 = 16.0; tau_si = 1.8875;
+        tau_winf = 0.07; w_infstar = 0.94;
+    }
+
+    real H = (u - theta_v > 0) ? 1.0 : 0.0;
+    real h_o = (u - theta_o > 0) ? 1.0 : 0.0;
+    real h_w = (u - theta_w > 0) ? 1.0 : 0.0;
+    real h_v_minus = (u - theta_vminus > 0) ? 1.0 : 0.0;
+
+    real tau_o = (1.0 - h_o) * tau_o1 + h_o * tau_o2;
+    real tau_so = tau_so1 + (tau_so2 - tau_so1) * (1.0 + tanh(k_so * (u - u_so))) * 0.5;
+    real tau_vminus = (1.0 - h_v_minus) * tau_v1minus + h_v_minus * tau_v2minus;
+
+    real J_fi = -v * H * (u - theta_v) * (u_u - u) / tau_fi;
+    real J_so = (u - u_o) * (1.0 - h_w) / tau_o + h_w / tau_so;
+    real J_si = -h_w * w * s / tau_si;
+
+    rDY_[0] = -(J_fi + J_so + J_si) + stim_current;
+
+    real v_inf = (u < theta_vminus) ? 1.0 : 0.0;
+    real tau_v_rl = (tau_vplus * tau_vminus) / (tau_vplus - tau_vplus * H + tau_vminus * H);
+    real v_inf_rl = (tau_vplus * v_inf * (1 - H)) / (tau_vplus - tau_vplus * H + tau_vminus * H);
+
+    if (tau_v_rl > 1e-10) {
+        rDY_[1] = v_inf_rl - (v_inf_rl - v) * exp(-dt / tau_v_rl);
+    } else {
+        rDY_[1] = dt * ((1.0 - H) * (v_inf - v) / tau_vminus - H * v / tau_vplus) + v;
+    }
+
+    real w_inf = (1.0 - h_o) * (1.0 - (u / tau_winf)) + h_o * w_infstar;
+    real tau_wminus = tau_w1minus + (tau_w2minus - tau_w1minus) * (1.0 + tanh(k_wminus * (u - u_wminus))) * 0.5;
+    real tau_w_rl = (tau_wplus * tau_wminus) / (tau_wplus - tau_wplus * h_w + tau_wminus * h_w);
+    real w_inf_rl = (tau_wplus * w_inf * (1 - h_w)) / (tau_wplus - tau_wplus * h_w + tau_wminus * h_w);
+
+    if (tau_w_rl > 1e-10) {
+        rDY_[2] = w_inf_rl - (w_inf_rl - w) * exp(-dt / tau_w_rl);
+    } else {
+        rDY_[2] = dt * ((1.0 - h_w) * (w_inf - w) / tau_wminus - h_w * w / tau_wplus) + w;
+    }
+
+    real tau_s = (1.0 - h_w) * tau_s1 + h_w * tau_s2;
+    real s_inf_rl = (1.0 + tanh(k_s * (u - u_s))) / 2;
+
+    if (tau_s > 1e-10) {
+        rDY_[3] = s_inf_rl - (s_inf_rl - s) * exp(-dt / tau_s);
+    } else {
+        rDY_[3] = dt * (((1.0 + tanh(k_s * (u - u_s))) * 0.5 - s) / tau_s) + s;
+    }
+}
\ No newline at end of file
diff --git a/src/models_library/minimal_model/minimal_model.h b/src/models_library/minimal_model/minimal_model.h
new file mode 100644
index 00000000..4a6d05e9
--- /dev/null
+++ b/src/models_library/minimal_model/minimal_model.h
@@ -0,0 +1,36 @@
+#ifndef MONOALG3D_MODEL_MINIMAL_MODEL_H
+#define MONOALG3D_MODEL_MINIMAL_MODEL_H
+
+#include "../model_common.h"
+#include "../../extra_data_library/helper_functions.h"
+
+/*
+Minimal model for human ventricular action potentials in tissue,
+Alfonso Bueno-Orovio, Elizabeth M Cherry, Flavio H Fenton,
+(2008)
+*/
+
+#define NEQ 4
+#define INITIAL_V (0.0f)
+
+#ifdef __CUDACC__
+
+#include "../../gpu_utils/gpu_utils.h"
+
+static __device__ size_t pitch;
+static size_t pitch_h;
+
+__global__ void kernel_set_model_inital_conditions(real *sv, real* ICs, int num_volumes);
+
+__global__ void solve_gpu(real dt, real *sv, real *stim_currents, uint32_t *cells_to_solve, uint32_t num_cells_to_solve,
+                          int num_steps,uint32_t *mapping);
+
+inline __device__ void RHS_gpu(real *sv_, real *rDY_, real stim_current, int threadID_, real dt,
+                                int type_cell);
+
+#endif
+
+void RHS_cpu(const real *sv, real *rDY_, real stim_current, real dt, int type_cell);
+void solve_model_ode_cpu(real dt, real *sv, real stim_current, int type_cell);
+
+#endif
\ No newline at end of file
diff --git a/src/monodomain/monodomain_solver.c b/src/monodomain/monodomain_solver.c
index d3be7adc..f1d885fe 100644
--- a/src/monodomain/monodomain_solver.c
+++ b/src/monodomain/monodomain_solver.c
@@ -653,7 +653,7 @@ int solve_monodomain(struct monodomain_solver *the_monodomain_solver, struct ode
 
         if(calc_ecg && (count % calc_ecg_rate == 0)) {
             start_stop_watch(&stop_watch);
-            ((calc_ecg_fn *)calc_ecg_config->main_function)(&time_info, calc_ecg_config, the_grid);
+            ((calc_ecg_fn *)calc_ecg_config->main_function)(&time_info, calc_ecg_config, the_grid, out_dir_name);
             total_ecg_time += stop_stop_watch(&stop_watch);
         }
 
@@ -1599,7 +1599,7 @@ void write_pmj_delay(struct grid *the_grid, struct config *config, struct termin
                         // fprintf(output_file,"ERROR! Probably there was a block on the anterograde direction!\n");
                         has_block = true;
                         cur_pulse = 0;
-                        // return;
+                        return;
                     }
                     mean_tissue_lat += activation_times_array_tissue[cur_pulse];
                     if(activation_times_array_tissue[cur_pulse] < min_tissue_lat) {