Initial commit of JDP titin addition

Joe Powers spent a considerable amount of effort incorporating titin into the model. Here’s the lion’s share of his work added in, starting from the point he initially pulled from the repo.

Author: cdw

multifil/af.py

@@ -190,6 +190,10 @@ def from_dict(self, tfd):
         self.binding_sites = [self.parent_thin.resolve_address(bsa) \
                               for bsa in tfd['binding_sites']]
 
+    def link_titin(self, titin_fil):
+        """Add a titin filament to this face"""
+        self.titin_fil = titin_fil
+
     def nearest(self, axial_location):
         """Where is the nearest binding site?
 

multifil/hs.py

@@ -16,6 +16,7 @@
 import scipy.optimize as opt
 from . import af
 from . import mf
+from . import ti
 
 class hs:
     """The half-sarcomere and ways to manage it"""
@@ -249,6 +250,56 @@ def __init__(self, lattice_spacing=None, z_line=None, poisson=None,
         self.actin_permissiveness = actin_permissiveness
         # Track how long we've been running
         self.current_timestep = 0
+        # Create the titin filaments and link them from thick
+        # faces to thin faces
+        # |--------------------------------------------------|
+        # |            Actin & titin around myosin           |
+        # |--------------------------------------------------|
+        # |           a1               a3                    |
+        # |                                                  |
+        # |  a0       t1      a2       t4       a0           |
+        # |       t0     t2        t3      t5                |
+        # |           M0               M1                    |
+        # |       t6     t8        t9      t11               |
+        # |  a4       t7      a6       t10      a4           |
+        # |                                                  |
+        # |           a5     t13       a7       t16      a5  |
+        # |               t12    t14        t15    t17       |
+        # |                   M2                M3           |
+        # |               t18    t20        t21    t23   a1  |
+        # |           a1      t19      a3       t22          |
+        # |                                                  |
+        # |                   a2                a0           |
+        # |--------------------------------------------------|
+        ## CHECK_JDP ## Link Thick filament to titin
+        ti_thick = lambda i, j: self.thick[i].thick_faces[j]
+        ti_thin = lambda i, j: self.thin[i].thin_faces[j]
+        self.titin = (
+            ti.Titin(self, 0, ti_thick(0, 0), ti_thin(0, 1)),
+            ti.Titin(self, 1, ti_thick(0, 1), ti_thin(1, 2)),
+            ti.Titin(self, 2, ti_thick(0, 2), ti_thin(2, 2)),
+            ti.Titin(self, 3, ti_thick(1, 0), ti_thin(2, 1)),
+            ti.Titin(self, 4, ti_thick(1, 1), ti_thin(3, 2)),
+            ti.Titin(self, 5, ti_thick(1, 2), ti_thin(0, 2)),
+            ti.Titin(self, 6, ti_thick(0, 5), ti_thin(4, 1)),
+            ti.Titin(self, 7, ti_thick(0, 4), ti_thin(5, 0)),
+            ti.Titin(self, 8, ti_thick(0, 3), ti_thin(6, 0)),
+            ti.Titin(self, 9, ti_thick(1, 5), ti_thin(6, 1)),
+            ti.Titin(self, 10, ti_thick(1, 4), ti_thin(7, 0)),
+            ti.Titin(self, 11, ti_thick(1, 3), ti_thin(4, 0)),
+            ti.Titin(self, 12, ti_thick(2, 0), ti_thin(5, 1)),
+            ti.Titin(self, 13, ti_thick(2, 1), ti_thin(6, 2)),
+            ti.Titin(self, 14, ti_thick(2, 2), ti_thin(7, 2)),
+            ti.Titin(self, 15, ti_thick(3, 0), ti_thin(7, 1)),
+            ti.Titin(self, 16, ti_thick(3, 1), ti_thin(4, 2)),
+            ti.Titin(self, 17, ti_thick(3, 2), ti_thin(5, 2)),
+            ti.Titin(self, 18, ti_thick(2, 5), ti_thin(1, 1)),
+            ti.Titin(self, 19, ti_thick(2, 4), ti_thin(2, 0)),
+            ti.Titin(self, 20, ti_thick(2, 3), ti_thin(3, 0)),
+            ti.Titin(self, 21, ti_thick(3, 5), ti_thin(3, 1)),
+            ti.Titin(self, 22, ti_thick(3, 4), ti_thin(0, 0)),
+            ti.Titin(self, 23, ti_thick(3, 3), ti_thin(1, 0)),
+        )
 
     def to_dict(self):
         """Create a JSON compatible representation of the thick filament
@@ -489,7 +540,8 @@ def radialtension(self):
 
     def radialforce(self):
         """The sum of the thick filaments' radial forces, as a (y,z) vector"""
-        return np.sum([t.radial_force_of_filament() for t in self.thick], 0)
+        return np.sum([t.radial_force_of_filament() for t in self.thick], 0)# + ...
+        #sum([titin.radialforce() for titin in self.titin]) #CHECK
 
     def _single_settle(self):
         """Settle down now, just a little bit"""
@@ -516,6 +568,25 @@ def _get_residual(self):
         mash = np.hstack([thick_f, thin_f])
         return mash
 
+    def length_perturbation(self, dist=None, n=None):
+        """Get the force response of the half-sarcomere to a length perturbation"""
+        if dist is None:
+            dist = 0.2; #take length perturbation steps of 'dist' nm
+        if n is None:
+            n = 10; #take 'n' length perturbation steps.
+        response = [];
+        stiffness = [];
+        initial_force = self.axialforce()
+        initial_zline = self.z_line
+        for i in range(n):
+            self.z_line += dist
+            self.settle()
+            response.append(((i+1) * dist, self.axialforce() - initial_force))
+            stiffness.append((self.axialforce() - initial_force) / ((i+1) * dist))
+        self.z_line = initial_zline
+        self.settle()
+        return np.mean(stiffness)
+
     def get_frac_in_states(self):
         """Calculate the fraction of cross-bridges in each state"""
         nested = [t.get_states() for t in self.thick]
@@ -524,6 +595,12 @@ def get_frac_in_states(self):
         frac_in_state = [n/float(len(xb_states)) for n in num_in_state]
         return frac_in_state
 
+    #ADDED BY JDP ON 2017-Aug-21
+    def get_31_trans(self):
+        """Calculate the number of xb's that transition from state 3 to 1"""
+        xb_trans = sum(sum(self.last_transitions,[]),[])
+        return xb_trans.count('31')
+
     def update_ls_from_poisson_ratio(self):
         """Update the lattice spacing consistant with the poisson ratio,
         initial lattice spacing, current z-line, and initial z-line

multifil/mf.py

@@ -7,6 +7,7 @@
 """
 
 from . import mh
+from . import ti
 import numpy as np
 
 
@@ -223,6 +224,10 @@ def to_dict(self):
                                      for xb in thickfaced['xb_by_crown']]
         return thickfaced
 
+    def link_titin(self, titin_fil):
+        """Add a titin filament to this face"""
+        self.titin_fil = titin_fil
+
     def from_dict(self, tfd):
         """ Load values from a thick face dict. Values read in correspond
         to the current output documented in to_dict.
@@ -258,6 +263,7 @@ def axialforce(self):
     def radialtension(self):
         """Sum of the absolute values of radial force for every myosin"""
         radial = [crossbridge.radialforce() for crossbridge in self.xb]
+        radial.append(self.titin.radialforce())
         return sum(radial)
 
     def radialforce(self):
@@ -563,7 +569,6 @@ def settle(self):
         forces = self.axialforce()
         # Individual displacements needed to balance force
         isolated = 0.95*forces/self.k
-        isolated[-1] *= 2 # Last node has spring on only one side
         # Cumulative displacements
         cumulative = np.cumsum(isolated)
         # New axial locations
@@ -648,11 +653,28 @@ def _axial_thick_filament_forces(self, axial_locations=None):
         # Find the distance from crown to crown, then the resulting forces
         dists = np.diff(axial_locations)
         spring_force = (dists - self.rests) * self.k
-        spring_force = np.hstack([spring_force, 0]) # Last node not connected
-        # Zero would be the force of titin, were it in the model
+        # Location zero is the force of titin
+        net_force_at_crown = np.diff(spring_force)
+        titin_force = self._normed_total_titin_force()
+        spring_force = np.hstack([spring_force, titin_force])
         net_force_at_crown = np.diff(spring_force)
         return net_force_at_crown
 
+    def _normed_total_titin_force(self):
+        """Settle expects to move nodes to satisfy springs of stiffness self.k.
+        Titin has a different stiffness. We express the force titin is
+        generating in terms of the thick fil stiffness in order to treat the
+        movement necessary to balance the node attached to titin the same as
+        the movement necessary to balance the force of the thick filaments.
+        """
+        normed_titin_forces = []
+        for thick_face in self.thick_faces:
+            titin = thick_face.titin_fil
+            titin_force = titin.axialforce()
+            normed = titin_force * titin.stiffness() / self.k
+            normed_titin_forces.append(normed)
+        return np.sum(normed_titin_forces)
+
     @staticmethod
     def _radial_force_to_components(crown_force, orientation):
         """Convert radial components of a crown's force into a y,z vector

multifil/ti.py

@@ -0,0 +1,116 @@
+#!/usr/bin/env python
+# encoding: utf-8
+"""
+ti.py - A titin filament, with variable compliance
+
+MORE ABOUT HOW TITIN WORKS
+
+Created by Joe Powers and Dave Williams on 2017-02-17
+"""
+
+
+import numpy as np
+import warnings
+import numpy.random as random
+random.seed() # Ensure proper seeding
+
+
+class Titin:
+    """This is all about the titin filament"""
+    def __init__(self, parent_lattice, index, thick_face, thin_face,
+        a=240, b=0.0045):
+        """Initialize the titin filament.
+
+        Parameters:
+            parent_lattice: calling half-sarcomere instance
+            index: which titin filament this is (0-23)
+            thick_face: List of thick filament faces' numerical orientation (0-5)
+            thin_face: List of thin filament faces' numerical orientation (0-2)
+
+        Returns:
+            None
+        """
+        # Name of the titin molecule
+        self.index = index
+        # A connection to the parent lattice
+        self.parent_lattice = parent_lattice
+        # Location of the titin filament relative thick filament
+        self.thick_face = thick_face
+        # Link titin to the thick filament face
+        self.thick_face.link_titin(self)
+        # location of the titin filament relative to the thin filament
+        self.thin_face = thin_face
+        # Link titin to that face of the thin filament
+        self.thin_face.link_titin(self)
+        ## And now we declare titin properties that will be used to
+        ## calculate forces
+        self.rest = 120 #nm, no citation
+        # Create the constants that determine force NOTE IMPROVE DOC
+        self.a = a
+        self.b = b
+
+    def angle(self):
+        """Caclulate the angle the titin makes relative to thick filament"""
+        act_loc = self.thin_face.parent_thin.parent_lattice.z_line
+        myo_loc = self.thick_face.get_axial_location(-1)
+        ls = self.parent_lattice.lattice_spacing
+        angle = np.arctan2(ls, act_loc-myo_loc)
+        return angle
+
+    def length(self):
+        """Calculate the length of the titin filament"""
+        act_loc = self.thin_face.parent_thin.parent_lattice.z_line
+        myo_loc = self.thick_face.get_axial_location(-1)
+        ls = self.parent_lattice.lattice_spacing
+        length = np.sqrt( (act_loc-myo_loc)**2 + ls**2 )
+        return length
+
+    def stiffness(self):
+        """Need instintanious stiffness at the current length to normalize
+        force for settling at each timestep. We get this as the derivitive of
+        force with respect to x. D[a*exp(b*x), x] = a*b*exp(b*x)
+        """
+        return self.force()*self.b
+
+    def force(self):
+        """Calculate the total force the titin filament exerts"""
+        length = self.length()
+        if length < self.rest:
+            return 0 # titin buckles
+        else:
+            #Exponential Model
+            x = length - self.rest
+            return self.a*np.exp(self.b*x)
+        #else:
+            #Cubic model
+            #x = length - self.rest
+            #a = 3.25e-5
+            #b = -5e-3
+            #c = 9e-1
+            #d = -7e-5
+            #return a*(x**3) + b*(x**2) + c*x + d #Linke et al. PNAS 1998
+        #else:
+            # Sawtooth Model
+            #extension = length - self.rest
+            #a = 1.2e-6
+            #b = 3
+            #c = 0.5e-6#*random.rand()
+            #d = 34
+            #e = 3
+            #return a*extension**b + c*np.fmod(extension,d)**e
+        #else:
+            #Hookean Model
+            #return self.k * (self.length() - self.rest)
+
+    def axialforce(self):
+        """Return the total force the titin filament exerts on the
+        thick filament's final node, (negate for the thin filament side)
+        """
+        return self.force() * np.cos(self.angle()) ## CHECK_JDP ##
+
+    def radialforce(self):
+        """Return the force in the radial direction (positive is compressive)
+        TODO: The 'positive is compressive' part needs to be double checked
+        """
+        warnings.warn("Check radial force direction in titin")
+        return self.force() * np.sin(self.angle())