Cross linkers (class CrossLinkedNetwork, DoubleEndedCrossLinkedNetwork(extends CrossLinkedNetwork))

class CrossLinkedNetwork.CrossLinkedNetwork(nFib, N, Nunisites, Lfib, kCL, rl, Dom, fibDisc, smoothForce, nThreads=1)

This class keeps track of cross linking information.

In our model, cross linkers form and break at \(N_u\) uniformly-spaced “sites” along the fiber centerline. The positions of these sites are obtained by resampling the Chebyshev interpolant \(X(s)\) of the fiber centerline. This abstract class is not aware of the dynamics of link formation and breakage at the sites (those are left to the child class DoubleEndedCrossLinkedNetwork). This class is, however, aware of the number of links, and contains members which list the site number which is the “head” of each link and the site number which is the “tail.” It also stores a list of periodic shifts associated with each link, which can also be thought of as keeping track of what periodic copy of fiber \(i\) is linked to what periodic copy of fiber \(j\).

This class also implements methods for computing force and stress, and information about the network morphology (number of bundles, etc.) As far as computing the force, we model each CL as a linear spring between the two sites.

__init__(nFib, N, Nunisites, Lfib, kCL, rl, Dom, fibDisc, smoothForce, nThreads=1)

Initialize the CrossLinkedNetwork.

Parameters:

  • nFib (int) – Number of fibers

  • N (int) – Number of Chebyshev collocation points for fiber POSITION (note that this differs from the number of tangent vectors, which we typically refer to using \(N\)).

  • Nunisites (int) – Number of uniformly placed CL binding sites per fiber

  • Lfib (double) – Length of each fiber

  • kCL (double) – Cross linking spring constant (for the linear spring model)

  • rl (double) – Rest length of the CLs (for the linear spring model)

  • Dom (Domain object) – Periodic domain on which the calculation is being carried out

  • fibDisc (FibCollocationDiscretization object) – The discretization of each fiber’s centerline

  • smoothForce (boolean) – Whether to smooth out the forcing or use an actual spring between the uniformly spaced points. See the method CLForce for details on what formulas this switches between

  • nThreads (int) – Number of threads for openMP calculation of CL forces and stress

sigmaFromNL(N, L)

In the case when we smooth the fiber cross linking force (see method CLForce), the forcing is smoothed over a Gaussian with standard deviation \(\sigma\). This method gives \(\sigma\) as a function of the number of Chebyshev points \(N\) and length \(L\).

Parameters:

  • N (int) – Number of Chebyshev collocation points for fiber POSITION (note that this differs from the number of tangent vectors, which we typically refer to using \(N\)).

  • L (double) – Length of the fibers

Returns:

The smoothing lengthscale \(\sigma\)

Return type:

double

CLForce(uniPoints, chebPoints, Dom, fibCollection)

This method computes the total cross linking FORCE (not force density) in one of two ways. To establish some notation, let the link be between uniform point \(p\) on fiber \(X^{(i)}\) and uniform point \(q\) on fiber \(X^{(j)}\). Let the distance between the uniform points be \(R\).

If the variable self._smoothForce is set to true, then we first compute a force density

\[f^{(i)}(s) = -K \left(1-\frac{\ell}{R}\right) \delta_h(s-s_i^*)\int_0^L \left(X^{(i)}(s)-X^{(j)}(s')\right) \delta_h(s'-s_j^*) ds'\]
on fiber \(i\), with the corresponding formula for fiber \(j\), where \(\delta_h\) is a Gaussian density with standard deviation \(\sigma\). We then obtain force by multiplying by the \(L^2\) inner product weights matrix \(\widetilde{W}\). More details on this are in Section 9.1.3.2 in Maxian’s PhD thesis.

If, on the other hand, self._smoothForce is set to false, then we compute forces at the uniform points using the standard formulas for a spring, then map those forces to the Chebyshev points by multiplying by the corresponding row of the uniform resampling matrix. The exact formulas can be found in Section 9.1.3.1 in Maxian’s PhD thesis.

Parameters:

  • uniPoints (array) – Uniform points for all fibers

  • chebPoints (array) – Chebyshev points on all fibers

  • Dom (Domain object) – The (potentially strained) periodic domain on which we due the calculation

Returns:

nPts*3 one-dimensional array of the forces at the Chebyshev points due to the CLs

Return type:

array

CLStress(fibCollection, chebPts, Dom)

Compute the cross-linker contribution to the stress. If a single cross linker connects fibers \(i\) and \(j\), the formula for the stress due to this CL is given by

\[\sigma_{CL} = \sum_p X^{(i)}_p F^{(i)}_p + X^{(j)}_p F^{(j)}_p\]
where the product between \(X\) and \(F\) is an outer product and the positions \(X\) are the same periodic copies of the fibers that are linked by the CL.

At present, this method is only implemented for nonsmooth forcing (actual springs between the fibers). As such, it returns an error if you try to call it with smooth forcing.

Parameters:

  • fibCollection (fiberCollection object) – The object that stores the collection of fibers (and is queried to get the uniform points needed to compute forces)

  • chebPts (array) – Chebyshev points of all fiber locations

  • Dom (Domain object) – The periodic domain; needed to compute volume

Returns:

The stress due to the cross linkers as a \(3 \times 3\) array.

Return type:

Array

Returns:

The number of links connecting distinct fibers

Return type:

int

mapSiteToFiber(sites)
Parameters:

sites (array of ints) – Indices of the uniform sites

Returns:

The fiber(s) associated with a set of uniform site(s)

Return type:

int array

mapFibToStartingSite(LinkedFibs)

The first uniform site on a given fiber(s)

Parameters:

LinkedFibs (array of ints) – The fibers we are interested in

Returns:

The index of the first uniform site associated with the fibers in question

Return type:

int array

numLinksOnEachFiber()
Returns:

One-dimensional array of length nFibers, which gives the number of links bound to each fiber. The total number of links in the system is thus half the sum of this array

Return type:

array

calcLinkStrains(uniPoints, Dom)

Method to compute the strain in each link. We define (signed) strain as \(r-\ell\), where \(r\) is the length of the link at present and \(\ell\) is the rest length. This statistic is useful to tell if the simulation is unstable.

Parameters:

  • uniPoints (array) – Uniform points, as an nUni*nFib x 3 two-dimensional numpy array

  • Dom (Domain object) – Periodic domain

Returns:

An array of nLinks size that has the signed strain for each link

Return type:

array

SpringCOM(uniPoints, Dom)

Method to compute the center of mass of each link.

Parameters:

  • uniPoints (array) – Uniform points, as an nUni*nFib x 3 two-dimensional numpy array

  • Dom (Domain object) – Periodic domain

Returns:

An array of nLinks x 3 size that has the center of each link (average of the two linked points)

Return type:

array

Sort the links so that the site with smaller index comes first

ConnectionMatrix(bundleDist=0)

Build an nFib x nFib matrix that gives the fibers connected to each other.

Parameters:

bundleDist (double) – If bundleDist = 0, it will give you the adjacency matrix with nonzero entries at (i,j) if there is one link between fibers i and j. If bundleDist > 0, the nonzero entries (i,j) have 2 links between them separated by a distance (in arclength) at least bundleDist

Return type:

The fiber adjacency matrix

FindBundles(bunddist)

Uses the adjacency matrix of the previous method to sort the fibers into bundles based on the connectivity.

Parameters:

bunddist (double) – Used in the previous method to make the connectivity matrix

Returns:

The method first returns the number of bundles in the system, followed by an integer array of length nFib that records what bundle each fiber is in

Return type:

(int, int array)

BundleOrderParameters(fibCollection, nBundles, whichBundlePerFib, minPerBundle=2)

Get the order parameter of each bundle. The order parameter for a bundle \(B\) of F filaments is defined as the maximum eigenvalue of the matrix

\[M = \frac{1}{F L} \sum_{i \in B} \int_0^L {\tau^{(i)}(s)\tau^{(i)}(s) ds}\]
where we discretize the integral by direct Clenshaw-Curtis quadrature.

Parameters:

  • fibCollecton (fiberCollection object) – The collection of fibers

  • nBundles (int) – The number of bundles in the system (obtained from previous method)

  • whichBundlePerFib (int-array) – An integer array of length nFib that records what bundle each fiber is in (obtained from previous method)

  • minPerBundle (int) – The minimum number of fibers that qualify as a “bundle.” The default is two, so that we don’t compute order parameters for “bundles” of one fiber.

Returns:

Three arrays are returned: the order parameters (1D array), the number of fibers per bundle (1D array), and the average tangent vector of the bundle (NBundle x 3 array). These arrays are resized to remove “bundles” that have less than minPerBundle fibers.

Return type:

(array, array, array)

avgBundleAlignment(BundleAlignmentParams, nPerBundle)

Calculate the weighted average of alignment across all bundles,

\[\bar{q}= \frac{\sum q_i n_i}{\sum n_i}\]
where \(q_i\) is the alignment parameter of bundle \(i\) and \(n_i\) is the number of fibers in bundle \(i\).

Parameters:

  • BundleAlignmentParams (array) – The alignment parameters for each bundle, obtained from previous method.

  • nPerBundle (int array) – The number of fibers in each bundle, obtained from previous method.

Returns:

The average bundle alignment parameter \(\bar{q}\) defined above.

Return type:

double

LocalOrientations(nEdges, fibCollection)

Calculate the orientiation parameter for every fiber within nEdges of a given fiber. The order parameter is defined as the maximum eigenvalue of the matrix

\[M^{(i)} = \frac{1}{F L} \sum_{j \in N^{(i)}} \int_0^L {\tau^{(j)}(s)\tau^{(j)}(s) ds}\]
where we discretize the integral by direct Clenshaw-Curtis quadrature. Here the sum occurs over all fibers that are within a certain number of “edges” (links) in the adjacency matrix from fiber \(i\). We define this as the set of “neighbors” \(N^{(i)}\).

Parameters:

  • nEdges (int) – Number of conections we allow fibers to be separated by. For instance, if nEdges=1, then the sum is over all fibers connected to fiber \(i\). If nEdges=2, then the sum is over all fibers connected directly to fiber \(i\), or connected to an intermediate fiber that is connected to fiber \(i\). And so forth.

  • fibCollection (fiberCollection object) – The collecton of fibers we are operating on.

Returns:

A one-dimensional array of length nFib, where the \(i\)th entry is the orientation parameter, the maximum eigenvalue of the matrix \(M^{(i)}\) defined above.

Return type:

array

SparseMatrixOfConnections()
Returns:

Sparse matrix of size Nsite x Nsite, whose \((i,j)\) entry is the number of connections (links) between sites \(i\) and \(j\)

Return type:

Sparse matrix

nLinksAllSites(PossiblePairs)

Return the number of links between each pair of sites in the rows of input PossiblePairs

Set the links from input vectors of iPts, jPts, and Shifts

Parameters:

  • iPts (int array) – Heads of links

  • jPts (int array) – Tails of links

  • Shifts (three-dimensional array) – The periodic shift associated with each link

setLinksFromFile(FileName)

Set the links from a file name. The file has a list of iPts, jPts (two ends of the links), and shift in zero strain coordinates.

Parameters:

FileName (string) – The name of the file

Write the links to a file object of.

getPairsThatCanBind(fiberCol, Dom)

Generate list of events for binding of a doubly-bound link to both sites. SpatialDatabase = database of uniform points (sites). This is a method for debugging (not called in the final production code)

class DoubleEndedCrossLinkedNetwork.DoubleEndedCrossLinkedNetwork(nFib, N, Nunisites, Lfib, kCL, rl, kon, koff, konsecond, koffsecond, CLseed, Dom, fibDisc, kT=0, nThreads=1, bindingSiteWidth=0, smoothForce=True)

This class is a child of CrossLinkedNetwork which implements a network with cross links where each end matters separately.

There are 6 reactions (the even ones are the reverse of the odds):

  1. Binding of a floating link to one site (rate _kon)

  2. Unbinding of a link that is bound to one site to become free (reverse of 1, rate _koff)

  3. Binding of a singly-bound link to another site to make a doubly-bound CL (rate _konSecond)

  4. Unbinding of a double bound link in one site to make a single-bound CL (rate _koffSecond)

  5. Binding of both ends of a CL (rate _kDoubleOn. This rate is set to zero at present.)

  6. Unbinding of both ends of a CL (rate _kDoubleOff. This rate is also set to zero.)

For more information on the reactions and their rates, see Section 9.1.1 of Maxian’s PhD thesis. We implement reaction 3, which is the most important because it controls the forces exerted on the fibers, by performing a neighbor search to find all pairs of uniform sites a certain distance apart. Note also that we do not explicitly keep track of unbound cross linkers. We instead assume that there are always enough CLs if an event is supposed to happen.

We use an event driven algorithm to simulate a time step. To efficiently organize the events, we use a heap queue data structure which we borrow from the simulation package SRBD (https://github.com/stochasticHydroTools/SRBD). The heap queue is implemented in the fortran code /Fortran/MinHeapModule.f90 in this repository.

__init__(nFib, N, Nunisites, Lfib, kCL, rl, kon, koff, konsecond, koffsecond, CLseed, Dom, fibDisc, kT=0, nThreads=1, bindingSiteWidth=0, smoothForce=True)
Parameters:

  • nFib (int) – Number of fibers

  • N (int) – Number of Chebyshev collocation points for fiber POSITION (note that this differs from the number of tangent vectors, which we typically refer to using \(N\)).

  • Nunisites (int) – Number of uniformly placed CL binding sites per fiber

  • Lfib (double) – Length of each fiber

  • kCL (double) – Cross linking spring constant (for the linear spring model)

  • rl (double) – Rest length of the CLs (for the linear spring model)

  • kon (double) – Rate (units 1/(length x time)) at which a single end of a CL binds to a site

  • koff (double) – Rate (units 1/time) at which a singly-bound CL comes off from a site

  • konsecond (double) – Rate (units 1/(length x time)) at which the second end of a singly-bound CL binds to a site (this is modified by an Arrhenuis factor later – see the method updateNetwork)

  • koffsecond (double) – Rate (units 1/time) at which one end of a doubly-bound CL comes off, leaving a singly-bound CL

  • CLseed (int) – The seed for cross linking calculations

  • Dom (Domain object) – Periodic domain on which the calculation is being carried out

  • fibDisc (FibCollocationDiscretization object) – The discretization of each fiber’s centerline

  • kT (double) – Thermal energy. If zero, the binding of a second end occurs with rate konsecond (no Arrhenius factor). If nonzero, then the rate is modified (see updateNetwork method for details)

  • smoothForce (boolean, defaults to true) – Whether to smooth out the forcing or use an actual spring between the uniformly spaced points. See the method CLForce (CrossLinkedNetwork.py) for details on what formulas this switches between

  • nThreads (int) – Number of threads for openMP calculation of CL forces and stress

  • bindingSiteWidth (double, optional (defaults to 0)) – If a binding site in biology has a width \(w\), then the number of links that can bind to one of our uniformly-spaced binding sites (spaced \(\Delta s_u\)) is given by \(N = \lceil \Delta s_u/w \rceil\). This in practice becomes a cap on the number of links per site.

sitesPerFib(iFib)
Parameters:

iFib (int) – Fiber index

Returns:

The indices of the uniform binding sites that correspond to the fiber iFib

Return type:

int array

updateNetwork(fiberCol, Dom, tstep)

Update the network using Kinetic Monte Carlo. The full procedure we use is documented in Section 9.1.1 of Maxian’s PhD thsis. Briefly, we consider six possible reactions, with the following rates in units 1/time:

  1. Binding of a floating link to one site (rate _kon \(\Delta s_u\))

  2. Unbinding of a link that is bound to one site to become free (reverse of 1, rate _koff)

  3. Binding of a singly-bound link to another site to make a doubly-bound CL. The rate of this reaction, assuming \(k_B T > 0\) is given by

    \[k_{on,s} = k_{on,s}^0 \exp{\left(\frac{-K_c (r-\ell)^2}{k_B T}\right)} \Delta s_u\]
    If \(k_B T = 0\), then the Arrhenius factor is ommitted. If \(k_B T > 0\), we assume possible connections can form if two sites are separated by \(2\sqrt{k_B T/K}\), where \(K\) is the CL stiffness (two standard deviations of the Gaussian above). If \(k_B T =0\) and the binding probability is uniform, we stick to one standard deviation.

  4. Unbinding of a double bound link in one site to make a single-bound CL (rate _koffSecond)

  5. Binding of both ends of a CL (rate _kDoubleOn. This rate is set to zero at present.)

  6. Unbinding of both ends of a CL (rate _kDoubleOff. This rate is also set to zero.)

Once the rates are computed, the time at which an event occurs is taken from an exponential distribution via sampling a time \(t=-\ln(1-u)/k\), where \(k\) is the event rate and \(u\) is a random draw from \(U(0,1)\). These times are organized into a heap queue which efficiently sorts them.

Parameters:

  • fiberCol (fiberCollection object) – The fibers we are simulating

  • Dom (Domain object) – The periodic (sheared?) domain we are simulating

  • tstep (double) – Time step we simulate over

Returns:

Nothing, but updates the state of the network internally. The actual update is done in C++. On the python side, we get the neighbors and filter them so that fibers cannot link to themselves. Then we pass the list of potential neighbors to C++ which updates its own chain. The Python side is then synced with the C++ side.

Return type:

Null

Set the links from input vectors of iPts, jPts, and Shifts

Parameters:

  • iPts (int array) – Heads of links

  • jPts (int array) – Tails of links

  • Shifts (three-dimensional array) – The periodic shift associated with each link

  • FreelyBound (array, optional) – This gives the number of singly-bound links bound to each of the uniform sites. If not supplied, it will set the number to zero for all sites.

setLinksFromFile(FileName, FreelyBound=None)

Set the links from a file name. The file has a list of iPts, jPts (two ends of the links), and shift in zero strain coordinates.

Parameters:

  • FileName (string) – The name of the file

  • FreelyBound (array, optional) – This gives the number of singly-bound links bound to each of the uniform sites. If not supplied, it will set the number to zero for all sites.

deleteLinksFromFibers(fibNums)

Removes all links (both double and singly bound) connected to a fiber

Parameters:

fibNums (list) – List of fiber numbers to remove links from

syncPythonAndCpp()

Synchronizes the python and C++ codes by copying the C++ over to python