# Exercise session nº 5
---
# Furrow Constriction in Animal Cell Cytokinesis

__*Sacha Ichbiah, 21/02/22, ENS Paris*__

This subject is extracted from : 
> Hervé Turlier et al., *Furrow Constriction in Animal Cell Cytokinesis*, Biophysical Journal, 2014. \
> https://doi.org/10.1016/j.bpj.2013.11.014

Cytokinesis is the process of physical cleavage at the end of cell division; it proceeds by ingression of an acto- myosin furrow at the equator of the cell. Its failure leads to multinucleated cells and is a possible cause of tumorigenesis. Despite its ubiquity in developmental biology, its precise description and understanding have challenged biologists and physicists.

In this paper, the authors propose a model based on a minimal geometry and scaling arguments that gives an interpretation of the process appearing during cytokinesis. It notably demonstrates that because of the cytoplasm incompressibility, the cytokinesis leads to a competition between the furrow line tension and the cell poles' surface tension. This competition sets a threshold for cytokinesis completion, and explains cytokinesis dynamics. 



During these session, we will derive the equations of this scaling model of furrow constriction, and will integrate these equation to see constriction dynamics. We will show that it allows to have a cytokinesis duration independant of cell size, which has been observed in C-Elegans.



---

## I - The Scaling Model

The geometry of the dividing cell is described by the apposition of two spherical caps, parametrized by an angle $\theta$ as shown on the left sketch.

The volume of a spherical cap (in blue) is : $\mathcal{V}_{sc}(r,h) = \dfrac{\pi}{3} h^2 (3r - h)$, and its area : $\mathcal{A}_{sc}(r,h)=2\pi r h $ (right sketch).



#### **Question 1 :**
> Noting that the cytoplasmic is an incompressible fluid, establish that $R_0 = R F(\theta)$

We define a dimensionless parameter $\kappa$ to express the competition between the mean contractile surface tension at the furrow and the tension at the cell poles: $\kappa = \dfrac{\gamma}{2R_0N^a_0}$.

The polar contractility tends to reduce the surface $A_p = 2\pi R^2 (1+\text{cos}(\theta))$ of each cell poles, whereas the line tension tends to reduce the contractile ring circumference $r_f$. These effects are captured by a simple mechanical energy $\mathcal{E} = 2\pi r_f \gamma + 2 A_p N^a_0$.

#### **Question 2 :**
> Rescale the energy $\mathcal{E}$ by an energy $\mathcal{E}_0 = 4 \pi R_0^2 N^a_0$ to make it only depend on $\theta$ and $\kappa$.


## II - Mechanical Equilibrium

The local minimum of the energy gives the equilibrium configuration of the cell. 

To find this minimum, we will use a library doing symbolic calculus in python called sympy. This will allow us to compute the derivatives effortlessly.



### Symbolic Computation with Sympy *(In french "Calcul Formel")*

We will use sympy, a library that allows to do symbolic computation. Analytical results are always the best, but sometimes the equations does not lead to beautiful simplifications. If we are interested in the numerical result of the equations, we can use sympy to work on the analytical expression directly, obtain derivatives, etc.. before evaluating them on real values. There are three main functions that we will use in sympy, that we will present briefly. If interested, the best symbolic calculus tools are present in Wolfram Mathematica, which is under license. 

#### a) Defining expression with symbols and trigonometric functions, and obtain derivatives : 


```python
!pip install sympy

import sympy
from sympy import symbols, diff,lambdify, simplify
from sympy import cos, sin
import numpy as np 
from scipy.optimize import minimize
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
```


```python
#Symbols are unknown variables, that we aim to replace with real values in the end. We define them with the function symbols : 
a,b,c = symbols("a b c")

#We can then define a symbolic expression, containing or not trigonometric functions (among many other possibilities !)
E = a**2 + a*b + cos(c)

#And obtain its derivatives with respect to any variables, eg a : 
First_derivative = E.diff(a)
Second_derivative = First_derivative.diff(a)
First_derivative, Second_derivative
```




    (2*a + b, 2)



#### b) Substituting variables and evaluating symbolic expressions : 


```python
# We can replace symbols with real variables with the method subs : 
print("c = pi gives :",E.subs([(c,np.pi)]))
print("Subs method : ",E.subs([(a,2),(b,1),(c,0)]) )

#We can also transform a symbolic expression into a lambda function
#This is a faster process than subs if we need to evaluate the function on many points : 
f = lambdify((a,b,c),E,"numpy")
print("Lambify method : ",f(2,1,0))

#We can combine both to replace certain variables before creating a lambda function with the remaining variables : 
g = lambdify(a,E.subs([(b,1),(c,0)]))
print("Subs and lambdify combined :",g(2))

#Short benchmarking 
from time import time
values_evaluated = np.linspace(0,np.pi,1000)
t1 = time()
g = lambdify(a,E.subs([(b,2),(c,0.2)]))
g(values_evaluated)
t2 = time()
for value in values_evaluated : 
    E.subs([(a,value),(b,2),(c,0.2)])
t3 = time()

print("Time with lambdify :",round((t2-t1),4))
print("Time with subs :",round((t3-t2),4))
```

    c = pi gives : a**2 + a*b - 1.0
    Subs method :  7
    Lambify method :  7.0
    Subs and lambdify combined : 7
    Time with lambdify : 0.0029
    Time with subs : 0.9684


### The equilibrium configuration during cytokinesis

Let's go back to the initial problem. Our goal is to study the properties of the normalized energy $\overline{\mathcal{E}} = \mathcal{E}/\mathcal{E}_0$.


```python
x, k = symbols("x k")
F = 1 +1.5*cos(x) - 0.5*(cos(x))**3
energy = k * sin(x)/(F**(1/3)) + (1+ cos(x))/(F**(2/3))

#We see that there is no simplification easily given by sympy : 
print(simplify(energy),'\n')

#We can replace the values of k and x : 
print(energy.subs([(x,np.pi/3),(k,.5)]))

energy
```

    (0.873580464736299*k*(-0.333333333333333*cos(x)**3 + cos(x) + 0.666666666666667)**0.666666666666667*sin(x) + 0.763142828368888*(cos(x) + 1)*(-0.333333333333333*cos(x)**3 + cos(x) + 0.666666666666667)**0.333333333333333)*(-0.333333333333333*cos(x)**3 + cos(x) + 0.666666666666667)**(-1.0) 
    
    1.42197524663604





$\displaystyle \frac{0.873580464736299 k \sin{\left(x \right)}}{\left(- 0.333333333333333 \cos^{3}{\left(x \right)} + \cos{\left(x \right)} + 0.666666666666667\right)^{0.333333333333333}} + \frac{0.763142828368888 \left(\cos{\left(x \right)} + 1\right)}{\left(- 0.333333333333333 \cos^{3}{\left(x \right)} + \cos{\left(x \right)} + 0.666666666666667\right)^{0.666666666666667}}$



Now that we have implemented our energy in sympy, we can automatically obtain the derivatives with the diff engine. We see that the analytical formulas are quite long and obtaining by hand the derivatives would be both painful and prone to errors.

#### **Question 3 :**
> Obtain the expression of the first and second derivatives of the energy $\overline{\mathcal{E}}$ with respect to theta (i.e x) with the diff function :


```python

```

#### **Question 4 :**
> Plot the energy profile $\overline{\mathcal{E}}(\theta), \theta \in [0,\dfrac{\pi}{2}]$ for $\kappa \in \{0.0,0.1,0.2,0.3,0.4,0.5\}$. What do you observe ?


```python
plt.figure(figsize = (15,25))
vals_theta = np.linspace(0,np.pi/2,10000)

for k_val in np.linspace(0,0.5,6):
    print("k_value :",k_val)
    ...
    plt.plot(vals_theta,...,label = k_val.round(2))
plt.legend()
```

#### **Question 5 :**
> Starting with $\theta = \pi/4$, find the angle giving the local energy minimum for each k in values_k. Plot the equilibrium angle $\theta_{min}$, and the value of the derivatives of the energy $\left. \dfrac{\partial \mathcal{E}}{\partial \theta}\right\rvert_{\theta_{min}}$, $\left. \dfrac{\partial^2 \mathcal{E}}{\partial \theta^2}\right\rvert_{\theta_{min}}$, $\left (\left. \dfrac{\partial^2 \mathcal{E}}{\partial \theta^2}\right\rvert_{\theta_{min}}\right)^2$ at this angle for each k. 

__Tip :__ : Lambdify the energy function into a function f and find its minimum with the function minimize from scipy (see bellow)


```python
e = lambdify((x,k), energy, "numpy")
npoints = 1001
values_k = np.linspace(0, 1, npoints)
...

for j,k_val in enumerate(values_k) : 
    
    f = lambda x : e(x,k_val)
    
    ...
    
    sols = minimize(fun = f,x0 = (np.pi/4),method = "SLSQP",bounds=[(0,np.pi/2)])
    assert (sols.success)
    min_theta = sols['x']
    
    ...

    
fig,ax = plt.subplots(1,4,figsize = (21,3))
ax[0].plot(values_k, ...)
ax[0].set_xlabel("k")
ax[0].set_ylabel("minimum found")
ax[1].plot(values_k, ...)
ax[1].set_xlabel("k")
ax[1].set_ylabel("derivative value at minimum")
ax[2].plot(values_k, ...)
ax[2].set_xlabel("k")
ax[2].set_ylabel("second derivative value at minimum")
ax[3].plot(values_k, ...)
ax[3].set_xlabel("k")
ax[3].set_ylabel("square of second derivative value at minimum")
```

#### **Question 6 :**
> Estimate the value $k_c$ where this local minimum disappears, and its associated angle $\theta_c$. Compute the values of the first two derivatives $\left. \dfrac{\partial \mathcal{E}}{\partial \theta}\right\rvert_{\theta_c} $ and  $\left. \dfrac{\partial^2 \mathcal{E}}{\partial \theta^2}\right\rvert_{\theta_c} $. What do you deduce on the order of the phase transition ? 

__Tip :__ Find the value where the second derivative evaluated at $\theta_c$ is equal to 0 (i.e the minimum of the square of the second derivative evaluated at $\theta_c$)

#### **Optional : Question 7 :**
> Plot the equilibrium angles and draw the shape of the cells for $\kappa \in \{0,0.2,0.4,0.6\}$


```python
fig,ax = plt.subplots(1,5,figsize =(25,5))
colors = ['tab:blue','tab:orange','tab:green','tab:red']
R0 = 1
for j,idx_val in enumerate([0,2000,4000,6000]) : 
    
    theta_sol = Solutions[idx_val]
    k_value = values_k[idx_val]
    R = R0/((F.subs(x,theta_sol[0]))**(1/3))
    e = lambdify((x), energy.subs(k,k_value), "numpy")
    
    theta_values_k = np.linspace(theta_sol, 2*np.pi-theta_sol,100)
    circle_x = R*np.cos(theta_values_k)
    circle_y = R*np.sin(theta_values_k)

    ax[0].plot(vals_theta*180/np.pi,e(vals_theta))
    ax[0].scatter(theta_sol*180/np.pi,e(theta_sol),s = 180)
    ax[0].set_ylabel("Energy")
    ax[0].set_xlabel("Angle value")
    ax[0].set_xlim(-5,95)

    
    ax[j+1].plot(circle_x-R*np.cos(theta_sol),circle_y,color=colors[j],linewidth=5)
    ax[j+1].plot(R*np.cos(theta_sol)-circle_x,circle_y,color=colors[j],linewidth=5)
    ax[j+1].set_title("Equilibrium angle value :" + str((theta_sol[0]*180/np.pi).round(2)))
    ax[j+1].set_xlim(-2,2)
    ax[j+1].set_ylim(-2,2)
    ax[j+1].set_aspect('equal')#, adjustable='box')
```

## III - Dynamics

We now want to study the furrow constriction dynamics, i.e the temporal evolution of $\dfrac{r_f}{R_0}$. We will establish these dynamics by expressing the derivative of this quantity with respect to $\theta$. As before, we will use the symbolic computation library sympy to evaluate numerical quantities.

To establish the dynamic equation, we note that the power of active effects is exactly dissipated by viscous cell deformations. The viscous dissipation is made of two contributions, the stretching of the poles and the constriction of the ring, which we estimate in scaling. The volume of acto-myosin in the poles is $V_p = 2A_p e_p$ and in the ring $V_f = 2\pi r_f w e_f$, where $w$ and $e_f$ are the width and thickness of the contractile ring. (We remind that the surface of each cell pole writes : $A_p = 2\pi R^2 (1+\text{cos}(\theta))$). The value $e_p \approx e_0$ and the ring thickness $e_f$ reach a steady-state value that depends on turnover. This yields the viscous dissipated power : 



$\begin{align}
P_d &= \dfrac{1}{2} \eta \left[ V_p \left(\dfrac{1}{R} \dfrac{dr_f}{dt} \right)^2 + V_f \left( \dfrac{1}{r_f} \dfrac{dr_f}{dt}\right)^2  \right] \newline
P_d &= \dfrac{1}{2} \eta \left[ e_f 4 \pi R^2 (1+\text{cos}\theta) \dfrac{1}{R^2}  \left(\dfrac{dr_f}{dt} \right)^2 + 2\pi w e_f r_f \dfrac{1}{r_f^2} \left(  \dfrac{dr_f}{dt}\right)^2  \right]  \newline
&\approx \dfrac{1}{2} \eta \left[ 4 e_0 \pi(1+\text{cos}\theta)  \left(\dfrac{dr_f}{dt} \right)^2 + 2\pi w e_f\dfrac{1}{r_f} \left( \dfrac{dr_f}{dt}\right)^2  \right]  \newline
&= \left(\dfrac{dr_f}{dt} \right)^2 \dfrac{1}{2} \eta \left[ 4 \pi e_0 (1+\text{cos}\theta)   + \dfrac{4\pi}{2} w e_f\dfrac{F(\theta)^{1/3}}{R_0 \sin \theta}  \right]  \newline
&= \left(\dfrac{dr_f}{dt} \right)^2 4 \pi e_0 \eta \left[ (1+\text{cos}\theta)   + \dfrac{1}{2 R_0} w \dfrac{e_f}{e_0} \dfrac{F(\theta)^{1/3}}{\sin \theta}  \right]  \newline
&= \left(\dfrac{dr_f}{dt} \right)^2 4 \pi e_0 \eta \left[ (1+\text{cos}\theta)   + \lambda \dfrac{F(\theta)^{1/3}}{\sin \theta}  \right]  \newline
\end{align}
$

The balance of mechanical and dissipated powers yields : 

$\dfrac{d \mathcal{E}}{dt} + P_d = 0$

Besides :

$
\dfrac{1}{\mathcal{E}_0} \dfrac{d\mathcal{E}}{dt} = \dfrac{\partial \mathcal{E}/\mathcal{E}_0}{\partial \theta} \dfrac{\partial \theta}{\partial r_f} \dfrac{d r_f}{d_t} = \dfrac{\partial \mathcal{E}/\mathcal{E}_0}{\partial \theta} \left(\dfrac{\partial r_f}{\partial \theta}\right)^{-1} \dfrac{d r_f}{d_t} 
$

And, with $T_a=\frac{\eta e_0}{N^a_0}$ : 

$\dfrac{1}{\mathcal{E}_0} P_d = \left(\dfrac{dr_f}{dt} \right)^2 \dfrac{4\pi e_0 \eta}{4 \pi R_0^2 N^a_0} \left[ (1+\text{cos}\theta)   + \dfrac{\lambda}{2} \dfrac{F(\theta)^{1/3}}{\sin \theta}  \right] = \left(\dfrac{dr_f}{dt} \right)^2 \dfrac{T_a}{ R_0^2} \left[ (1+\text{cos}\theta)   + \dfrac{\lambda}{2} \dfrac{F(\theta)^{1/3}}{\sin \theta}  \right] $

We have thus from $\dfrac{d \mathcal{E}/\mathcal{E}_0}{dt} = - P_d/\mathcal{E}_0$ : 

$
\begin{align}
\dfrac{\partial \mathcal{E}/\mathcal{E}_0}{\partial \theta} \left(\dfrac{\partial r_f}{\partial \theta}\right)^{-1} \dfrac{d r_f}{d_t}  &= - \left(\dfrac{dr_f}{dt} \right)^2 \dfrac{T_a}{ R_0^2} \left[ (1+\text{cos}\theta)   + \dfrac{\lambda}{2} \dfrac{F(\theta)^{1/3}}{\sin \theta}  \right] 
\newline
\dfrac{\partial \mathcal{E}/\mathcal{E}_0}{\partial \theta} \left(\dfrac{\partial r_f/R_0}{\partial \theta}\right)^{-1} &= - \dfrac{dr_f}{dt} \dfrac{T_a}{ R_0} \left[ (1+\text{cos}\theta)   + \dfrac{\lambda}{2} \dfrac{F(\theta)^{1/3}}{\sin \theta}  \right] \newline 
\dfrac{dr_f}{dt} \dfrac{T_a}{ R_0} &= - \dfrac{\partial \mathcal{E}/\mathcal{E}_0}{\partial \theta} \left(\dfrac{\partial r_f/R_0}{\partial \theta}\right)^{-1} \left[ (1+\text{cos}\theta)   + \dfrac{\lambda}{2} \dfrac{F(\theta)^{1/3}}{\sin \theta}  \right]^{-1} = -\mathcal{H}(\theta,\kappa,\lambda)
\end{align}
$

We will compute numerical the values of this function $\mathcal{H}$ to obtain the evolution of the furrow radius $r_f$

#### **Question 8 :**
> From the last equation, express the angle temporal variation $\dot \theta$. 

#### **Question 9 :**
> Compute numerically $\dot \theta$ with sympy, and integrate the evolution of $\theta$, $r_f$ in time with a forward-euler-scheme for $t \in [0,15]$, starting with $\theta(0)=\pi/2$, with $\lambda = 1$ and $\kappa \in \{0.1,0.25,0.4,0.5,0.75,1\}$. Check that it is compatible with the previous results obtained from the static analysis.



```python
x, k, l = symbols("x k l")
F = 1 +1.5*cos(x) - 0.5*(cos(x))**3
energy = k * sin(x)/(F**(1/3)) + (1+ cos(x))/(F**(2/3))

r_f = sin(x)*(F**(-1/3))

dr_f = diff(r_f,x)
first_derivative = diff(energy,x)

H = ((1+cos(x)) + l*(F**(1/3))/sin(x))**(-1) * first_derivative/dr_f

dtheta = - H/dr_f
dtheta = dtheta.subs(l,0.1)

#We can already lambdify the function r_f as it does not depend on k : 
func_rf = lambdify(x,r_f,"numpy")
```


```python
theta0 = np.pi/2
R0 = 1

for k_value in [0.1,0.25,0.4,0.5,0.75,1]:

    ...

    npoints = 1000
    timepoints = np.linspace(0,15,npoints)
    dt = timepoints[1]-timepoints[0]
    
    ...


    for j,t in enumerate(timepoints) :

        ...

    plt.plot(timepoints,...,label = k_value)
    
plt.title("Constriction completion or failure")
plt.xlabel("time t/Ta")
plt.ylabel("furrow radius r_f/R_0")
plt.legend()
```

#### **Question 10 :**
> Determine the cytokinesis duration with $R_0 \in \{0.5,1,2,4\}$. Show that in case of cytokinesis completion $\lambda = 1, \kappa = 0.75$, the initial cell radius R0 has no impact on the cytokinesis time.


```python
theta0 = np.pi/2
k_value = 0.75

for R0 in [0.5,1,2,4]:

    ...

    npoints = 1000
    timepoints = np.linspace(0,6,npoints)
    dt = timepoints[1]-timepoints[0]
    
    ...

    for j,t in enumerate(timepoints) :

        ...


    plt.plot(timepoints,...,label ="R0: " +str(R0))
    
plt.title("Cytokinesis time is independant of the cell initial radius")
plt.xlabel("time t/Ta")
plt.ylabel("furrow radius r_f/R_0")
plt.legend()
```

981121f584e2e30ca2961318f5cf173510a6e0ef
ipynb
Jupyter Notebook
Ichbiah/TD_5-Cytokinesis/TD_5_Cytokinesis.ipynb
hturlier/M2ICFP
1d91ff837b05a6058ee34a03fdc8062893287c6e
MIT
2022-02-14T10:17:11.000Z
2022-03-22T21:16:42.000Z
Ichbiah/TD_5-Cytokinesis/.ipynb_checkpoints/TD_5_Cytokinesis-checkpoint.ipynb
hturlier/M2ICFP
1d91ff837b05a6058ee34a03fdc8062893287c6e
MIT
Ichbiah/TD_5-Cytokinesis/.ipynb_checkpoints/TD_5_Cytokinesis-checkpoint.ipynb
hturlier/M2ICFP
1d91ff837b05a6058ee34a03fdc8062893287c6e
MIT
2022-01-24T15:08:21.000Z
2022-02-14T10:17:00.000Z
Qwen/Qwen-72B
1. YES
2. YES