Laboratoire RDP » Biophysics Team » Cell Division Rule

#include "constants.h"

#include "growth_and_chemistry.h"

#include <gsl/gsl_errno.h>

#include <math.h>

#include <algorithm>

#include <gsl/gsl_matrix.h>

#include <gsl/gsl_blas.h>

#include <gsl/gsl_eigen.h>

#include <gsl/gsl_math.h>

#include <gsl/gsl_cblas.h>

/* this function computes the equations of chemistry and growth */

int growth_and_chemistry(double t, const double x[], double dx[], void *parameters){

struct parameters_list2 {int cell_number; int *cells_vertices; int *vertices_number_in_each_cell; double *shapes; double *lengths; double *cell_areas; double *barycenters;

                         int *cell_walls; int wall_number; int *boundary; int *wall_vertices; double *vertices; double *growth_rate; double *noise;};

/* diagonalize the deformation tensor in order to compute growth */

double delta_shapes[3*(((parameters_list2*)parameters)->cell_number)];

for (int i=0; i<(((parameters_list2*)parameters)->cell_number);i++){

    gsl_matrix *gsl_delta_shapes=gsl_matrix_calloc(2,2);

    gsl_matrix_set(gsl_delta_shapes,0,0,(((parameters_list2*)parameters)->shapes)[3*i]-x[3*i]); gsl_matrix_set(gsl_delta_shapes,1,1,(((parameters_list2*)parameters)->shapes)[3*i+1]-x[3*i+1]);

    gsl_matrix_set(gsl_delta_shapes,1,0,(((parameters_list2*)parameters)->shapes)[3*i+2]-x[3*i+2]); gsl_matrix_set(gsl_delta_shapes,0,1,(((parameters_list2*)parameters)->shapes)[3*i+2]-x[3*i+2]);

    gsl_eigen_symmv_workspace *workspace=gsl_eigen_symmv_alloc(2); gsl_vector *eigenvalues=gsl_vector_alloc(2); gsl_matrix *eigenvectors=gsl_matrix_alloc(2,2); gsl_eigen_symmv(gsl_delta_shapes,eigenvalues,eigenvectors,workspace);

    if (gsl_vector_get(eigenvalues,0)<0.){ gsl_vector_set(eigenvalues,0,0.);}

    if (gsl_vector_get(eigenvalues,1)<0.){ gsl_vector_set(eigenvalues,1,0.);}

    delta_shapes[3*i]=gsl_vector_get(eigenvalues,0)*gsl_matrix_get(eigenvectors,0,0)*gsl_matrix_get(eigenvectors,0,0)+gsl_vector_get(eigenvalues,1)*gsl_matrix_get(eigenvectors,0,1)*gsl_matrix_get(eigenvectors,0,1);

    delta_shapes[3*i+1]=gsl_vector_get(eigenvalues,0)*gsl_matrix_get(eigenvectors,1,0)*gsl_matrix_get(eigenvectors,1,0)+gsl_vector_get(eigenvalues,1)*gsl_matrix_get(eigenvectors,1,1)*gsl_matrix_get(eigenvectors,1,1);

    delta_shapes[3*i+2]=gsl_vector_get(eigenvalues,0)*gsl_matrix_get(eigenvectors,0,0)*gsl_matrix_get(eigenvectors,1,0)+gsl_vector_get(eigenvalues,1)*gsl_matrix_get(eigenvectors,0,1)*gsl_matrix_get(eigenvectors,1,1);

    gsl_vector_free(eigenvalues); gsl_matrix_free(eigenvectors); gsl_matrix_free(gsl_delta_shapes); gsl_eigen_symmv_free(workspace);

}

    /* computation of the growth rate:

     the growth rate of periclinal walls is directly given by the list "growth_rate", and here we add the noise

     the growth of anticlinal walls is the average of the growth of the two neighbouring periclinal walls */

double cell_growth[(((parameters_list2*)parameters)->cell_number)];

double wall_growth[(((parameters_list2*)parameters)->wall_number)];

for (int i=0;i<(((parameters_list2*)parameters)->wall_number);i++){

    wall_growth[i]=0.;

}

for (int i=0;i<(((parameters_list2*)parameters)->cell_number);i++){

    //cell_growth[i]=1.;

    /* with noise */

    cell_growth[i]=(((parameters_list2*)parameters)->growth_rate)[i]*(((parameters_list2*)parameters)->noise)[i];

    for (int j=(((parameters_list2*)parameters)->vertices_number_in_each_cell)[i];j<(((parameters_list2*)parameters)->vertices_number_in_each_cell)[i+1];j++){

        wall_growth[(((parameters_list2*)parameters)->cell_walls)[j]]+=0.5*cell_growth[i]*(((parameters_list2*)parameters)->noise)[i];

        if ((((parameters_list2*)parameters)->boundary)[j]==1){wall_growth[(((parameters_list2*)parameters)->cell_walls)[j]]*=2.;}

    }

}

    /* growth of the periclinal walls, proportionnal to the deformation, with possiblity to set negative growth, corresponding to compression, to zero

     this was not used in the end since it made no difference */

for (int i=0;i<(((parameters_list2*)parameters)->cell_number);i++){

    dx[3*i]=cell_growth[i]*(1./(tau1*(((parameters_list2*)parameters)->cell_number)))*

               //std::max(0.,(((parameters_list2*)parameters)->shapes)[3*i]-x[3*i]);

               //(((parameters_list2*)parameters)->shapes)[3*i]-x[3*i];

               delta_shapes[3*i];

    dx[3*i+1]=cell_growth[i]*(1./(tau1*(((parameters_list2*)parameters)->cell_number)))*

               // std::max(0.,(((parameters_list2*)parameters)->shapes)[3*i+1]-x[3*i+1]);

               // (((parameters_list2*)parameters)->shapes)[3*i+1]-x[3*i+1];

              delta_shapes[3*i+1];

    dx[3*i+2]=cell_growth[i]*(1./(tau1*(((parameters_list2*)parameters)->cell_number)))*

              //  std::max(0.,(((parameters_list2*)parameters)->shapes)[3*i+2]-x[3*i+2]);

             //   (((parameters_list2*)parameters)->shapes)[3*i+2]-x[3*i+2];

               delta_shapes[3*i+2];

}

    /* growth of the anticlinal walls */

for (int i=0;i<(((parameters_list2*)parameters)->wall_number);i++){

        dx[3*(((parameters_list2*)parameters)->cell_number)+i]=wall_growth[i]*(1./(tau1*(((parameters_list2*)parameters)->cell_number)))*

                std::max((((parameters_list2*)parameters)->lengths)[i]-x[3*(((parameters_list2*)parameters)->cell_number)+i],0.);

}

return GSL_SUCCESS;

}