single_cell_minc_model.py

from shapely.geometry import Polygon, LineString, MultiLineString
from matplotlib import pyplot as plt
import numpy as np

from scipy import interpolate


def get_minc_division_prediction(boundaryVerts, beta=1, numMTs=100,
                                 spindleLen='default', plot=1):
    """ Implement "Minc" spindle orientation model.
    Params:

    boundaryVerts (list): A list of (x,y) coords representing the cell boundary pixels.

    beta (double): Force-length power law scaling for MT forces: f = l^beta, where l=length of MT.

    numMTs (int): How many MTs are placed around the cell (uniformly distributed).

    spindleLen (double/string): Size of the spindle. (distance between centrosomes). Can manually set a double, or pass "default" to use sqrt(area/pi).

    plot (bool): To plot results, or not.
    """

    # Get the list of angles to distribute MTs
    spindleAngles = np.linspace(-np.pi/2, np.pi/2, 100)
    # Define how far the initial MT lines should extend (just needs to cross the
    # cortex so we can see the distance to crossing.
    MTexpansionLen = 100
    numMTs = 100
    # # Make some generic cell shapes for testing:
    # # polygon = Polygon([(0, 0), (1, 0), (1.5, 0.5), (1, 1), (0, 1)])
    # # polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]) # Square
    # polygon = Polygon([(0, 0), (2, 0), (2, 1), (0, 1)]) # Rectangle
    # polygon = Polygon([(0, 0), (2, 0), (1, 1)]) # Triangle
    # # Hexagons
    # # Stretched
    # xs = 1 + np.cos(np.linspace(0, 2*np.pi, 100))
    # ys = 1 + 0.5*np.sin(np.linspace(0, 2*np.pi, 100))
    # # Regular
    # xs = np.cos((np.pi * np.array([0, 1, 2, 3, 4, 5]))/3)
    # ys = np.sin((np.pi * np.array([0, 1, 2, 3, 4, 5]))/3)
    # polygon = Polygon([i for i in zip(xs, ys)])

    # Polygon geometry:
    polygon = Polygon(boundaryVerts)
    centroid = polygon.centroid.coords[:][0]
    area = polygon.area

    if spindleLen == 'default':
        spindleLen = np.sqrt(area/np.pi)  # Set spindle len for circle
    # spindleLen = 1.95

    torqueList = []
    for spindleAngle in spindleAngles:

        # Centrosome position
        cSome1 = (centroid[0] - 0.5*spindleLen*np.cos(spindleAngle),
                  centroid[1] - 0.5*spindleLen*np.sin(spindleAngle))
        cSome2 = (centroid[0] + 0.5*spindleLen*np.cos(spindleAngle),
                  centroid[1] + 0.5*spindleLen*np.sin(spindleAngle))

        # PLOTTING
        if plot:
            fig = plt.figure(1)
            ax = fig.add_subplot(111)
            ax.cla()

            # Centrosomes
            ax.plot([cSome1[0], cSome2[0]], [cSome1[1], cSome2[1]], color='tomato',
                    alpha=0.8, linewidth=3, solid_capstyle='round', zorder=2)
            ax.plot(cSome1[0], cSome1[1], 'o', color='tomato',
                    alpha=0.8, linewidth=3, solid_capstyle='round', zorder=3)
            ax.plot(cSome2[0], cSome2[1], 'o', color='tomato',
                    alpha=0.8, linewidth=3, solid_capstyle='round', zorder=3)
            # Outline
            x, y = polygon.exterior.xy
            ax.plot(x, y, color='#6699cc', alpha=0.7,
                    linewidth=3, solid_capstyle='round', zorder=2)
            ax.set_title('Single Cell Spindle')

        # TODO() Csome1 needs to be to the left of cSome2
        # Loop over a range of angles to get sum of MT forces at each angle.
        torque1, torque2 = 0, 0
        MTangleList = np.linspace(-np.pi/2, np.pi/2, numMTs)
        for phi in MTangleList:

            angle = phi + spindleAngle
            #############################
            # Create a MT
            MT1 = LineString([cSome1, (cSome1[0] - MTexpansionLen*np.cos(angle),
                                       cSome1[1] - MTexpansionLen*np.sin(angle))])
            # Get the length of the microtubule from centrosome to the cortex.
            intersectLine = MT1.intersection(polygon)
            # If the cell is highly convex, there may be multiple crosses
            # but we just wat the first. So check and take the first if so.
            _multiString = MultiLineString()
            if type(intersectLine) == type(_multiString):
                intersectLine = intersectLine[0]

            # Get the torque, if there was a line.
            if not intersectLine.is_empty:
                MTLen = intersectLine.length

                # Store the torque generated by the MT
                torque1 += 0.5*spindleLen*np.sin(phi)*(MTLen**beta)

                if plot:
                    # Plot MTs
                    x1, y1 = intersectLine.xy
                    ax.plot(x1, y1, alpha=0.7, color='limegreen',
                            linewidth=1, solid_capstyle='round', zorder=2)
            else:
                torque1 += np.nan

            #############################

            # Create a MT
            MT2 = LineString([cSome2, (cSome2[0] + MTexpansionLen*np.cos(angle),
                                       cSome2[1] + MTexpansionLen*np.sin(angle))])
            # Get the length of the microtubule from centrosome to the cortex.
            intersectLine = MT2.intersection(polygon)
            _multiString = MultiLineString()
            if type(intersectLine) == type(_multiString):
                intersectLine = intersectLine[0]

            if not intersectLine.is_empty:
                MTLen = intersectLine.length

                # Store the torque generated by the MT
                torque2 += 0.5*spindleLen*np.sin(phi)*(MTLen**beta)

                if plot:
                    x1, y1 = intersectLine.xy
                    ax.plot(x1, y1, alpha=0.7, color='limegreen',
                            linewidth=1, solid_capstyle='round', zorder=2)
            else:
                torque2 += np.nan

        # Update the torque
        torque = torque1+torque2
        torqueList.append(torque)

        if plot:
            plt.axis('equal')
            plt.tight_layout()

    # Plot the torques vs angle
    if plot:
        fig = plt.figure(2)
        ax = fig.add_subplot(111)
        ax.plot(spindleAngles, torqueList)
        ax.set_xlabel('Spindle Angle')
        ax.set_ylabel('Torque')
        ax.set_title("Torque vs Spindle Angle")
        plt.tight_layout()

        fig = plt.figure(3)
        ax = fig.add_subplot(111)

    torqueList = np.array(torqueList)

    # Add nan to the MTangleList so we can remove same way as torques
    MTangleList[np.isnan(torqueList)] = np.nan
    # Split arrays up between nan entries
    torqueListSplit = np.split(torqueList, np.where(np.isnan(torqueList))[0])
    MTangleListSplit = np.split(MTangleList, np.where(np.isnan(MTangleList))[0])
    # removing NaN entries
    torqueListSplit = [ev[~np.isnan(ev)] for ev in torqueListSplit
                       if not np.isnan(ev).all()]
    MTangleListSplit = [ev[~np.isnan(ev)] for ev in MTangleListSplit
                        if not np.isnan(ev).all()]
    # removing empty DataFrames
    torqueListSplit = [ev for ev in torqueListSplit if ev.size > 1]
    MTangleListSplit = [ev for ev in MTangleListSplit if ev.size > 1]

    # For each sublist, get the minima.
    minima = []
    for i in range(0, len(torqueListSplit)):
        if len(torqueListSplit[i]) > 2:
            # INterpolate a spline
            spline = interpolate.InterpolatedUnivariateSpline(
                MTangleListSplit[i], torqueListSplit[i])
            # Get the roots
            xp = spline.roots()
            # minima.extend(xp)
            # Check if the torque was decreasing, indicating a minimum.
            tempMinima = []
            for root in xp:
                if spline(root-.01) > spline(root+.01):
                    tempMinima.append(root)

            # Get the derivativeself
            deriv = spline.derivative()
            # If there were multiple roots, choose the smallest minimum
            if len(tempMinima) > 1:
                # Get values of deriv:
                derivVals = deriv(tempMinima)
                # Select the one with smallest deriv.
                minima.append(tempMinima[np.argmin(derivVals)])
            elif len(tempMinima) == 1:
                minima.append(tempMinima[0])

            if plot:
                xs = np.linspace(MTangleListSplit[i][0],
                                 MTangleListSplit[0][-1], 1000)
                # plt.plot(MTangleListSplit[i], torqueListSplit[i])
                plt.plot(xs, spline(xs))
                plt.plot(xp, [0]*xp.size, 'o')
                plt.plot(xs, deriv(xs))

    if plot:
        plt.tight_layout()

    # If there were no minima, just get the smallest torque
    if len(minima) == 0:
        MTangleList = MTangleList[~np.isnan(torqueList)]
        torqueList = torqueList[~np.isnan(torqueList)]
        index = np.argmin(np.abs(torqueList))
        minima = [MTangleList[index]]
    # If there were multiple minima, then just choose the smallest?
    if len(minima) > 1:
        print("Multiple minima")
    minimum = minima[0]

    if plot:
        plt.show()

    # This has found the angle of the spindle, now we return the predicted
    # Divison angle.
    minimum = np.degrees(minimum)
    print(minimum)

    return minimum
    # divAngle = minimum + 90
    # # if minimum < 0:
    # #     minimum += 90
    # # else:
    # #     minimum -= 90
    #
    # return divAngle