Fixing demo for length_tension

code/FiberCpp/half_sarcomere.cpp

@@ -2602,6 +2602,44 @@ void half_sarcomere::thin_filament_kinetics(double time_step, double Ca_conc)
             }
         }
 
+        // Try the k_on
+        if (!gsl_isnan(gsl_vector_get(p_fs_model->a_k_on_t_force, 0)))
+        {
+            a_k_on_t_force_factor = 1.0 +
+                gsl_vector_get(p_fs_model->a_k_on_t_force, 0) *
+                ((gsl_vector_get(p_fs_model->a_k_on_t_force, 1) * hs_titin_force) +
+                    (gsl_vector_get(p_fs_model->a_k_on_t_force, 2) * t_force_across_M) +
+                    (gsl_vector_get(p_fs_model->a_k_on_t_force, 3) * t_force_across_Z));
+
+            // Limit
+            if (!gsl_isnan(gsl_vector_get(p_fs_model->a_k_on_t_force, 4)))
+                a_k_on_t_force_factor = GSL_MAX(a_k_on_t_force_factor,
+                    gsl_vector_get(p_fs_model->a_k_on_t_force, 4));
+
+            if (!gsl_isnan(gsl_vector_get(p_fs_model->a_k_on_t_force, 5)))
+                a_k_on_t_force_factor = GSL_MIN(a_k_on_t_force_factor,
+                    gsl_vector_get(p_fs_model->a_k_on_t_force, 5));
+        }
+
+        // Try the k_off
+        if (!gsl_isnan(gsl_vector_get(p_fs_model->a_k_off_t_force, 0)))
+        {
+            a_k_off_t_force_factor = 1.0 +
+                gsl_vector_get(p_fs_model->a_k_off_t_force, 0) *
+                ((gsl_vector_get(p_fs_model->a_k_off_t_force, 1) * hs_titin_force) +
+                    (gsl_vector_get(p_fs_model->a_k_off_t_force, 2) * t_force_across_M) +
+                    (gsl_vector_get(p_fs_model->a_k_off_t_force, 3) * t_force_across_Z));
+
+            // Limit
+            if (!gsl_isnan(gsl_vector_get(p_fs_model->a_k_off_t_force, 4)))
+                a_k_off_t_force_factor = GSL_MAX(a_k_off_t_force_factor,
+                    gsl_vector_get(p_fs_model->a_k_off_t_force, 4));
+
+            if (!gsl_isnan(gsl_vector_get(p_fs_model->a_k_off_t_force, 5)))
+                a_k_off_t_force_factor = GSL_MIN(a_k_off_t_force_factor,
+                    gsl_vector_get(p_fs_model->a_k_off_t_force, 5));
+        }
+
         // Try the coop
         if (!gsl_isnan(gsl_vector_get(p_fs_model->a_k_coop_t_force, 0)))
         {
@@ -2683,7 +2721,7 @@ void half_sarcomere::thin_filament_kinetics(double time_step, double Ca_conc)
                         coop_boost = a_gamma_coop *
                             (double)gsl_vector_short_get(p_af[a_counter]->active_neighbors, unit_counter);
 
-                        rate = a_k_on * Ca_conc * (1.0 + (a_k_coop_t_force_factor * coop_boost));
+                        rate = (a_k_on_t_force_factor * a_k_on) * Ca_conc * (1.0 + (a_k_coop_t_force_factor * coop_boost));
 
                         // Test event with a random number
                         rand_double = gsl_rng_uniform(rand_generator);
@@ -2713,7 +2751,7 @@ void half_sarcomere::thin_filament_kinetics(double time_step, double Ca_conc)
                             coop_boost = a_gamma_coop *
                                 (2.0 - (double)gsl_vector_short_get(p_af[a_counter]->active_neighbors, unit_counter));
 
-                            rate = a_k_off * (1.0 + (a_k_coop_t_force_factor * coop_boost));
+                            rate = (a_k_off_t_force_factor * a_k_off) * (1.0 + (a_k_coop_t_force_factor * coop_boost));
 
                             // Test event with a random number
                             rand_double = gsl_rng_uniform(rand_generator);

code/FiberPy/FiberPy/package/modules/display/myofibril.py

@@ -8,6 +8,7 @@
 import os
 import json
 import cv2
+import math
 
 import numpy as np
 import pandas as pd
@@ -24,13 +25,14 @@ def draw_myofibril_frame(df, frame=1, image_file_string=[]):
     """ Draws a myofibril """
 
     # Variables
+
+    hsl_region = [0.45, 0.55]
+
     fig_hs_width = 0.5
     fig_hs_x_padding = 0.5
     fig_hs_height = 0.5
     fig_hs_y_padding = 0.5
 
-    ax_myofibril = 0
-
     z_width = 100
     z_line_color = 'black'
 
@@ -45,21 +47,31 @@ def draw_myofibril_frame(df, frame=1, image_file_string=[]):
     thick_width = 0.3
     thick_color = 'turquoise'
 
+    no_of_rows = 2
+    no_of_cols = 1
+
+    ax_myofibril = 0
+    ax_sl = 1
 
     # Work out how many half-sarcomeres there are
     command_length_cols = [x for x in d.columns.to_list() if x.endswith('command_length')]
     no_of_half_sarcomeres = len(command_length_cols)
 
     # Work out fig dimensions
     fig = plt.figure(constrained_layout=False)
-    gs = fig.add_gridspec(nrows=1, ncols=1)
+    gs = fig.add_gridspec(nrows=no_of_rows, ncols=no_of_cols)
 
     fig_width = (2 * fig_hs_x_padding) + (no_of_half_sarcomeres * fig_hs_width)
     fig_height = (2 * fig_hs_y_padding) + fig_hs_height
 
     fig.set_size_inches([fig_width, fig_height])
+
     ax = []
-    ax.append(fig.add_subplot(gs[0,0]))
+    for i in range(no_of_rows):
+        for j in range(no_of_cols):
+            ax.append(fig.add_subplot(gs[i,j]))
+
+
 
     # Plot half_sarcomere
     x_anchor = 0
@@ -102,9 +114,24 @@ def draw_myofibril_frame(df, frame=1, image_file_string=[]):
         # Move the anchor
         if (hs_polarity == 1):        
             x_anchor = x_anchor + hs_length
+
+    # Get half-sarcomere length in center
+    hs_left = math.floor(hsl_region[0] * no_of_half_sarcomeres)
+    hs_right = math.ceil(hsl_region[1] * no_of_half_sarcomeres)
+
+    # Make a new column
+    d['hsl_center'] = 0
+    for i in range(hs_left, hs_right):
+        hsl_string = 'hs_%i_length' % i
+        d['hsl_center'] = d['hsl_center'] + d[hsl_string]
+
+    d['hsl_center'] = (1 / float(hs_right - hs_left)) * d['hsl_center']
+
+
+    ax[1].plot(d['time'], d['hsl_center'], 'b-')
+
 
-    plt.xlim(0, 70000)
-    plt.ylim(-1, 2)
+    ax[ax_myofibril].set_xlim(0, 70000)
 
 
 
@@ -137,7 +164,7 @@ def create_movie_from_image_folder(image_folder, movie_file_string):
 
 
 if __name__ == "__main__":
-    fs = 'd:/ken/github/campbellmusclelab/publications/paper_fibersim_sequential_relaxation/simulations/kens_playground/f/sim_data/sim_output/1/sim_pCa_59_s_0_r1.txt'
+    fs = 'd:/ken/github/campbellmusclelab/publications/paper_fibersim_sequential_relaxation/simulations/kens_playground/g/sim_data/sim_output/3/sim_pCa_59_s_0_r1.txt'
 
     d = pd.read_csv(fs, sep='\t')
 

docs/pages/demos/length_tension/steady_state/steady_state.md

@@ -30,12 +30,12 @@ If you need help with these step, check the [installation instructions](../../..
 + Change directory to `<FiberSim_repo>/code/FiberPy/FiberPy`
 + Run the command
 ```
- python FiberPy.py characterize "../../../demo_files/isotonic/force_velocity/base/setup.json"
+ python FiberPy.py characterize "../../../demo_files/length_tension/steady_state/base/setup.json"
  ```
 
 ### Viewing the results
 
-All of the results from the isotonic simulations are written to files in `<FiberSim_repo>/demo_files/isotonic/force_velocity/sim_data/isotonic/sim_output`
+All of the results from the isotonic simulations are written to files in `<FiberSim_repo>/demo_files/length_tension/steady_state/sim_output`
 
 The file `superposed_traces.png` shows pCa, length, force per cross-sectional area (stress), and thick and thin filamnt properties plotted against time.
 
@@ -45,28 +45,25 @@ The file `rates.png` summarizes the kinetic scheme.
 
 <img src="images/rates.png" width="50%">
 
-The file `fv_traces_1.png` shows the loaded shortening segment of the simulation in more details.
+The file `length_tension.png` shows the active, passive, and total length-tension curves.
 
-<img src="images/fv_traces_1.png" width="50%">
+<img src="images/length_tension.png" width="50%">
 
-The file `fv_and_power.png` shows the force-velocity and force-power curves.
 
-<img src="images/fv_and_power.png" width="50%">
 
-Finally, data derived from the simulations are tabulated in `fv_analysis.xlsx`.
-
-<img src="images/fv_analysis.png" width="50%">
+### How this worked
 
+The setup file follows the normal template. The experimental protocols are defined by the `characterization` element.
 
+The new feature in this demo is the `hs_lengths` option in the characterization structure. As you can see from the output figures, the code simulates a pCa 4.5 activation at half-sarcomere lengths ranging from 700 to 1800 nm.
 
+The `post_sim_Python_call` generates the length-tension figure which is unique to this demo.
 
-### How this worked
-
-The setup file follows the normal template. The experimental protocols are defined by the `characterization` element. See the table below for more details.
 
 ```text
 {
-  "FiberSim_setup": {
+  "FiberSim_setup":
+  {
     "FiberCpp_exe": {
       "relative_to": "this_file",
       "exe_file": "../../../../bin/FiberCpp.exe"
@@ -77,40 +74,22 @@ The setup file follows the normal template. The experimental protocols are defin
       "model_files": ["model.json"]
     },
     "characterization": [
-      {
-        "type": "force_velocity",
-        "pCa": 4.5,
-        "hs_lengths": [1100],
-        "m_n": 4,
-        "randomized_repeats": 1,
-        "length_fit_mode": "exponential",
-        "time_step_s": 0.001,
-        "sim_duration_s": 0.5,
-        "sim_release_s": 0.4,
-        "rel_isotonic_forces": [0.03, 0.04, 0.075, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.85, 0.875, 0.9, 0.925, 0.95],
-        "fit_time_s": [ 0.41, 0.49 ],
-        "relative_to": "this_file",
-        "sim_folder": "../sim_data",
-        "output_image_formats": [ "png" ],
-        "figures_only": "False",
-        "trace_figures_on": "False"
-      }
+        {
+            "type": "pCa_length_control",
+            "relative_to": "this_file",
+            "sim_folder": "../sim_data",
+            "hs_lengths": [700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800],
+            "m_n": 4,
+            "pCa_values": [4.5],
+            "sim_duration_s": 1,
+            "time_step_s": 0.001,
+            "pCa_step_up_s": 0.1,
+            "output_image_formats": [ "png" ],
+            "figures_only": "False",
+            "trace_figures_on": "False",
+            "post_sim_Python_call": "../Python_code/plot_length_tension.py"
+        }
     ]
   }
 }
 ```
-
-| Parameter | Value | Comments |
-| ---- | ---- | ---- |
-| type | force_velocity | |
-| pCa | float | The pCa value for the tests (as in this example) |
-|  | "pCa_50" | Alternative mode, which runs isotonic shortening at pCa<sub>50</sub>. If this mode is selected, an additional parameter `pCa_values` must also be provided. FiberPy will run isometric simulations with these values to determine the pCa<sub>50</sub>. Example, `"pCa_values": [9, 6.5, 6.3, 6.1, 5.9, 5.7, 5.5, 5.3, 4.5]` |
-| hs_lengths (optional)| array of floats | array of hs_lengths to run the simulations at. FiberPy uses the hs_length in the model file if not provided |
-| randomized_repeats (optional) | integer | number of repeats to run at each length with different seeds |
-| length_fit_mode | "linear" | fits shortening velocity with a straight line |
-|  | "exponential" | fits shortening velocity with an exponential and extrapolates to the onset of shortening to establish fastest shortening (useful for traces where shortening slows progressively during shortening) |
-| time_step_s | float | time step in seconds for simulation |
-| sim_duration_s | float | total duration in seconds of each simulation |
-| sim_release_s | float | time in seconds at which muscle starts to shorten isotonically |
-| rel_isotonic_forces | array of floats | the forces (relative to isometric) which the muscle shortens against |
-| fit_time_s | [2 element array of floats] | time range (in seconds) over which to fit shortening velocity |