Laboratoire RDP » Biophysics Team » Cell Division Rule

#include <math.h>

#include "cell_division.h"

#include "constants.h"

#include <stdio.h>

#include <stdlib.h>

#include <omp.h>

/* this function:   finds the division direction

                    computes the position of the new vertices

                    saves the indices and the vertices of the two daughter cells and of the two walls which are cut by the new one (these data are usefull for cell_division2.cpp) */

void cell_division(double *vertices, int vertices_number, int *cells_vertices, int *vertices_number_in_each_cell, int cell_number,

          double *barycenters, int dividing_cell, int *boundary, double *stress, double *shapes,

          int *adjacent_cells_number, int *adjacent_cell1, int *adjacent_cell2, int *daughter_cells_vertices_number, int *cell1, int *cell2, int *wall1, int *wall2){

///bias division

//double width=3; double dilation=1.;

int j;

/* vector of the division line */

double division_direction[2];

/* positions of new vertices */

double new_vertex1[2]={0}; double new_vertex2[2]={0};

double eigenvalue=0; double norm;

if (threshold_division>0.0){ /* probabilistic division: the probability to divide a cell following the vector n is proportionnal to

                        exp(n^t D n / D0) where D is the deviatoric stress (S-tr(S)/2) and D0 the threshold over which we start felling the stress */

    double anisotropic_stress[2];

    anisotropic_stress[0]=(0.5*stress[3*dividing_cell]-0.5*stress[3*dividing_cell+1])/threshold_division;

    anisotropic_stress[1]=0.5*stress[3*dividing_cell+2]/threshold_division;

    double probability[20];

    for (int i=0;i<20;i++){

        probability[i]=exp(cos(i*pi/19.0)*anisotropic_stress[0]+sin(i*pi/19.0)*anisotropic_stress[1]);

        if (probability[i]==HUGE_VAL){ printf("division_h\n"); goto division_highest_tension;}

    }

    probability[0]/=2.; probability[19]/=2.;

    for (int i=1;i<20;i++){probability[i]+=probability[i-1];}

    for (int i=0;i<20;i++){probability[i]/=probability[19];}

    double random=(double)rand()/(double)RAND_MAX; //printf("random %f\n",random);

    int random2=rand()%2;

    if (random<=probability[0]){ //printf("%i\n",0);

        division_direction[0]=1;

        division_direction[1]=0;}

    for (int i=1;i<20;i++){//printf("%f %f\n",probability[i-1],probability[i]);

        if (random<=probability[i] && random>probability[i-1]){

            if(random2==0){//printf("%i %i\n",i,0);

                division_direction[0]=cos(i*0.5*pi/19.0);

                division_direction[1]=sin(i*0.5*pi/19.0);}

            if(random2==1){//printf("%i %i\n",i,1);

                division_direction[0]=cos(i*0.5*pi/19.0);

                division_direction[1]=-sin(i*0.5*pi/19.0);}

        }

    }

}

if (threshold_division==0.0){ /* division following the maximal stress = particular case of the previous rule */

    division_highest_tension:; //printf("%f %d %d\n",stress[3*dividing_cell+2],dividing_cell,cell_number);

    if (stress[3*dividing_cell+2]!=0){

    eigenvalue=(stress[3*dividing_cell]+stress[3*dividing_cell+1])/2.;

    eigenvalue+=0.5*sqrt((stress[3*dividing_cell]-stress[3*dividing_cell+1])*(stress[3*dividing_cell]-stress[3*dividing_cell+1])+

                    4.*stress[3*dividing_cell+2]*stress[3*dividing_cell+2]);

    division_direction[0]=stress[3*dividing_cell+2];

    division_direction[1]=-stress[3*dividing_cell]+eigenvalue;

    norm=sqrt(division_direction[0]*division_direction[0]+division_direction[1]*division_direction[1]);

    division_direction[0]/=norm;

    division_direction[1]/=norm;}

    else if (stress[3*dividing_cell]>stress[3*dividing_cell+1]){

    division_direction[0]=1;

    division_direction[1]=0;}

    else if (stress[3*dividing_cell]<stress[3*dividing_cell+1]){

    division_direction[0]=0;

    division_direction[1]=1;}

    else if (stress[3*dividing_cell]==stress[3*dividing_cell+1]){

    division_direction[0]=rand()/(double)RAND_MAX;

    division_direction[1]=sqrt(1.-division_direction[0]*division_direction[0]);}

}

double temp_shape[3];

temp_shape[0]=shapes[3*dividing_cell];

temp_shape[1]=shapes[3*dividing_cell+1];

temp_shape[2]=shapes[3*dividing_cell+2];

//if ((barycenters[2*dividing_cell]*barycenters[2*dividing_cell])<(width*width)){ temp_shape[0]*=dilation*dilation; temp_shape[2]*=dilation;}

if (threshold_division<0.0){ /* geometrical rule: division following the shortest direction */

    double shape_direction[2];

    eigenvalue=(temp_shape[0]+temp_shape[1])/2.;

    eigenvalue+=0.5*sqrt((temp_shape[0]-temp_shape[1])*(temp_shape[0]-temp_shape[1])+

                    4.*temp_shape[2]*temp_shape[2]);

    shape_direction[0]=temp_shape[2];

    shape_direction[1]=eigenvalue-temp_shape[0];

    norm=sqrt(shape_direction[0]*shape_direction[0]+shape_direction[1]*shape_direction[1]);

    shape_direction[0]/=norm;

    shape_direction[1]/=norm;

    division_direction[0]=shape_direction[1];

    division_direction[1]=-shape_direction[0];

}

/* computation of the new vertices as the intersection between the walls of the cell and the division direction */

int l=0,m=0,max=vertices_number_in_each_cell[dividing_cell+1],min=vertices_number_in_each_cell[dividing_cell],previous;

double test;

for (int i=min;i<max;i++){ /* loop over the vertices of the dividing cell */

    previous=i-1; if (i==min){previous=max-1;}

    test=-division_direction[1]*(vertices[2*cells_vertices[previous]]-barycenters[2*dividing_cell]);

    test+=division_direction[0]*(vertices[2*cells_vertices[previous]+1]-barycenters[2*dividing_cell+1]);

    test/=-division_direction[1]*(vertices[2*cells_vertices[previous]]-vertices[2*cells_vertices[i]])+

            division_direction[0]*(vertices[2*cells_vertices[previous]+1]-vertices[2*cells_vertices[i]+1]);

        /* test if the considered point of the wall is on the division direction */

    if ((test<1. && test>=0.) && (m==0 || m==1) && l<2){

        //if (test>0.9){test=0.9;} if (test<0.1){test=0.1;} // if new vertices should not be too closed of the older ones for stability reasons

        l++;

        if (l==1){

            new_vertex1[0]=test*(vertices[2*cells_vertices[i]]-vertices[2*cells_vertices[previous]])+vertices[2*cells_vertices[previous]];

            new_vertex1[1]=test*(vertices[2*cells_vertices[i]+1]-vertices[2*cells_vertices[previous]+1])+vertices[2*cells_vertices[previous]+1];

            cell1[daughter_cells_vertices_number[0]]=vertices_number-2; daughter_cells_vertices_number[0]+=1;

            cell2[daughter_cells_vertices_number[1]]=vertices_number-2; daughter_cells_vertices_number[1]+=1;

            wall1[0]=cells_vertices[i];

            wall1[1]=cells_vertices[previous];}

        if (l==2){

            new_vertex2[0]=test*(vertices[2*cells_vertices[i]]-vertices[2*cells_vertices[previous]])+vertices[2*cells_vertices[previous]];

            new_vertex2[1]=test*(vertices[2*cells_vertices[i]+1]-vertices[2*cells_vertices[previous]+1])+vertices[2*cells_vertices[previous]+1];

            cell1[daughter_cells_vertices_number[0]]=vertices_number-1; daughter_cells_vertices_number[0]+=1;

            cell2[daughter_cells_vertices_number[1]]=vertices_number-1; daughter_cells_vertices_number[1]+=1;

            wall2[0]=cells_vertices[i];

            wall2[1]=cells_vertices[previous];}

            m=2;

        }

    else if ((test<1. && test>=0.) && (m==0 || m==2) && l<2){

        //if (test>0.9){test=0.9;} if (test<0.1){test=0.1;}

        l++;

        if (l==1){

            new_vertex1[0]=test*(vertices[2*cells_vertices[i]]-vertices[2*cells_vertices[previous]])+vertices[2*cells_vertices[previous]];

            new_vertex1[1]=test*(vertices[2*cells_vertices[i]+1]-vertices[2*cells_vertices[previous]+1])+vertices[2*cells_vertices[previous]+1];

            cell1[daughter_cells_vertices_number[0]]=vertices_number-2; daughter_cells_vertices_number[0]+=1;

            cell2[daughter_cells_vertices_number[1]]=vertices_number-2; daughter_cells_vertices_number[1]+=1;

            wall1[0]=cells_vertices[i];

            wall1[1]=cells_vertices[previous];}

        if (l==2){

            new_vertex2[0]=test*(vertices[2*cells_vertices[i]]-vertices[2*cells_vertices[previous]])+vertices[2*cells_vertices[previous]];

            new_vertex2[1]=test*(vertices[2*cells_vertices[i]+1]-vertices[2*cells_vertices[previous]+1])+vertices[2*cells_vertices[previous]+1];

            cell1[daughter_cells_vertices_number[0]]=vertices_number-1; daughter_cells_vertices_number[0]+=1;

            cell2[daughter_cells_vertices_number[1]]=vertices_number-1; daughter_cells_vertices_number[1]+=1;

            wall2[0]=cells_vertices[i];

            wall2[1]=cells_vertices[previous];}

            m=1;

        }

    /* write the old vertices in the two daughter cells */

if (l==0 || l==2){cell1[daughter_cells_vertices_number[0]]=cells_vertices[i]; daughter_cells_vertices_number[0]+=1;}

if (l==1){cell2[daughter_cells_vertices_number[1]]=cells_vertices[i]; daughter_cells_vertices_number[1]+=1;}

}

/* if l<2, there has been a problem in the computation of the new vertices */

if (l!=2){printf("\n"); printf("division de %d ok,recherche des voisines\n",dividing_cell); exit(EXIT_FAILURE);}

/* the new vertices are saved in the list "vertices" */

vertices[2*(vertices_number-2)]=new_vertex1[0]; vertices[2*(vertices_number-2)+1]=new_vertex1[1];

vertices[2*(vertices_number-1)]=new_vertex2[0]; vertices[2*(vertices_number-1)+1]=new_vertex2[1];

/* find the cells which are adjacent to the dividing one and have to be modified to include the new vertices */

adjacent_cell1[0]=-1; adjacent_cell2[0]=-1; adjacent_cells_number[0]=0;

for (int i=0;i<cell_number-1;i++){ if (i==dividing_cell){goto nextcell;}

    min=vertices_number_in_each_cell[i]; max=vertices_number_in_each_cell[i+1]; previous=max-min-1;

    for (j=0;j<max-min;j++){

        if ((cells_vertices[min+j]==wall1[1] && cells_vertices[min+previous]==wall1[0])

            ||(cells_vertices[min+j]==wall2[1] && cells_vertices[min+previous]==wall2[0])

            ){

            if (adjacent_cells_number[0]==0){adjacent_cell1[0]=i;}

            if (adjacent_cells_number[0]==1){adjacent_cell2[0]=i;}

            adjacent_cells_number[0]+=1; goto nextcell;

            }

        previous=j;}

    nextcell:;

}

//printf("%d voisines trouvee(s): %d %d\n",adjacent_cells_number[0],adjacent_cell1[0],adjacent_cell2[0]);

}