/*************************************************************/ /* Search and Capture with KT evolving after NEB */ /* Authors: R.Paul (IACS,Kolkata, India & NPB UC Davis) */ /* A. Mogilner (NYU, USA); 2015 */ /*************************************************************/ #include #include #include #include //Global Memory Declaration, Random number /**********dyanamic memory allocation********/ double *vector ( int nrl, int nrh); double **matrix ( int nrl, int nrh, int ncl , int nch ); double ***tensor3 ( int nxl, int nxh, int nyl, int nyh, int nzl, int nzh ); int *ivector ( int nrl, int nrh ); int **imatrix ( int nrl, int nrh, int ncl, int nch ); int ***itensor3 ( int nxl, int nxh, int nyl, int nyh, int nzl, int nzh ); double *vector( int nl, int nh) { double *v; v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double)); if (!v) printf("allocation failure in vector()\n\n"); return v-nl; } double **matrix ( int nrl, int nrh, int ncl, int nch ) { int i; double **m; m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*)); if (!m) printf("allocation failure 1 in matrix()\n\n"); m -= nrl; for(i=nrl;i<=nrh;i++) { m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double)); if (!m[i]) printf("allocation failure 2 in matrix()\n\n"); m[i] -= ncl; } return m; } double ***tensor3 ( int nxl, int nxh, int nyl, int nyh, int nzl, int nzh ) { int i, j; double ***m; m=(double ***) malloc((unsigned) (nxh-nxl+1)*sizeof(double*)); if (!m) printf("allocation failure 1 in tensor3()\n\n"); m -= nxl; for(i=nxl;i<=nxh;i++) { m[i]=(double **) malloc((unsigned) (nyh-nyl+1)*sizeof(double*)); if (!m[i]) printf("allocation failure 2 in tensor3()\n\n"); m[i] -= nyl; }; for(i=nxl;i<=nxh;i++) { for(j=nyl;j<=nyh;j++) { m[i][j]=(double *) malloc((unsigned) (nzh-nzl+1)*sizeof(double)); if (!m[i][j]) printf("allocation failure 3 in tensor3()\n\n"); m[i][j] -= nzl; } }; return m; } int *ivector( int nl, int nh) { int *v; v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int)); if (!v) printf("allocation failure in ivector()\n\n"); return v-nl; } int **imatrix ( int nrl, int nrh, int ncl, int nch ) { int i; int **m; m=(int **) malloc((unsigned) (nrh-nrl+1)*sizeof(int*)); if (!m) printf("allocation failure 1 in matrix()\n\n"); m -= nrl; for(i=nrl;i<=nrh;i++) { m[i]=(int *) malloc((unsigned) (nch-ncl+1)*sizeof(int)); if (!m[i]) printf("allocation failure 2 in matrix()\n\n"); m[i] -= ncl; } return m; } int ***itensor3 ( int nxl, int nxh, int nyl, int nyh, int nzl, int nzh ) { int i, j; int ***m; m=(int ***) malloc((unsigned) (nxh-nxl+1)*sizeof(int*)); if (!m) printf("allocation failure 1 in tensor3()\n\n"); m -= nxl; for(i=nxl;i<=nxh;i++) { m[i]=(int **) malloc((unsigned) (nyh-nyl+1)*sizeof(int*)); if (!m[i]) printf("allocation failure 2 in tensor3()\n\n"); m[i] -= nyl; }; for(i=nxl;i<=nxh;i++) { for(j=nyl;j<=nyh;j++) { m[i][j]=(int *) malloc((unsigned) (nzh-nzl+1)*sizeof(int)); if (!m[i][j]) printf("allocation failure 3 in tensor3()\n\n"); m[i][j] -= nzl; } }; return m; } /************random number and varialble declaration*******/ #define N_BITS 31 #define IRND1 250 #define IRND2 103 #define MAX_RN 0X7FFFFFFF #define MAX_R2 0X40000000 int ir[256], irndx1[256], irndx2[256]; int Pw(int n); double ran2(int *idum); int Rand(); void rndini(); void Place_Chromosome_Kinetochore(); double Rot_matrix_ch(int i, int j, double phi, double theta, double psi, double x0, double y0, double z0); double Rot_matrix_kt(int i, int j, double phi, double theta, double psi, double x0, double y0, double z0); void Generate_MT(int t); void Generate_MT_DirCos(int i, int j); void Grow_MT(int t, int n_ch); void Shrink_MT(int t, int n_ch); int Check_MT_Capture(int Ni, int ve,double x, double y, double z, double xold, double yold, double zold, int t); int Check_MT_hit_CH(int Ni, int ve, double x, double y, double z, double xold, double yold, double zold, int t); int Check_all_kt_capture(int t); double Shortest_Distance(double x1, double y1, double z1, double x2, double y2, double z2, double x, double y, double z); double Dot_Product(double x1,double y1,double z1, double x2,double y2,double z2, double x3,double y3,double z3, double x4,double y4,double z4); void Visualization_povray(int t, int n_ch); void Output_ch_kt_coordinate(); void Check_length_MT(); void Output_cat_rescue(int ensem, int eff_ensem); void Output_cap_time(int n_ch, int eff_ensem); double DistP1P2(double x1, double y1, double z1, double x2, double y2, double z2); void Rotate_after_first_attach(int t); void Rotate(int i, double psi, double theta, double phi); void Translate(int i, double x0, double y0, double z0); void Trans_matrix_ch(int i, int j, double x0, double y0, double z0); void Trans_matrix_kt(int i, int j, double x0, double y0, double z0); int Check_if_outside(double x, double y, double z); int IntoTheSphere(double x,double y,double z, double r); void RotateSingleChrom(int t, int i); void MoveCentrosome(int t, int f_v); int OutOfSpheroid(double x,double y,double z,double a, double b, double c); void DisolveSyntMero(int t); void Disolve_Synt(int i); void Disolve_Mero(int i, int x, int y); void Disolve_Mono(int i, int x, int y); void DetachBondswithAging(int t); void Output_Errors(int n_ch, int ensem); int Check_MT_Capture_Monotelic(int Ni, int ve,double x, double y, double z, double xold, double yold, double zold, int t); void MoveChromToCentro(int t, int i, int cnt); void MoveMonotelic(int t); void ChangeKTShapeCylinder(int t); void ChangeKTShapeCrescent(int t); double ShiftOfCrescentKT(); void LaterallyAttachChromosome(int t); void FindLateralAngleKT(int ensem, int input_flag); int Lx,Ly, Lz,N, counter = 0, seed = 12357, idum = 12357,size,tag, MC_max, MAX_ensem=1, N_ch,N_kt,t_max, t_max_init, Nmt_cap_max, Nmt_cap_l,Nmt_cap_r,Nmt_cap_u,Nmt_cap_d, MAX_KT_t=1, NMT_max, NCentro_max, Nch_for_hist, MAX_disolv_t, t_NEB, T_separate = 50*60,Delta_t, Tdelay_mean, Tdelay, Tmono, segT,N_ch_min, matureKT=0, label_seg=0, T_evolve,Rot_flag = 0; int *nn, *nnn, *node, *state, *label_kt_l, *label_kt_r, **labelMT, *NMT, cap_time,N_rescue, *Nkt_cap_l, *Nkt_cap_r,*T_cap, Nkt_cap_max, N_ch_max, t_count_l, t_count_r, t_count_u, t_count_d, *num_succ, *num_fail, **KT_l_pole, **KT_r_pole, *Err_KT_l, **KT_l_time,**KT_r_time, *Err_KT_r,*Err_mono, *Err_noatch, *Err_synt, *Err_mero, *Err_none, *KT_off_l, *KT_off_r, *KT_t_l, *KT_t_r, *Lktcapby,*Rktcapby, *rotlab_l,*rotlab_r, *NMTcap, *SgnCent, *KT_t_syn, *KT_t_mero, *KT_t_mono, *Count_ErNone,*Count_ErMero,*Count_ErSyn,*Count_ErMon,*Count_ErNoat,*Capture_KT_NUM, *label_kt_l_cir, *label_kt_r_cir; double pi=3.1415926535897932, l_ch=1.0,l_kt=l_ch/4.0, R_ch=l_ch/10.0, R_kt=R_ch, l_tub=R_kt, l_gr, l_sr,d_cent=10.0, R_nuc = d_cent/1.5, X_C, Y_C, Z_C, P_grow = 0.8, P_rescue, X_var=d_cent/3., Y_var=d_cent/3., v_g, v_s,l_av_mt_l,l_av_mt_r, l_av_mt_u, l_av_mt_d, f_move = 0.5,Ax,Ay,Az, D_complete,R_cent, Dcent_ini, l_gr_ini,ax,ay,az, koff_mono, koff_synt, koff_amphi, koff_mero, V_centro,Rnuc_eff, d_overlap, rnd_shift, f_cat, f_res, P_cat, P_res, Tau, Run_time, R_kt_cir,l_kt_cir,d_kt_shift, Ax_out, Ax_in, Az_in, Az_out, Ay_in, Ay_out, hkt, wkt,MP_width_x=2.0,MP_width_y=2.0,MP_width_z=2.0,arc_l,arc_l_fin,psi_lateral, R_kt_cir_fin,d_kt_shift_fin; double *PCHu_x,*PCHu_y,*PCHu_z,*PCHd_x,*PCHd_y,*PCHd_z,*PCHm_x,*PCHm_y,*PCHm_z, *PKTl_x,*PKTl_y,*PKTl_z, *PKTr_x,*PKTr_y,*PKTr_z, **PMT_x, **PMT_y,**PMT_z, *PKTl0_x,*PKTl0_y,*PKTl0_z, *PKTr0_x,*PKTr0_y,*PKTr0_z, *PKTuL_x,*PKTuL_y,*PKTuL_z,*PKTdL_x,*PKTdL_y,*PKTdL_z, *PKTuR_x,*PKTuR_y,*PKTuR_z,*PKTdR_x,*PKTdR_y,*PKTdR_z, *PKTmL_x,*PKTmL_y,*PKTmL_z,*PKTmR_x,*PKTmR_y,*PKTmR_z, **Alpha, **Beta, **Gamma,**c_Alpha, **c_Beta, **c_Gamma, **Centro_x,**Centro_y,**Centro_z,*x_tr, *y_tr, *z_tr, *N_cat, *x_cap_l, *y_cap_l, *z_cap_l,*x_cap_r, *y_cap_r, *z_cap_r, *L_av, **PMTcap_x, **PMTcap_y, **PMTcap_z, *psi_n,*theta_n,*phi_n; /************ Calculation of 2 ***********************/ int Pw(int n) { return((int)(pow(2.0,(double)n) + 0.1)); } void nrerror0(const char error_text[]) { printf("Numerical Recipes run-time error...\n"); printf("%s\n",error_text); printf("...now exiting to system...\n"); exit(1); } /***** random number generator ran2 ***********************/ #include #define M 714025 #define IA 1366 #define IC 150889 double ran2(int *idum) { static long iy,irn[98]; static int iff=0; int j; void nrerror(); if (*idum < 0 || iff == 0) { iff=1; if ((*idum=(IC-(*idum)) % M) < 0) *idum = -(*idum); for (j=1;j<=97;j++) { *idum=(IA*(*idum)+IC) % M; irn[j]=(*idum); } *idum=(IA*(*idum)+IC) % M; iy=(*idum); } j=1 +(int)( 97.0*iy/M); if (j > 97 || j < 1) nrerror0("RAN2: This cannot happen."); iy=irn[j]; *idum=(IA*(*idum)+IC) % M; irn[j]=(*idum); return (double)iy/M; } #undef M #undef IA #undef IC /****************************************************************/ /*************** Sub Routines ************************/ double norm(double x, double y, double z){ return(sqrt(x*x+y*y+z*z)); } double dotprod(double x1, double y1, double z1, double x2, double y2, double z2){ return(x1*x2+y1*y2+z1*z2); } int Check_if_outside(double x, double y, double z){ double lngth = sqrt(x*x+y*y+z*z); if(lngth0.0){ return(d_kt_shift); }else{ printf("d_kt_shift is out of range....\n"); exit(-1); }; } //============================= void Place_Chromosome_Kinetochore(int ensem, int input_flag){ int i,j; double x0,y0,z0, d=5.0*R_nuc, phi,theta,psi; double x,y,z; char fn[70]; FILE *fp; /*origin (0,0,0) = center of chromosome (ch) */ for(i=1;i<=N_ch;i++){ do{ int flag_within_vol = 0; do{ x=-MP_width_x*Ax_out; y=MP_width_y*Ay_out; //metaphase plate width MP_width (0-2) z=MP_width_z*Az_out; int pointisin=OutOfSpheroid(x,y,z,Ax_out,Ay_out,Az_out); int pointisout=OutOfSpheroid(x,y,z,Ax_in,Ay_in,Az_in); if((pointisin==0)&&(pointisout==1))flag_within_vol = 1; }while(flag_within_vol==0); //center of the chromosome x_tr[i]=x; y_tr[i]=y; z_tr[i]=z; //distribute the CH/KT such that they have minimal overlap d=5.0*R_nuc;//large reference distance for(j=1;j<=i-1;j++){ double d1 = norm((x_tr[i]-x_tr[j]),(y_tr[i]-y_tr[j]),(z_tr[i]-z_tr[j])); if(d1dist_min)break; }; psi -= pi/2.0; //For this value of psi lateral attachment is achieved psi_lateral = psi; //psi=theta=phi=0 is chosen such that KT is perpendicular to pole-pole axis psi = -pi/2.0;theta = 0.0; phi = 0.0; //initial position of CH before translation PCHu_x[i] = 0.;PCHu_y[i] = 0.;PCHu_z[i] = l_ch/2; PCHd_x[i] = 0.;PCHd_y[i] = 0.;PCHd_z[i] = -l_ch/2; PCHm_x[i] = 0.;PCHm_y[i] = 0.;PCHm_z[i] = 0.; //kinetochores (kt) are perpendicular to the ch axis. //the left kt extends from PKTl0_x --> PKTl_x //the right kt extends from PKTr0_x --> PKTr_x*/ //Always keep PKTl0 coincide with the center of chromosome //this is important to determine the size of the ring KT PKTl_x[i]=0.; PKTl_y[i]=-l_kt-R_ch; PKTl_z[i]=0.; PKTl0_x[i]=0.; PKTl0_y[i]=0.0; PKTl0_z[i]=0.; PKTr_x[i]=0.; PKTr_y[i]=l_kt+R_ch; PKTr_z[i]=0.; PKTr0_x[i]=0.; PKTr0_y[i]=0.0; PKTr0_z[i]=0.; //Find the center of the left and right ring KT double dist_kt_tip = norm(PKTl_x[i]-PKTl0_x[i],PKTl_y[i]-PKTl0_y[i],PKTl_z[i]-PKTl0_z[i]); PKTmL_x[i] = PCHm_x[i] + d_kt_shift*(PKTl_x[i]-PKTl0_x[i])/dist_kt_tip; PKTmL_y[i] = PCHm_y[i] + d_kt_shift*(PKTl_y[i]-PKTl0_y[i])/dist_kt_tip; PKTmL_z[i] = PCHm_z[i] + d_kt_shift*(PKTl_z[i]-PKTl0_z[i])/dist_kt_tip; dist_kt_tip = norm(PKTr_x[i]-PKTr0_x[i],PKTr_y[i]-PKTr0_y[i],PKTr_z[i]-PKTr0_z[i]); PKTmR_x[i] = PCHm_x[i] + d_kt_shift*(PKTr_x[i]-PKTr0_x[i])/dist_kt_tip; PKTmR_y[i] = PCHm_y[i] + d_kt_shift*(PKTr_y[i]-PKTr0_y[i])/dist_kt_tip; PKTmR_z[i] = PCHm_z[i] + d_kt_shift*(PKTr_z[i]-PKTr0_z[i])/dist_kt_tip; PKTuL_x[i]=PKTmL_x[i]; PKTuL_y[i]=PKTmL_y[i]; PKTuL_z[i]=PKTmL_z[i]+l_kt_cir/2.0; PKTdL_x[i]=PKTmL_x[i]; PKTdL_y[i]=PKTmL_y[i]; PKTdL_z[i]=PKTmL_z[i]-l_kt_cir/2.0; PKTuR_x[i]=PKTmR_x[i]; PKTuR_y[i]=PKTmR_y[i]; PKTuR_z[i]=PKTmR_z[i]+l_kt_cir/2.0; PKTdR_x[i]=PKTmR_x[i]; PKTdR_y[i]=PKTmR_y[i]; PKTdR_z[i]=PKTmR_z[i]-l_kt_cir/2.0; /*********** Initialization ***********/ //initial label of the KT label_kt_l[i] = 0; label_kt_r[i] = 0; label_kt_l_cir[i]=0; label_kt_r_cir[i]=0; //presegregated KT if(label_seg==1){ label_kt_l_cir[i]=1; label_kt_r_cir[i]=1; }; //Number of captured MT by KT Nkt_cap_l[i] = Nkt_cap_r[i] = 0; //From which pole the KT is captured for(j=0;j<=Nkt_cap_max;j++){ KT_l_pole[i][j] = KT_r_pole[i][j] = 0; KT_l_time[i][j] = KT_r_time[i][j] = 0; }; //switch on attachment KT_off_l[i] = KT_off_r[i] = 0; // attachment errors is zero Err_KT_l[i] = Err_KT_r[i] = 0 ; //time after capture is zero KT_t_l[i] = KT_t_r[i] = 0; //duration of wrong attachment KT_t_mero[i] = KT_t_syn[i] = KT_t_mono[i] = 0; //Lkt must be captured from left pole //Rkt must be captured from right pole //but initially they can be captured by any pole Lktcapby[i]=Rktcapby[i]=0; rotlab_l[i] = rotlab_r[i] = 0; psi_n[i] = psi; theta_n[i] = theta; phi_n[i] = phi; Rot_matrix_ch(i,1,phi,theta,psi,x0,y0,z0); Rot_matrix_ch(i,-1,phi,theta,psi,x0,y0,z0); Rot_matrix_ch(i,0,phi,theta,psi,x0,y0,z0); Rot_matrix_kt(i,-1,phi,theta,psi,x0,y0,z0); Rot_matrix_kt(i,1,phi,theta,psi,x0,y0,z0); Rot_matrix_kt(i,-2,phi,theta,psi,x0,y0,z0); Rot_matrix_kt(i,2,phi,theta,psi,x0,y0,z0); }; } //============================= void FindLateralAngleKT(int ensem, int input_flag){ int i,j; double x0,y0,z0, d=5.0*R_nuc, phi,theta,psi; double x,y,z; char fn[70]; FILE *fp; /*origin (0,0,0) = center of chromosome (ch). place the ch along the z axis with two ends at (0,0,l/2) and (0,0,-l/2)*/ for(i=1;i<=N_ch;i++){ do{ int flag_within_vol = 0; do{ x=-MP_width_x*Ax_out; y=MP_width_y*Ay_out; //metaphase plate width MP_width (0-2) z=MP_width_z*Az_out; int pointisin=OutOfSpheroid(x,y,z,Ax_out,Ay_out,Az_out); int pointisout=OutOfSpheroid(x,y,z,Ax_in,Ay_in,Az_in); if((pointisin==0)&&(pointisout==1))flag_within_vol = 1; }while(flag_within_vol==0); //center of the chromosome x_tr[i]=x; y_tr[i]=y; z_tr[i]=z; //distribute the CH/KT such that they have minimal overlap d=5.0*R_nuc;//large reference distance for(j=1;j<=i-1;j++){ double d1 = norm((x_tr[i]-x_tr[j]),(y_tr[i]-y_tr[j]),(z_tr[i]-z_tr[j])); if(d1dist_min)break; }; psi -= pi/2.0; //For this value of psi lateral attachment is achieved psi_lateral = psi; }; } void LaterallyAttachChromosome(){ int i,j; double psi, theta, phi,x0,y0,z0; i = 1; psi_n[i] = psi_lateral; psi = psi_n[i]; theta=0.0; phi=0.0; x0 = x_tr[i]; y0 = y_tr[i]; z0 = z_tr[i]; //Check if the captured KT rotates in the direction of the pole it is captured from //initial position of CH before translation PCHu_x[i] = 0.;PCHu_y[i] = 0.;PCHu_z[i] = l_ch/2; PCHd_x[i] = 0.;PCHd_y[i] = 0.;PCHd_z[i] = -l_ch/2; PCHm_x[i] = 0.;PCHm_y[i] = 0.;PCHm_z[i] = 0.; PKTl_x[i]=0.; PKTl_y[i]=-l_kt-R_ch; PKTl_z[i]=0.; PKTl0_x[i]=0.; PKTl0_y[i]=0.0; PKTl0_z[i]=0.; PKTr_x[i]=0.; PKTr_y[i]=l_kt+R_ch; PKTr_z[i]=0.; PKTr0_x[i]=0.; PKTr0_y[i]=0.0; PKTr0_z[i]=0.; //Find the center of the left and right ring KT double dist_kt_tip = norm(PKTl_x[i]-PKTl0_x[i],PKTl_y[i]-PKTl0_y[i],PKTl_z[i]-PKTl0_z[i]); PKTmL_x[i] = PCHm_x[i] + d_kt_shift*(PKTl_x[i]-PKTl0_x[i])/dist_kt_tip; PKTmL_y[i] = PCHm_y[i] + d_kt_shift*(PKTl_y[i]-PKTl0_y[i])/dist_kt_tip; PKTmL_z[i] = PCHm_z[i] + d_kt_shift*(PKTl_z[i]-PKTl0_z[i])/dist_kt_tip; dist_kt_tip = norm(PKTr_x[i]-PKTr0_x[i],PKTr_y[i]-PKTr0_y[i],PKTr_z[i]-PKTr0_z[i]); PKTmR_x[i] = PCHm_x[i] + d_kt_shift*(PKTr_x[i]-PKTr0_x[i])/dist_kt_tip; PKTmR_y[i] = PCHm_y[i] + d_kt_shift*(PKTr_y[i]-PKTr0_y[i])/dist_kt_tip; PKTmR_z[i] = PCHm_z[i] + d_kt_shift*(PKTr_z[i]-PKTr0_z[i])/dist_kt_tip; PKTuL_x[i]=PKTmL_x[i]; PKTuL_y[i]=PKTmL_y[i]; PKTuL_z[i]=PKTmL_z[i]+l_kt_cir/2.0; PKTdL_x[i]=PKTmL_x[i]; PKTdL_y[i]=PKTmL_y[i]; PKTdL_z[i]=PKTmL_z[i]-l_kt_cir/2.0; PKTuR_x[i]=PKTmR_x[i]; PKTuR_y[i]=PKTmR_y[i]; PKTuR_z[i]=PKTmR_z[i]+l_kt_cir/2.0; PKTdR_x[i]=PKTmR_x[i]; PKTdR_y[i]=PKTmR_y[i]; PKTdR_z[i]=PKTmR_z[i]-l_kt_cir/2.0; //exit(-1); Rot_matrix_ch(i,1,phi,theta,psi,x0,y0,z0); Rot_matrix_ch(i,-1,phi,theta,psi,x0,y0,z0); Rot_matrix_ch(i,0,phi,theta,psi,x0,y0,z0); Rot_matrix_kt(i,-1,phi,theta,psi,x0,y0,z0); Rot_matrix_kt(i,1,phi,theta,psi,x0,y0,z0); Rot_matrix_kt(i,-2,phi,theta,psi,x0,y0,z0); Rot_matrix_kt(i,2,phi,theta,psi,x0,y0,z0); Translate(i,x0,y0,z0); } double Rot_matrix_ch(int i,int j, double phi, double theta, double psi, double x0, double y0, double z0){ double x,y,z,x1,y1,z1; /* i = choromosome number; j = 1, -1, 0 labels for up down and middle coordinate x,y,z = coordinate of the before rotation; x1,y1,z1 = coordinates after rotation; (x0,y0,z0) = translation vector. phi = 1st rotation angle about z axis; theta = 2nd rotation angle about new x=x' axis; psi = 3rd rotation angle about new z=z' axis; all rotations are anti-clockwise. */ if(j==1){ x=PCHu_x[i]; y=PCHu_y[i]; z=PCHu_z[i]; x1=(cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi))*x+ (cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi))*y+ (sin(psi)*sin(theta))*z; y1=(-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi))*x+ (-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi))*y+ (cos(psi)*sin(theta))*z; z1=(sin(theta)*sin(phi))*x+ (-sin(theta)*cos(phi))*y+ (cos(theta))*z; //after rotation and translation PCHu_x[i]=x1+x0;PCHu_y[i]=y1+y0;PCHu_z[i]=z1+z0; }; if(j==-1){ x=PCHd_x[i]; y=PCHd_y[i]; z=PCHd_z[i]; x1=(cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi))*x+ (cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi))*y+ (sin(psi)*sin(theta))*z; y1=(-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi))*x+ (-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi))*y+ (cos(psi)*sin(theta))*z; z1=(sin(theta)*sin(phi))*x+ (-sin(theta)*cos(phi))*y+ (cos(theta))*z; //after rotation and translation PCHd_x[i]=x1+x0;PCHd_y[i]=y1+y0;PCHd_z[i]=z1+z0; }; if(j==0){ x=PCHm_x[i]; y=PCHm_y[i]; z=PCHm_z[i]; x1=(cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi))*x+ (cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi))*y+ (sin(psi)*sin(theta))*z; y1=(-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi))*x+ (-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi))*y+ (cos(psi)*sin(theta))*z; z1=(sin(theta)*sin(phi))*x+ (-sin(theta)*cos(phi))*y+ (cos(theta))*z; //after rotation and translation PCHm_x[i]=x1+x0;PCHm_y[i]=y1+y0;PCHm_z[i]=z1+z0; }; } double Rot_matrix_kt(int i, int j,double phi, double theta, double psi, double x0, double y0, double z0){ double x,y,z,x1,y1,z1; /* i = choromosome number; j = 1, -1 labels for left and right coordinate x,y,z = coordinate before rotation; x1,y1,z1 = coordinates after rotation; (x0,y0,z0) = translation vector. phi = 1st rotation angle about z axis; theta = 2nd rotation angle about new x=x' axis; psi = 3rd rotation angle about new z=z' axis; all rotations are anti-clockwise. */ if(j==-1){ x=PKTl_x[i];y=PKTl_y[i];z=PKTl_z[i]; x1=(cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi))*x+ (cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi))*y+ (sin(psi)*sin(theta))*z; y1=(-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi))*x+ (-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi))*y+ (cos(psi)*sin(theta))*z; z1=(sin(theta)*sin(phi))*x+ (-sin(theta)*cos(phi))*y+ (cos(theta))*z; PKTl_x[i]=x1+x0;PKTl_y[i]=y1+y0;PKTl_z[i]=z1+z0; x=PKTl0_x[i];y=PKTl0_y[i];z=PKTl0_z[i]; x1=(cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi))*x+ (cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi))*y+ (sin(psi)*sin(theta))*z; y1=(-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi))*x+ (-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi))*y+ (cos(psi)*sin(theta))*z; z1=(sin(theta)*sin(phi))*x+ (-sin(theta)*cos(phi))*y+ (cos(theta))*z; PKTl0_x[i]=x1+x0;PKTl0_y[i]=y1+y0;PKTl0_z[i]=z1+z0; }; if(j==1){ x=PKTr_x[i];y=PKTr_y[i];z=PKTr_z[i]; x1=(cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi))*x+ (cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi))*y+ (sin(psi)*sin(theta))*z; y1=(-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi))*x+ (-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi))*y+ (cos(psi)*sin(theta))*z; z1=(sin(theta)*sin(phi))*x+ (-sin(theta)*cos(phi))*y+ (cos(theta))*z; PKTr_x[i]=x1+x0;PKTr_y[i]=y1+y0;PKTr_z[i]=z1+z0; x=PKTr0_x[i];y=PKTr0_y[i];z=PKTr0_z[i]; x1=(cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi))*x+ (cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi))*y+ (sin(psi)*sin(theta))*z; y1=(-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi))*x+ (-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi))*y+ (cos(psi)*sin(theta))*z; z1=(sin(theta)*sin(phi))*x+ (-sin(theta)*cos(phi))*y+ (cos(theta))*z; PKTr0_x[i]=x1+x0;PKTr0_y[i]=y1+y0;PKTr0_z[i]=z1+z0; }; if(j==-2){ x=PKTuL_x[i];y=PKTuL_y[i];z=PKTuL_z[i]; x1=(cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi))*x+ (cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi))*y+ (sin(psi)*sin(theta))*z; y1=(-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi))*x+ (-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi))*y+ (cos(psi)*sin(theta))*z; z1=(sin(theta)*sin(phi))*x+ (-sin(theta)*cos(phi))*y+ (cos(theta))*z; PKTuL_x[i]=x1+x0;PKTuL_y[i]=y1+y0;PKTuL_z[i]=z1+z0; x=PKTdL_x[i];y=PKTdL_y[i];z=PKTdL_z[i]; x1=(cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi))*x+ (cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi))*y+ (sin(psi)*sin(theta))*z; y1=(-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi))*x+ (-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi))*y+ (cos(psi)*sin(theta))*z; z1=(sin(theta)*sin(phi))*x+ (-sin(theta)*cos(phi))*y+ (cos(theta))*z; PKTdL_x[i]=x1+x0;PKTdL_y[i]=y1+y0;PKTdL_z[i]=z1+z0; x=PKTmL_x[i];y=PKTmL_y[i];z=PKTmL_z[i]; x1=(cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi))*x+ (cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi))*y+ (sin(psi)*sin(theta))*z; y1=(-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi))*x+ (-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi))*y+ (cos(psi)*sin(theta))*z; z1=(sin(theta)*sin(phi))*x+ (-sin(theta)*cos(phi))*y+ (cos(theta))*z; PKTmL_x[i]=x1+x0;PKTmL_y[i]=y1+y0;PKTmL_z[i]=z1+z0; }; if(j==2){ x=PKTuR_x[i];y=PKTuR_y[i];z=PKTuR_z[i]; x1=(cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi))*x+ (cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi))*y+ (sin(psi)*sin(theta))*z; y1=(-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi))*x+ (-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi))*y+ (cos(psi)*sin(theta))*z; z1=(sin(theta)*sin(phi))*x+ (-sin(theta)*cos(phi))*y+ (cos(theta))*z; PKTuR_x[i]=x1+x0;PKTuR_y[i]=y1+y0;PKTuR_z[i]=z1+z0; x=PKTdR_x[i];y=PKTdR_y[i];z=PKTdR_z[i]; x1=(cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi))*x+ (cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi))*y+ (sin(psi)*sin(theta))*z; y1=(-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi))*x+ (-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi))*y+ (cos(psi)*sin(theta))*z; z1=(sin(theta)*sin(phi))*x+ (-sin(theta)*cos(phi))*y+ (cos(theta))*z; PKTdR_x[i]=x1+x0;PKTdR_y[i]=y1+y0;PKTdR_z[i]=z1+z0; x=PKTmR_x[i];y=PKTmR_y[i];z=PKTmR_z[i]; x1=(cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi))*x+ (cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi))*y+ (sin(psi)*sin(theta))*z; y1=(-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi))*x+ (-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi))*y+ (cos(psi)*sin(theta))*z; z1=(sin(theta)*sin(phi))*x+ (-sin(theta)*cos(phi))*y+ (cos(theta))*z; PKTmR_x[i]=x1+x0;PKTmR_y[i]=y1+y0;PKTmR_z[i]=z1+z0; }; } double Dot_Product(double x1,double y1,double z1, double x2,double y2,double z2, double x3,double y3,double z3, double x4,double y4,double z4){ //P1=(x1,y1,z1); P2=(x2,y2,z2); P3=(x3,y3,z3); P4=(x4,y4,z4); //find dot product of P1P2(.)P3P4 double P1P2x = x2-x1; double P1P2y = y2-y1; double P1P2z = z2-z1; double P3P4x = x4-x3; double P3P4y = y4-y3; double P3P4z = z4-z3; double P1P2_P3P4 = P1P2x*P3P4x+P1P2y*P3P4y+P1P2z*P3P4z; return(P1P2_P3P4); } void Generate_MT(int t){ int i, flag, count,ncnt,ncnt_check; double a1,a2,a3,r; double x,y,z; //starting point of MT are centrosomes. I.e., //place centrosomes (complete separation) Centro_x[1][t] = 0.0; Centro_y[1][t] = -d_cent; Centro_z[1][t] = 0.0; SgnCent[1] = 1; Centro_x[2][t] = 0.0; Centro_y[2][t] = d_cent; Centro_z[2][t] = 0.0 ; SgnCent[2] = 2; //initialize labels, growth=1, shrink=-1, captured=0 for(i=1;i<=NMT[1];i++)labelMT[1][i] = 0; for(i=1;i<=NMT[2];i++)labelMT[2][i] = 0; for(ncnt=1;ncnt<=NCentro_max;ncnt++){ for(i=1;i<=NMT[ncnt];i++){ do{ //generate random unit vector for the direction of MT a1 = (0.5-ran2(&idum)); a2 = (0.5-ran2(&idum)); a3 = (0.5-ran2(&idum)); r = norm(a1,a2,a3); a1 /=r; a2 /=r; a3 /=r; labelMT[ncnt][i]=1; //label_l = 1 (grow) PMT_x[ncnt][i] = l_gr_ini*a1+Centro_x[ncnt][t]; PMT_y[ncnt][i] = l_gr_ini*a2+Centro_y[ncnt][t]; PMT_z[ncnt][i] = l_gr_ini*a3+Centro_z[ncnt][t]; }while(OutOfSpheroidMT(PMT_x[ncnt][i],PMT_y[ncnt][i],PMT_z[ncnt][i])==1); }; }; } void Generate_MT_DirCos(int ncnt, int i){ double x,y,z,r; //ncnt is the number as well as the sign of the pole x = 0.5-ran2(&idum); y = 0.5-ran2(&idum); z = 0.5-ran2(&idum); r = sqrt(x*x+y*y+z*z); c_Alpha[ncnt][i] = x/r; c_Beta[ncnt][i] = y/r; c_Gamma[ncnt][i] = z/r; } void Generate_Single_MT(int ncnt, int i, int t){ double a1,a2,a3,r; int ncnt_check; do{ //generate random unit vector for the direction of MT a1 = (0.5-ran2(&idum)); a2 = (0.5-ran2(&idum)); a3 = (0.5-ran2(&idum)); r = norm(a1,a2,a3); a1 /=r; a2 /=r; a3 /=r; labelMT[ncnt][i]=1; //label_l = 1 (grow) PMT_x[ncnt][i] = l_gr_ini*a1+Centro_x[ncnt][t]; PMT_y[ncnt][i] = l_gr_ini*a2+Centro_y[ncnt][t]; PMT_z[ncnt][i] = l_gr_ini*a3+Centro_z[ncnt][t]; }while(OutOfSpheroidMT(PMT_x[ncnt][i],PMT_y[ncnt][i],PMT_z[ncnt][i])==1); } void Grow_MT(int t, int n_ch){ int i,j, ncnt, flag_kt, flag_ch, flag_outside, flag_otherCentro, flag_catastrophe, count=0; double length_mt,px,py,pz,p0x,p0y,p0z,a1,a2,a3,r; double l_gr_mic; /*******************/ for(ncnt=1;ncnt<=NCentro_max;ncnt++){ count = 0; //counts how many MT is attached to kt for(i=1;i<=NMT[ncnt];i++){ if(labelMT[ncnt][i]==1){ double rnd = ran2(&idum); if(rnd>P_cat){ //then grow //generate unit vector in the direction of the MT //source of the MT p0x = Centro_x[ncnt][t]; p0y = Centro_y[ncnt][t]; p0z = Centro_z[ncnt][t]; //tip of the MT px = PMT_x[ncnt][i]; //x coordinate before growth py = PMT_y[ncnt][i]; //y coordinate before growth pz = PMT_z[ncnt][i]; //z coordinate before growth a1 = px-p0x; a2 = py-p0y; a3 = pz-p0z; r = norm(a1,a2,a3); //unit vector in the direction of the MT growth a1 /= r; a2 /= r; a3 /= r; //Grow MT in smaller steps //l_gr is the growth velocity per sec. //l_gr_mic = microsocipic growth (could be even at monomer level) l_gr_mic = l_gr/50.0; for(j=1;j<=50;j++){ if(labelMT[ncnt][i]==1){ PMT_x[ncnt][i] += l_gr_mic*a1; //x coordinate after growth PMT_y[ncnt][i] += l_gr_mic*a2; //y coordinate after growth PMT_z[ncnt][i] += l_gr_mic*a3; //z coordinate after growth length_mt = norm((PMT_x[ncnt][i]-Centro_x[ncnt][t]), (PMT_y[ncnt][i]-Centro_y[ncnt][t]), (PMT_z[ncnt][i]-Centro_z[ncnt][t])); //Check if MT tip is within the KT cylinder, -1= MTs from left pole flag_kt = Check_MT_Capture(i,ncnt,PMT_x[ncnt][i],PMT_y[ncnt][i],PMT_z[ncnt][i],px,py,pz,t); if(flag_kt == 1){//MT is captured by KT //store the location of the capt. MTs. Generate_Single_MT(ncnt,i,t); count++; }; if(flag_kt==-1){//sufficient MT is captured, no place for new KT: shrink MT labelMT[ncnt][i]=-1; break; }; //check if mt end has hit Ch then catastrophes flag_ch = Check_MT_hit_CH(i,ncnt,PMT_x[ncnt][i],PMT_y[ncnt][i],PMT_z[ncnt][i],px,py,pz,t); if(flag_ch == 1){//MT hit CH labelMT[ncnt][i]=-1; break; }; //Check if mt tip is out of nuclear sphere flag_outside = OutOfSpheroidMT(PMT_x[ncnt][i],PMT_y[ncnt][i],PMT_z[ncnt][i]); if(flag_outside==1){ labelMT[ncnt][i]=-1; //shrink otherwise break; }; }; }; }; }; }; }; } void Shrink_MT(int t, int n_ch, int ensem){ int i, ncnt; double rnd, length_mt,px,py,pz,p0x,p0y,p0z,a1,a2,a3,r; for(ncnt=1;ncnt<=NCentro_max;ncnt++){ /*******************/ for(i=1;i<=NMT[ncnt];i++){ if(labelMT[ncnt][i]==-1){ length_mt = norm((PMT_x[ncnt][i]-Centro_x[ncnt][t]), (PMT_y[ncnt][i]-Centro_y[ncnt][t]), (PMT_z[ncnt][i]-Centro_z[ncnt][t])); if(length_mt>1.5*l_sr){ //generate unit vector in the direction of the MT //source of the MT p0x = Centro_x[ncnt][t]; p0y = Centro_y[ncnt][t]; p0z = Centro_z[ncnt][t]; //tip of the MT px = PMT_x[ncnt][i]; //x coordinate before growth py = PMT_y[ncnt][i]; //y coordinate before growth pz = PMT_z[ncnt][i]; //z coordinate before growth a1 = px-p0x; a2 = py-p0y; a3 = pz-p0z; r = norm(a1,a2,a3); //unit vector in the direction of the MT growth a1 /= r; a2 /= r; a3 /= r; PMT_x[ncnt][i] -= l_sr*a1; PMT_y[ncnt][i] -= l_sr*a2; PMT_z[ncnt][i] -= l_sr*a3; rnd = ran2(&idum); if(rnd<=P_res){ labelMT[ncnt][i]=1; //rescue the MT N_rescue += 1; }; } else { Generate_Single_MT(ncnt,i,t); }; }; }; }; } void Rotate(int i, double psi, double theta, double phi){ Rot_matrix_ch(i,1,phi,theta,psi,0.,0.,0.); Rot_matrix_ch(i,-1,phi,theta,psi,0.,0.,0.); Rot_matrix_ch(i,0,phi,theta,psi,0.,0.,0.); Rot_matrix_kt(i,-1,phi,theta,psi,0.,0.,0.); Rot_matrix_kt(i,1,phi,theta,psi,0.,0.,0.); Rot_matrix_kt(i,0,phi,theta,psi,0.,0.,0.); } void Translate(int i, double x0, double y0, double z0){ double x1,y1,z1, del_x,del_y,del_z; //present center of mass of the Ch.; x1 = PCHm_x[i]; y1 = PCHm_y[i]; z1 = PCHm_z[i]; //vector joining present COM to new COM (x1-->x0) del_x = x0-x1; del_y = y0-y1; del_z = z0-z1; //translation of Ch. within the nuclear sphere PCHu_x[i] += del_x; PCHu_y[i] += del_y; PCHu_z[i] += del_z; PCHd_x[i] += del_x; PCHd_y[i] += del_y; PCHd_z[i] += del_z; PCHm_x[i] += del_x; PCHm_y[i] += del_y; PCHm_z[i] += del_z; //translation of KT within the nuclear sphere PKTl_x[i] += del_x; PKTl_y[i] += del_y; PKTl_z[i] += del_z; PKTl0_x[i] += del_x; PKTl0_y[i] += del_y; PKTl0_z[i] += del_z; PKTr_x[i] += del_x; PKTr_y[i] += del_y; PKTr_z[i] += del_z; PKTr0_x[i] += del_x; PKTr0_y[i] += del_y; PKTr0_z[i] += del_z; PKTmL_x[i] += del_x; PKTmL_y[i] += del_y; PKTmL_z[i] += del_z; PKTuL_x[i] += del_x; PKTuL_y[i] += del_y; PKTuL_z[i] += del_z; PKTdL_x[i] += del_x; PKTdL_y[i] += del_y; PKTdL_z[i] += del_z; PKTmR_x[i] += del_x; PKTmR_y[i] += del_y; PKTmR_z[i] += del_z; PKTuR_x[i] += del_x; PKTuR_y[i] += del_y; PKTuR_z[i] += del_z; PKTdR_x[i] += del_x; PKTdR_y[i] += del_y; PKTdR_z[i] += del_z; } void Move_chromosome(int t){ int i,j; double x,y,z,x0,y0,z0, d=5.0*R_nuc, phi,theta,psi; for(i=1;i<=N_ch;i++){ do{ int flag_within_vol = 0; do{ x=MP_width_x*(0.5-ran2(&idum))*Ax_out; y=MP_width_y*(0.5-ran2(&idum))*Ay_out; //metaphase plate width MP_width (0-2) z=MP_width_z*(0.5-ran2(&idum))*Az_out; int pointisin=OutOfSpheroid(x,y,z,Ax_out,Ay_out,Az_out); int pointisout=OutOfSpheroid(x,y,z,Ax_in,Ay_in,Az_in); if((pointisin==0)&&(pointisout==1))flag_within_vol = 1; }while(flag_within_vol==0); //center of the chromosome x_tr[i]=x; y_tr[i]=y; z_tr[i]=z; //distribute the CH/KT such that they have minimal overlap d=5.0*R_nuc;//large reference distance for(j=1;j<=i-1;j++){ double d1 = norm((x_tr[i]-x_tr[j]),(y_tr[i]-y_tr[j]),(z_tr[i]-z_tr[j])); if(d1=1000){ //try 1000 times to fit them flag_within_vol = 1; }; }while(flag_within_vol==0); //distribute the CH/KT such that they have minimal overlap d=5.0*R_nuc;//large reference distance for(j=1;j<=i-1;j++){ double d1 = norm((x-x_tr[j]),(y-y_tr[j]),(z-z_tr[j])); if(d1= 1000){ d=d_overlap; }; }while(d=0 and (P-P2).(P1-P2)>=0 and shortest distance of P from kt axis < R_kt If P is not within CH vol., check if the line PoldP has gone through the CH. Check if P */ double pole_pole_dist = norm(Centro_x[1][t]-Centro_x[2][t], Centro_y[1][t]-Centro_y[2][t], Centro_z[1][t]-Centro_z[2][t]); double prob_fn = 1.0; //attachement probability int flag_capture =0; for(i=1;i<=N_kt;i++){ if(labelMT[ncnt][Ni]==1){ //Check for all left KT if(label_kt_l_cir[i]==0){ //If KT is ring Like x1 = PKTuL_x[i]; y1 = PKTuL_y[i]; z1=PKTuL_z[i]; x2 = PKTdL_x[i]; y2 = PKTdL_y[i]; z2=PKTdL_z[i]; //If (P-P1).(P2-P1)>0, i.e. MT tip is between two ends of the KT PP1_P2P1 = (x-x1)*(x2-x1)+(y-y1)*(y2-y1)+(z-z1)*(z2-z1); PP2_P1P2 = (x-x2)*(x1-x2)+(y-y2)*(y1-y2)+(z-z2)*(z1-z2); if((PP1_P2P1>=0.0) && (PP2_P1P2>=0.0)){ //Find the shortest dist. from P to the kt axis. d = Shortest_Distance(x1,y1,z1,x2,y2,z2,x,y,z); if((d<=R_kt_cir)){// && (d>=R_ch)){ //i.e. MT end is within KT cylinder; if((KT_off_l[i] == 0) && (Nkt_cap_l[i]0, i.e. MT tip is between two ends of the KT PP1_P2P1 = (x-x1)*(x2-x1)+(y-y1)*(y2-y1)+(z-z1)*(z2-z1); PP2_P1P2 = (x-x2)*(x1-x2)+(y-y2)*(y1-y2)+(z-z2)*(z1-z2); if((PP1_P2P1>=0.0) && (PP2_P1P2>=0.0)){ //find the shortest dist. from P to the kt axis. d = Shortest_Distance(x1,y1,z1,x2,y2,z2,x,y,z); if(d<=R_kt){ //i.e. MT end is within KT cylinder; if((KT_off_l[i] == 0) && (Nkt_cap_l[i]0, i.e. MT tip is between two ends of the KT PP1_P2P1 = (x-x1)*(x2-x1)+(y-y1)*(y2-y1)+(z-z1)*(z2-z1); PP2_P1P2 = (x-x2)*(x1-x2)+(y-y2)*(y1-y2)+(z-z2)*(z1-z2); if((PP1_P2P1>=0.0) && (PP2_P1P2>=0.0)){ //Find the shortest dist. from P to the kt axis. d = Shortest_Distance(x1,y1,z1,x2,y2,z2,x,y,z); if((d<=R_kt_cir)){// && (d>=R_ch)){ //i.e. MT end is within KT cylinder; if((KT_off_r[i] == 0) && (Nkt_cap_r[i]=0.0) && (PP2_P1P2>=0.0)){ //find the shortest dist. from P to the kt axis. d = Shortest_Distance(x1,y1,z1,x2,y2,z2,x,y,z); if(d0 and (P-P2).(P1-P2)> and shortest distance of P from CH axis < R_ch If P is not within CH vol., check if the line PoldP has gone through the CH. Check if P */ int flag_capture =0; for(i=1;i<=N_ch;i++){ //check for all chromosomes x1 = PCHu_x[i]; y1 = PCHu_y[i]; z1=PCHu_z[i]; x2 = PCHd_x[i]; y2 = PCHd_y[i]; z2=PCHd_z[i]; x3 = PCHm_x[i]; y3 = PCHm_y[i]; z3=PCHm_z[i]; //find if (P-P1).(P2-P1)>0 double PP1_P2P1 = (x-x1)*(x2-x1)+(y-y1)*(y2-y1)+(z-z1)*(z2-z1); double PP2_P1P2 = (x-x2)*(x1-x2)+(y-y2)*(y1-y2)+(z-z2)*(z1-z2); if((PP1_P2P1>=0.0) && (PP2_P1P2>=0.0)){ //find the shortest dist. from P to the kt axis. d = Shortest_Distance(x1,y1,z1,x2,y2,z2,x,y,z); if(d<=R_ch){ //MT hit CH - should shrink flag_capture = 1; break; }; } }; return(flag_capture); } double DistP1P2(double x1, double y1, double z1, double x2, double y2, double z2){ double dst = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2)); return(dst); } int Check_all_kt_capture(int t){ int i, sum=0, cap_label =0; for (i=1;i<=N_kt;i++){ sum += label_kt_l[i]+label_kt_r[i]; //printf("sum+%d\n",sum); }; if(sum == 2*N_kt) cap_label = 1; return(cap_label); } void Calc_time_after_first_attach(int t){ int i; for(i=1;i<=N_kt;i++){ if(label_kt_l[i] == 1)KT_t_l[i] += 1; if(KT_t_l[i]>=MAX_KT_t) { rotlab_l[i] = 1; }; if(label_kt_r[i] == 1)KT_t_r[i] += 1; if(KT_t_r[i]>=MAX_KT_t) { rotlab_r[i] = 1; }; }; } void ChangeKTShapeCylinder(int t){ int i,j; double x0,y0,z0,phi,theta,psi; int T_KT_max; //change the shape of the kT when the time after capture is bigger than segT for(i=1;i<=N_ch;i++){ T_KT_max = KT_t_l[i]; if(T_KT_max == segT){ //set KT segregation lable to 1 label_kt_l_cir[i]=1; }; T_KT_max = KT_t_r[i]; if(T_KT_max == segT){ //set KT segregation lable to 1 label_kt_r_cir[i]=1; }; }; } void Rotate_after_first_attach(int t, int n_ch){ int i; double x,y,z; for(i=1;i<=N_ch;i++){ RotateSingleChrom(t,i); }; } void RotateSingleChrom(int t, int i){ // this is a geometric algorithm to rotate the KT // in the direction of the captured pole. int j,x,y,LR; double x0,y0,z0, d=5.0*R_nuc, phi,theta,psi; double p0x,p0y,p0z,p1x,p1y,p1z,p2x,p2y,p2z,p3x,p3y,p3z,p4x,p4y,p4z; double a1,a2,a3,b1,b2,b3,c1,c2,c3,r; double x_ori, y_ori, z_ori; double del_x, del_y, del_z; double pCx,pCy,pCz; // p0 // \ // \ p5 // \ | | // \p1 | | // ----| |---- // ----|p2--p3|---- p4 // | | // | | // p6 x=0; y=0; if((rotlab_l[i] == 1)){ //check if not a mixed attachment x = 1; //not a mixed one y = KT_l_pole[i][1]; //sign of the attached pole(1,2,3,4) for(j=2;j<=Nkt_cap_l[i];j++){ if(KT_l_pole[i][1] != KT_l_pole[i][j]) { x=2; //mixed y = 0; }; }; if(x==1){ for(j=1;j<=Nkt_cap_r[i];j++){ if(KT_l_pole[i][1] == KT_r_pole[i][j]){ x=2; //mixed y = 0; }; }; }; if(x==1){ //that is captured by a single pole Rot_flag = 1; //If Rot_flag=0, then rotation has not happened //find unit vector from p_0 to center of chromosome (pC) p0x = Centro_x[y][t]; p0y = Centro_y[y][t]; p0z = Centro_z[y][t]; pCx = PCHm_x[i]; pCy = PCHm_y[i]; pCz = PCHm_z[i]; //draw unit vector p0-->pC double p0pCx = pCx-p0x; double p0pCy = pCy-p0y; double p0pCz = pCz-p0z; r = norm(p0pCx,p0pCy,p0pCz); //unit vector a1 = p0pCx/r; a2 = p0pCy/r; a3 = p0pCz/r; //rotate p1-->p2 //distance between p1-->p2 p1x = PKTl_x[i]; p1y = PKTl_y[i]; p1z = PKTl_z[i]; p2x = PKTl0_x[i]; p2y = PKTl0_y[i]; p2z = PKTl0_z[i]; double p1p2x= p2x-p1x; double p1p2y= p2y-p1y; double p1p2z= p2z-p1z; r = norm(p1p2x,p1p2y,p1p2z); PKTl0_x[i]=PKTl_x[i]+r*a1; PKTl0_y[i]=PKTl_y[i]+r*a2; PKTl0_z[i]=PKTl_z[i]+r*a3; //rotate p2-->p3 //distance between p2-->p3 p3x = PKTr0_x[i]; p3y = PKTr0_y[i]; p3z = PKTr0_z[i]; double p2p3x= p3x-p2x; double p2p3y = p3y-p2y; double p2p3z = p3z-p2z; r = norm(p2p3x,p2p3y,p2p3z); PKTr0_x[i]=PKTl0_x[i]+r*a1;PKTr0_y[i]=PKTl0_y[i]+r*a2;PKTr0_z[i]=PKTl0_z[i]+r*a3; //rotate p3-->p4 //distance between p2-->p3 p4x = PKTr_x[i]; p4y = PKTr_y[i]; p4z = PKTr_z[i]; double p3p4x= p4x-p3x; double p3p4y = p4y-p3y; double p3p4z = p4z-p3z; r = norm(p3p4x,p3p4y,p3p4z); PKTr_x[i]=PKTr0_x[i]+r*a1;PKTr_y[i]=PKTr0_y[i]+r*a2;PKTr_z[i]=PKTr0_z[i]+r*a3; //ROTATE CH and RING KT //find unit normal vector to a1,a2,a3 c1 = 0.5+0.001*ran2(&idum); c2 = 0.5+0.001*ran2(&idum); c3 = -(a1*c1+a2*c2)/a3; r = norm(c1,c2,c3); b1 = c1/r; b2=c2/r; b3=c3/r; //mid point of CH after rotation; PCHm_x[i]=(PKTl0_x[i]+PKTr0_x[i])/2.0; PCHm_y[i]=(PKTl0_y[i]+PKTr0_y[i])/2.0; PCHm_z[i]=(PKTl0_z[i]+PKTr0_z[i])/2.0; //up-point of CH after rotation; PCHu_x[i]= PCHm_x[i]+ b1*l_ch/2.0; PCHu_y[i]= PCHm_y[i]+ b2*l_ch/2.0; PCHu_z[i]= PCHm_z[i]+ b3*l_ch/2.0; //down-point of CH after rotation; PCHd_x[i]= PCHm_x[i]-b1*l_ch/2.0; PCHd_y[i]= PCHm_y[i]-b2*l_ch/2.0; PCHd_z[i]= PCHm_z[i]-b3*l_ch/2.0; //Find the center of the left and right ring KT double dist_kt_tip = norm(PKTl_x[i]-PKTl0_x[i],PKTl_y[i]-PKTl0_y[i],PKTl_z[i]-PKTl0_z[i]); PKTmL_x[i] = PCHm_x[i] + d_kt_shift*(PKTl_x[i]-PKTl0_x[i])/dist_kt_tip; PKTmL_y[i] = PCHm_y[i] + d_kt_shift*(PKTl_y[i]-PKTl0_y[i])/dist_kt_tip; PKTmL_z[i] = PCHm_z[i] + d_kt_shift*(PKTl_z[i]-PKTl0_z[i])/dist_kt_tip; dist_kt_tip = norm(PKTr_x[i]-PKTr0_x[i],PKTr_y[i]-PKTr0_y[i],PKTr_z[i]-PKTr0_z[i]); PKTmR_x[i] = PCHm_x[i] + d_kt_shift*(PKTr_x[i]-PKTr0_x[i])/dist_kt_tip; PKTmR_y[i] = PCHm_y[i] + d_kt_shift*(PKTr_y[i]-PKTr0_y[i])/dist_kt_tip; PKTmR_z[i] = PCHm_z[i] + d_kt_shift*(PKTr_z[i]-PKTr0_z[i])/dist_kt_tip; //left circular KT rotation //up-point of KT after rotation; PKTuL_x[i]= PKTmL_x[i]+b1*l_kt_cir/2.0; //check the middle position PKTuL_y[i]= PKTmL_y[i]+b2*l_kt_cir/2.0; PKTuL_z[i]= PKTmL_z[i]+b3*l_kt_cir/2.0; //down-point of KT after rotation; PKTdL_x[i]= PKTmL_x[i]-b1*l_kt_cir/2.0; PKTdL_y[i]= PKTmL_y[i]-b2*l_kt_cir/2.0; PKTdL_z[i]= PKTmL_z[i]-b3*l_kt_cir/2.0; //right circular KT rotation //up-point of KT after rotation; PKTuR_x[i]= PKTmR_x[i]+b1*l_kt_cir/2.0; //check the middle position PKTuR_y[i]= PKTmR_y[i]+b2*l_kt_cir/2.0; PKTuR_z[i]= PKTmR_z[i]+b3*l_kt_cir/2.0; //down-point of KT after rotation; PKTdR_x[i]= PKTmR_x[i]-b1*l_kt_cir/2.0; PKTdR_y[i]= PKTmR_y[i]-b2*l_kt_cir/2.0; PKTdR_z[i]= PKTmR_z[i]-b3*l_kt_cir/2.0; //now translate everything back to the original position del_x=x_tr[i]-PCHm_x[i]; del_y=y_tr[i]-PCHm_y[i]; del_z=z_tr[i]-PCHm_z[i]; PKTl_x[i] += del_x;PKTl_y[i] += del_y;PKTl_z[i] += del_z; PKTl0_x[i] += del_x;PKTl0_y[i] += del_y;PKTl0_z[i] += del_z; PKTr0_x[i] += del_x;PKTr0_y[i] += del_y;PKTr0_z[i] += del_z; PKTr_x[i] += del_x;PKTr_y[i] += del_y;PKTr_z[i] += del_z; PCHm_x[i] += del_x; PCHm_y[i] += del_y; PCHm_z[i] += del_z; PCHu_x[i] += del_x; PCHu_y[i] += del_y; PCHu_z[i] += del_z; PCHd_x[i] += del_x; PCHd_y[i] += del_y; PCHd_z[i] += del_z; PKTmL_x[i] += del_x; PKTmL_y[i] += del_y; PKTmL_z[i] += del_z; PKTuL_x[i] += del_x; PKTuL_y[i] += del_y; PKTuL_z[i] += del_z; PKTdL_x[i] += del_x; PKTdL_y[i] += del_y; PKTdL_z[i] += del_z; PKTmR_x[i] += del_x; PKTmR_y[i] += del_y; PKTmR_z[i] += del_z; PKTuR_x[i] += del_x; PKTuR_y[i] += del_y; PKTuR_z[i] += del_z; PKTdR_x[i] += del_x; PKTdR_y[i] += del_y; PKTdR_z[i] += del_z; }; }; if((rotlab_r[i]==1)){ //check if not a mixed attachment x = 1; //not a mixed one y = KT_r_pole[i][1]; //sign of the attached pole for(j=2;j<=Nkt_cap_r[i];j++){ if(KT_r_pole[i][1] != KT_r_pole[i][j]){ x=2; //mixed y = 0; }; }; if(x==1){ for(j=1;j<=Nkt_cap_l[i];j++){ if(KT_r_pole[i][1] == KT_l_pole[i][j]){ x=2; //mixed y = 0; }; }; }; if(x==1){ Rot_flag = 1; //If Rot_flag=0, then rotation has not happened //find unit vector from p_0 to center of chromosome (pC) p0x = Centro_x[y][t]; p0y = Centro_y[y][t]; p0z = Centro_z[y][t]; pCx = PCHm_x[i]; pCy = PCHm_y[i]; pCz = PCHm_z[i]; //draw unit vector p0-->pC double p0pCx = pCx-p0x; double p0pCy = pCy-p0y; double p0pCz = pCz-p0z; r = norm(p0pCx,p0pCy,p0pCz); //unit vector a1 = p0pCx/r; a2 = p0pCy/r; a3 = p0pCz/r; //rotate p1-->p2 //distance between p1-->p2 p1x = PKTr_x[i]; p1y = PKTr_y[i]; p1z = PKTr_z[i]; p2x = PKTr0_x[i]; p2y = PKTr0_y[i]; p2z = PKTr0_z[i]; double p1p2x= p2x-p1x; double p1p2y = p2y-p1y; double p1p2z = p2z-p1z; r = norm(p1p2x,p1p2y,p1p2z); PKTr0_x[i] = PKTr_x[i] + r*a1; PKTr0_y[i] = PKTr_y[i] + r*a2; PKTr0_z[i] = PKTr_z[i] + r*a3; //rotate p2-->p3 //distance between p2-->p3 p3x = PKTl0_x[i]; p3y = PKTl0_y[i]; p3z = PKTl0_z[i]; double p2p3x= p3x-p2x; double p2p3y = p3y-p2y; double p2p3z = p3z-p2z; r = norm(p2p3x,p2p3y,p2p3z); PKTl0_x[i] = PKTr0_x[i] + r*a1;PKTl0_y[i] = PKTr0_y[i] + r*a2;PKTl0_z[i] = PKTr0_z[i] + r*a3; //rotate p3-->p4 //distance between p2-->p3 p4x = PKTl_x[i]; p4y = PKTl_y[i]; p4z = PKTl_z[i]; double p3p4x= p4x-p3x; double p3p4y = p4y-p3y; double p3p4z = p4z-p3z; r = norm(p3p4x,p3p4y,p3p4z); PKTl_x[i] = PKTl0_x[i] + r*a1;PKTl_y[i] = PKTl0_y[i] + r*a2;PKTl_z[i] = PKTl0_z[i] + r*a3; //ROTATE CH //find unit normal vector to a1,a2,a3 c1 = 0.5+0.001*ran2(&idum); c2 = 0.5+0.001*ran2(&idum); c3 = -(a1*c1+a2*c2)/a3; r = norm(c1,c2,c3); b1 = c1/r; b2=c2/r; b3=c3/r; //mid point of CH after rotation; PCHm_x[i] = (PKTl0_x[i] + PKTr0_x[i])/2.0; PCHm_y[i] = (PKTl0_y[i] + PKTr0_y[i])/2.0; PCHm_z[i] = (PKTl0_z[i] + PKTr0_z[i])/2.0; //up-point of CH after rotation; PCHu_x[i] = PCHm_x[i] + b1*l_ch/2.0; PCHu_y[i] = PCHm_y[i] + b2*l_ch/2.0; PCHu_z[i] = PCHm_z[i] + b3*l_ch/2.0; //down-point of CH after rotation; PCHd_x[i] = PCHm_x[i] - b1*l_ch/2.0; PCHd_y[i] = PCHm_y[i] - b2*l_ch/2.0; PCHd_z[i] = PCHm_z[i] - b3*l_ch/2.0; //Find the center of the left and right ring KT double dist_kt_tip = norm(PKTl_x[i]-PKTl0_x[i],PKTl_y[i]-PKTl0_y[i],PKTl_z[i]-PKTl0_z[i]); PKTmL_x[i] = PCHm_x[i] + d_kt_shift*(PKTl_x[i]-PKTl0_x[i])/dist_kt_tip; PKTmL_y[i] = PCHm_y[i] + d_kt_shift*(PKTl_y[i]-PKTl0_y[i])/dist_kt_tip; PKTmL_z[i] = PCHm_z[i] + d_kt_shift*(PKTl_z[i]-PKTl0_z[i])/dist_kt_tip; dist_kt_tip = norm(PKTr_x[i]-PKTr0_x[i],PKTr_y[i]-PKTr0_y[i],PKTr_z[i]-PKTr0_z[i]); PKTmR_x[i] = PCHm_x[i] + d_kt_shift*(PKTr_x[i]-PKTr0_x[i])/dist_kt_tip; PKTmR_y[i] = PCHm_y[i] + d_kt_shift*(PKTr_y[i]-PKTr0_y[i])/dist_kt_tip; PKTmR_z[i] = PCHm_z[i] + d_kt_shift*(PKTr_z[i]-PKTr0_z[i])/dist_kt_tip; //left circular KT rotation //up-point of KT after rotation; PKTuL_x[i]= PKTmL_x[i]+b1*l_kt_cir/2.0; //check the middle position PKTuL_y[i]= PKTmL_y[i]+b2*l_kt_cir/2.0; PKTuL_z[i]= PKTmL_z[i]+b3*l_kt_cir/2.0; //down-point of KT after rotation; PKTdL_x[i]= PKTmL_x[i]-b1*l_kt_cir/2.0; PKTdL_y[i]= PKTmL_y[i]-b2*l_kt_cir/2.0; PKTdL_z[i]= PKTmL_z[i]-b3*l_kt_cir/2.0; //right circular KT rotation //up-point of KT after rotation; PKTuR_x[i]= PKTmR_x[i]+b1*l_kt_cir/2.0; //check the middle position PKTuR_y[i]= PKTmR_y[i]+b2*l_kt_cir/2.0; PKTuR_z[i]= PKTmR_z[i]+b3*l_kt_cir/2.0; //down-point of KT after rotation; PKTdR_x[i]= PKTmR_x[i]-b1*l_kt_cir/2.0; PKTdR_y[i]= PKTmR_y[i]-b2*l_kt_cir/2.0; PKTdR_z[i]= PKTmR_z[i]-b3*l_kt_cir/2.0; //now translate everything back to the original position del_x=x_tr[i]-PCHm_x[i]; del_y=y_tr[i]-PCHm_y[i]; del_z=z_tr[i]-PCHm_z[i]; PKTl_x[i] += del_x;PKTl_y[i] += del_y;PKTl_z[i] += del_z; PKTl0_x[i] += del_x;PKTl0_y[i] += del_y;PKTl0_z[i] += del_z; PKTr0_x[i] += del_x;PKTr0_y[i] += del_y;PKTr0_z[i] += del_z; PKTr_x[i] += del_x;PKTr_y[i] += del_y;PKTr_z[i] += del_z; PCHm_x[i] += del_x; PCHm_y[i] += del_y; PCHm_z[i] += del_z; PCHu_x[i] += del_x; PCHu_y[i] += del_y; PCHu_z[i] += del_z; PCHd_x[i] += del_x; PCHd_y[i] += del_y; PCHd_z[i] += del_z; PKTmL_x[i] += del_x; PKTmL_y[i] += del_y; PKTmL_z[i] += del_z; PKTuL_x[i] += del_x; PKTuL_y[i] += del_y; PKTuL_z[i] += del_z; PKTdL_x[i] += del_x; PKTdL_y[i] += del_y; PKTdL_z[i] += del_z; PKTmR_x[i] += del_x; PKTmR_y[i] += del_y; PKTmR_z[i] += del_z; PKTuR_x[i] += del_x; PKTuR_y[i] += del_y; PKTuR_z[i] += del_z; PKTdR_x[i] += del_x; PKTdR_y[i] += del_y; PKTdR_z[i] += del_z; }; }; } void MoveCentrosome(int t, int f_v){ double x,y,z; //f_v is the velocity factor, =1 if moving else 0. Centro_x[1][t] = Centro_x[1][0]; Centro_y[1][t] = Centro_y[1][0]; Centro_z[1][t] = Centro_z[1][0]; SgnCent[1] = 1; Centro_x[2][t] = Centro_x[2][0]; Centro_y[2][t] = Centro_y[2][0]; Centro_z[2][t] = Centro_z[2][0]; SgnCent[2] = 2; } void Prevent_amphitelic(int t){ int i,j,x,y; for(i=1;i<=N_kt;i++){ //check left KT if(Nkt_cap_l[i] == 0){ x = 0; //monotelic, i.e. no attachment }else{ x = 1; //amphitelic, i.e., connected to one pole for(j=2;j<=Nkt_cap_l[i];j++) if(KT_l_pole[i][1] != KT_l_pole[i][j]){ x=2; //merotelic, i.e. connected to other poles }; }; //check right KT if(Nkt_cap_r[i] == 0){ y = 0; //monotelic, i.e. no attachment }else{ y = 1; //amphitelic, i.e., connected to one pole for(j=2;j<=Nkt_cap_r[i];j++) if(KT_r_pole[i][1] != KT_r_pole[i][j]){ y=2; //merotelic, i.e. connected to other poles }; }; //no error if((x==1) && (y==1)){ if(KT_l_pole[i][1] != KT_r_pole[i][1]){ //printf("x=%d y=%d t=%d\n",x,y,t); KT_off_l[i]=KT_off_r[i] = 1; }; }; }; } void Disolve_Mono(int i, int x, int y){ int j; if(x==1){ //left KT for(j=1;j<=Nkt_cap_l[i];j++)KT_l_pole[i][j] = 0; KT_off_l[i] = 0; Nkt_cap_l[i]=0; label_kt_l[i] = 0; KT_t_l[i] = 0; }; if(y==1){ //right KT for(j=1;j<=Nkt_cap_r[i];j++)KT_r_pole[i][j] = 0; KT_off_r[i] = 0; Nkt_cap_r[i]=0; label_kt_r[i] = 0; KT_t_r[i] = 0; }; } void Disolve_Synt(int i){ int j; if(ran2(&idum)<0.5){ //left KT for(j=1;j<=Nkt_cap_l[i];j++)KT_l_pole[i][j] = 0; KT_off_l[i] = 0; Nkt_cap_l[i]=0; label_kt_l[i] = 0; KT_t_l[i] = 0; }else{ //right KT for(j=1;j<=Nkt_cap_r[i];j++)KT_r_pole[i][j] = 0; KT_off_r[i] = 0; Nkt_cap_r[i]=0; label_kt_r[i] = 0; KT_t_r[i] = 0; }; } void Disolve_Mero(int i, int x, int y){ int j; if((x==2) && (y<2)){ for(j=1;j<=Nkt_cap_l[i];j++)KT_l_pole[i][j] = 0; KT_off_l[i] = 0; Nkt_cap_l[i]=0; label_kt_l[i] = 0; KT_t_l[i] = 0; }; if((x<2) && (y==2)){ for(j=1;j<=Nkt_cap_r[i];j++)KT_r_pole[i][j] = 0; KT_off_r[i] = 0; Nkt_cap_r[i]=0; label_kt_r[i] = 0; KT_t_r[i] = 0; }; if((x==2) && (y==2)){ if(ran2(&idum)<0.5){ for(j=1;j<=Nkt_cap_l[i];j++)KT_l_pole[i][j] = 0; KT_off_l[i] = 0; Nkt_cap_l[i]=0; label_kt_l[i] = 0; KT_t_l[i] = 0; }else{ for(j=1;j<=Nkt_cap_r[i];j++)KT_r_pole[i][j] = 0; KT_off_r[i] = 0; Nkt_cap_r[i]=0; label_kt_r[i] = 0; KT_t_r[i] = 0; }; }; } void Disolve_Amphi(int i){ int j; if(ran2(&idum)<0.5){ //left KT for(j=1;j<=Nkt_cap_l[i];j++)KT_l_pole[i][j] = 0; KT_off_l[i] = 0; Nkt_cap_l[i]=0; label_kt_l[i] = 0; KT_t_l[i] = 0; }else{ //right KT for(j=1;j<=Nkt_cap_r[i];j++)KT_r_pole[i][j] = 0; KT_off_r[i] = 0; Nkt_cap_r[i]=0; label_kt_r[i] = 0; KT_t_r[i] = 0; }; } void DetachBondswithRates(int t){ int i,j,x,y,z1,z2,sum=0; //calculate attachment errors for the left KTs // when a MT hits KT, KT_l_pole or KT_r_pole is added // the value +1 or -1 w.r.t. given pole for(i=1;i<=N_kt;i++){ //check left KT if(Nkt_cap_l[i] == 0){ x = 0; //monotelic, i.e. no attachment z1 = 0; }else{ x = 1; //amphitelic, i.e., connected to one pole z1 = KT_l_pole[i][1]; //sign of the attached pole(-1,1,2,-2) for(j=2;j<=Nkt_cap_l[i];j++) if(KT_l_pole[i][1] != KT_l_pole[i][j]){ x=2; //merotelic, i.e. connected to other poles z1 = 3; //mixed }; }; //check right KT if(Nkt_cap_r[i] == 0){ y = 0; //monotelic, i.e. no attachment z2=0; }else{ y = 1; //amphitelic, i.e., connected to one pole z2 = KT_r_pole[i][1]; //sign of the attached pole(-1,1,2,-2) for(j=2;j<=Nkt_cap_r[i];j++) if(KT_r_pole[i][1] != KT_r_pole[i][j]){ y=2; //merotelic, i.e. connected to other poles z2 = 3; //mixed }; }; //Disassemble monotelic if(((x==0) && (y==1)) || ((x==1) && (y==0))) { if(koff_mono>=ran2(&idum)){ Disolve_Mono(i,x,y); }; }; //Disassemble syntelic if((x*y==1) && (z1==z2)){ if(koff_synt>=ran2(&idum)){ Disolve_Synt(i); }; }; double pole_pole_dist = norm(Centro_x[1][t]-Centro_x[2][t], Centro_y[1][t]-Centro_y[2][t], Centro_z[1][t]-Centro_z[2][t]); double dst = 3.0 + 1.0*(0.5-ran2(&idum)); //Disassemble amphitelic if(pole_pole_dist=ran2(&idum)){ Disolve_Amphi(i); }; }; }; }; //Disassemble merotelic if(pole_pole_dist=ran2(&idum)){ Disolve_Mero(i,x,y); }; }; }; }; } void DisolveSyntMero(int t){ int i,j,x,y,z1,z2,sum=0; //calculate attachment errors for the left KTs // when a MT hits KT, KT_l_pole or KT_r_pole is added // the value +1 or -1 w.r.t. given pole for(i=1;i<=N_kt;i++){ //check left KT if(Nkt_cap_l[i] == 0){ x = 0; //monotelic, i.e. no attachment z1 = 0; }else{ x = 1; //amphitelic, i.e., connected to one pole z1 = KT_l_pole[i][1]; //sign of the attached pole(-1,1,2,-2) for(j=2;j<=Nkt_cap_l[i];j++) if(KT_l_pole[i][1] != KT_l_pole[i][j]){ x=2; //merotelic, i.e. connected to other poles z1 = 3; //mixed }; }; //check right KT if(Nkt_cap_r[i] == 0){ y = 0; //monotelic, i.e. no attachment z2=0; }else{ y = 1; //amphitelic, i.e., connected to one pole z2 = KT_r_pole[i][1]; //sign of the attached pole(-1,1,2,-2) for(j=2;j<=Nkt_cap_r[i];j++) if(KT_r_pole[i][1] != KT_r_pole[i][j]){ y=2; //merotelic, i.e. connected to other poles z2 = 3; //mixed }; }; //monotelic correction if(((x==0) && (y==1)) || ((x==1) && (y==0))) { KT_t_mono[i] += 1; if((KT_t_mono[i]>=MAX_disolv_t)) { KT_t_mono[i] = 0; Disolve_Mono(i,x,y); }; }; //syntelic correction if((x*y==1) && (z1==z2)){ KT_t_syn[i] += 1; if((KT_t_syn[i]>=MAX_disolv_t)) { KT_t_syn[i] = 0; Disolve_Synt(i); }; }; //merotelic correction if((x*y==2) || (x*y==4) || (((x==0) && (y==2))) || ((x==2) && (y==0))){ KT_t_mero[i] += 1; if((KT_t_mero[i]>=MAX_disolv_t)){ KT_t_mero[i] = 0; Disolve_Mero(i,x,y); }; }; }; } void CheckVisibility(int t){ int i,j,x,y,z1,z2,sum=0; int count=0,count_amphi=0,count_mero=0; //if chromosomes are visible from both poles then //they become either merotelic or amphitelic //Thus, visibility requires to find these two numbers. for(i=1;i<=N_kt;i++){ //check left KT if(Nkt_cap_l[i] == 0){ x = 0; //monotelic, i.e. no attachment z1 = 0; }else{ x = 1; //amphitelic, i.e., connected to one pole z1 = KT_l_pole[i][1]; //sign of the attached pole(-1,1,2,-2) for(j=2;j<=Nkt_cap_l[i];j++) if(KT_l_pole[i][1] != KT_l_pole[i][j]){ x=2; //merotelic, i.e. connected to other poles z1 = 3; //mixed }; }; //check right KT if(Nkt_cap_r[i] == 0){ y = 0; //monotelic, i.e. no attachment z2=0; }else{ y = 1; //amphitelic, i.e., connected to one pole z2 = KT_r_pole[i][1]; //sign of the attached pole(-1,1,2,-2) for(j=2;j<=Nkt_cap_r[i];j++) if(KT_r_pole[i][1] != KT_r_pole[i][j]){ y=2; //merotelic, i.e. connected to other poles z2 = 3; //mixed }; }; //find amphitelic if((x*y==1) && (z1 != z2)){ count_amphi++; //printf("i amphi=%d\n",i); }; //find merotelic if((x*y==2) || (x*y==4) || (((x==0) && (y==2))) || ((x==2) && (y==0))){ count_mero++; //printf("i mero=%d\n",i); }; }; count = count_amphi+count_mero; Capture_KT_NUM[t] += count; } void Calc_attach_error(int n_ch, int ensem){ /************* file handling ***********/ char fn[100],fn1[100],fn2[100],fn3[100],fn4[100]; FILE *fp,*fp1,*fp2,*fp3,*fp4; if(n_ch==Nch_for_hist){ sprintf(fn,"HistAmph2cnt_%dNch_%dNmt_%1.2fRkt_%1.2flkt_%1.2fRch_%1.2flch_%dtrot_%1.5fVcent.dat",n_ch,NMT[1],R_kt,l_kt,R_ch,l_ch, MAX_KT_t,V_centro); if(ensem==1) { fp = fopen(fn,"w"); }else fp = fopen(fn,"a"); sprintf(fn1,"HistMer2cnt_%dNch_%dNmt_%1.2fRkt_%1.2flkt_%1.2fRch_%1.2flch_%dtrot_%1.5fVcent.dat",n_ch,NMT[1],R_kt,l_kt,R_ch,l_ch, MAX_KT_t,V_centro); if(ensem==1) { fp1 = fopen(fn1,"w"); }else fp1 = fopen(fn1,"a"); sprintf(fn2,"HistSyn2cnt_%dNch_%dNmt_%1.2fRkt_%1.2flkt_%1.2fRch_%1.2flch_%dtrot_%1.5fVcent.dat",n_ch,NMT[1],R_kt,l_kt,R_ch,l_ch, MAX_KT_t,V_centro); if(ensem==1) { fp2 = fopen(fn2,"w"); }else fp2 = fopen(fn2,"a"); sprintf(fn3,"HistMon2cnt_%dNch_%dNmt_%1.2fRkt_%1.2flkt_%1.2fRch_%1.2flch_%dtrot_%1.5fVcent.dat",n_ch,NMT[1],R_kt,l_kt,R_ch,l_ch, MAX_KT_t,V_centro); if(ensem==1) { fp3 = fopen(fn3,"w"); }else fp3 = fopen(fn3,"a"); sprintf(fn4,"HistNoat2cnt_%dNch_%dNmt_%1.2fRkt_%1.2flkt_%1.2fRch_%1.2flch_%dtrot_%1.5fVcent.dat",n_ch,NMT[1],R_kt,l_kt,R_ch,l_ch, MAX_KT_t,V_centro); if(ensem==1) { fp4 = fopen(fn4,"w"); }else fp4 = fopen(fn4,"a"); }; /***********************/ int ErNone=0, ErMero=0, ErSyn=0, ErMon=0, ErNoat=0; int i,j,x,y,z1,z2,sum=0; //calculate attachment errors for the left KTs // when a MT hits KT, KT_l_pole or KT_r_pole is added // the value +1 or -1 w.r.t. given pole for(i=1;i<=N_kt;i++){ //check left KT if(Nkt_cap_l[i] == 0){ x = 0; //monotelic, i.e. no attachment z1 = 0; }else{ x = 1; //amphitelic, i.e., connected to one pole z1 = KT_l_pole[i][1]; //sign of the attached pole(-1,1,2,-2) for(j=2;j<=Nkt_cap_l[i];j++) if(KT_l_pole[i][1] != KT_l_pole[i][j]){ x=2; //merotelic, i.e. connected to other poles z1 = 3; //mixed }; }; //check right KT if(Nkt_cap_r[i] == 0){ y = 0; //monotelic, i.e. no attachment z2=0; }else{ y = 1; //amphitelic, i.e., connected to one pole z2 = KT_r_pole[i][1]; //sign of the attached pole(-1,1,2,-2) for(j=2;j<=Nkt_cap_r[i];j++) if(KT_r_pole[i][1] != KT_r_pole[i][j]){ y=2; //merotelic, i.e. connected to other poles z2 = 3; //mixed }; }; //no attachment if((x==0) && (y==0)){ Err_noatch[n_ch] += 1; ErNoat ++; }; //monotelic if((x==0) && (y==1)){ Err_mono[n_ch] += 1; ErMon ++; }; if((x==1) && (y==0)){ Err_mono[n_ch] += 1; ErMon ++; }; //syntelic if((x*y==1) && (z1==z2)){ Err_synt[n_ch] += 1; ErSyn ++; }; //merotelic if(x*y==2){ Err_mero[n_ch] += 1; ErMero ++; }; if(x*y==4){ Err_mero[n_ch] += 1; ErMero ++; }; if((x==0) && (y==2)){ Err_mero[n_ch] += 1; ErMero ++; }; if((x==2) && (y==0)){ Err_mero[n_ch] += 1; ErMero ++; }; //no error (amphitelic) if((x*y==1) && (z1 != z2)){ Err_none[n_ch] += 1; ErNone ++; }; }; if(n_ch==Nch_for_hist){ if(ErNone != 0){fprintf(fp,"%d\n",ErNone);}; if(ErMero != 0){fprintf(fp1,"%d\n",ErMero);}; if(ErSyn != 0){fprintf(fp2,"%d\n",ErSyn);}; if(ErMon != 0){fprintf(fp3,"%d\n",ErMon);}; if(ErNoat != 0){fprintf(fp4,"%d\n",ErNoat);} fclose(fp); fclose(fp1); fclose(fp2); fclose(fp3); fclose(fp4); }; } void CalculateFinalErr(int n_ch){ char fn[100]; FILE *fp; /***********************/ int ErNone=0, ErMero=0, ErSyn=0, ErMon=0, ErNoat=0; int i,j,x,y,z1,z2,sum=0; //calculate attachment errors for the left KTs // when a MT hits KT, KT_l_pole or KT_r_pole is added // the value +1 or -1 w.r.t. given pole for(i=1;i<=N_kt;i++){ //check left KT if(Nkt_cap_l[i] == 0){ x = 0; //monotelic, i.e. no attachment z1 = 0; }else{ x = 1; //amphitelic, i.e., connected to one pole z1 = KT_l_pole[i][1]; //sign of the attached pole(-1,1,2,-2) for(j=2;j<=Nkt_cap_l[i];j++) if(KT_l_pole[i][1] != KT_l_pole[i][j]){ x=2; //merotelic, i.e. connected to other poles z1 = 3; //mixed }; }; //check right KT if(Nkt_cap_r[i] == 0){ y = 0; //monotelic, i.e. no attachment z2=0; }else{ y = 1; //amphitelic, i.e., connected to one pole z2 = KT_r_pole[i][1]; //sign of the attached pole(-1,1,2,-2) for(j=2;j<=Nkt_cap_r[i];j++) if(KT_r_pole[i][1] != KT_r_pole[i][j]){ y=2; //merotelic, i.e. connected to other poles z2 = 3; //mixed }; }; //no attachment if((x==0) && (y==0)){ Err_noatch[n_ch] += 1; }; //monotelic if((x==0) && (y==1)){ Err_mono[n_ch] += 1; }; if((x==1) && (y==0)){ Err_mono[n_ch] += 1; }; //syntelic if((x*y==1) && (z1==z2)){ Err_synt[n_ch] += 1; }; //merotelic if(x*y==2){ Err_mero[n_ch] += 1; }; if(x*y==4){ Err_mero[n_ch] += 1; }; if((x==0) && (y==2)){ Err_mero[n_ch] += 1; }; if((x==2) && (y==0)){ Err_mero[n_ch] += 1; }; //no error (amphitelic) if((x*y==1) && (z1 != z2)){ Err_none[n_ch] += 1; }; }; } void ErrFreq_with_time(int n_ch, int t){ char fn[100]; FILE *fp; /***********************/ int ErNone=0, ErMero=0, ErSyn=0, ErMon=0, ErNoat=0; int i,j,x,y,z1,z2,sum=0; //calculate attachment errors for the left KTs // when a MT hits KT, KT_l_pole or KT_r_pole is added // the value +1 or -1 w.r.t. given pole for(i=1;i<=N_kt;i++){ //check left KT if(Nkt_cap_l[i] == 0){ x = 0; //monotelic, i.e. no attachment z1 = 0; }else{ x = 1; //amphitelic, i.e., connected to one pole z1 = KT_l_pole[i][1]; //sign of the attached pole(-1,1,2,-2) for(j=2;j<=Nkt_cap_l[i];j++) if(KT_l_pole[i][1] != KT_l_pole[i][j]){ x=2; //merotelic, i.e. connected to other poles z1 = 3; //mixed }; }; //check right KT if(Nkt_cap_r[i] == 0){ y = 0; //monotelic, i.e. no attachment z2=0; }else{ y = 1; //amphitelic, i.e., connected to one pole z2 = KT_r_pole[i][1]; //sign of the attached pole(-1,1,2,-2) for(j=2;j<=Nkt_cap_r[i];j++) if(KT_r_pole[i][1] != KT_r_pole[i][j]){ y=2; //merotelic, i.e. connected to other poles z2 = 3; //mixed }; }; //no attachment if((x==0) && (y==0)){ ErNoat ++; }; //monotelic if((x==0) && (y==1)){ ErMon ++; }; if((x==1) && (y==0)){ ErMon ++; }; //syntelic if((x*y==1) && (z1==z2)){ ErSyn ++; }; //merotelic if(x*y==2){ ErMero ++; }; if(x*y==4){ ErMero ++; }; if((x==0) && (y==2)){ ErMero ++; }; if((x==2) && (y==0)){ ErMero ++; }; //no error (amphitelic) if((x*y==1) && (z1 != z2)){ ErNone ++; }; }; if(n_ch==Nch_for_hist){ if(ErNone != 0){Count_ErNone[t] += ErNone;}; if(ErMero != 0){Count_ErMero[t] += ErMero;}; if(ErSyn != 0){Count_ErSyn[t] += ErSyn;}; if(ErMon != 0){Count_ErMon[t] += ErMon;}; if(ErNoat != 0){Count_ErNoat[t] += ErNoat;} }; } void Output_Errors(int n_ch, int ensem){ char fn[100]; FILE *fp; int i; double denom; sprintf(fn,"ErrFinal_%dNch_%dTcond_%dMT_%1.2fhkt_%1.2fwkt_%1.2fGap_%1.2fRktcir_%1.2fXoff_%dCflg_%dtrot.dat",n_ch,segT,NMT[1],hkt,wkt,arc_l_fin,R_kt_cir_fin,MP_width_x*Ax_out,label_seg,MAX_KT_t); if(n_ch==N_ch_min) { fp = fopen(fn,"w"); }else fp = fopen(fn,"a"); denom = (double)(ensem*n_ch)/100.0; fprintf(fp,"#Nch Amphi Mero Synt Synt+Mero\n"); fprintf(fp,"%d %f %f %f %f\n",n_ch,(double)Err_none[n_ch]/denom,(double)Err_mero[n_ch]/denom,(double)Err_synt[n_ch]/denom,(double)(Err_mero[n_ch]+Err_synt[n_ch])/denom); fclose(fp); } void Visualization_povray(int tt, int n_ch){ char fn[100]; FILE *fp; int i,j,s,x,y,t,ncnt; printf("povray -D %dsnp2cnt_%dNch_%dTcond_%dMT_%1.2fhkt_%1.2fwkt_%1.2fGap_%1.2fRktcir_%1.2fXoff_%dCflg_%dtrot.pov\n",100000+tt,n_ch,segT,NMT[1],hkt,wkt,arc_l_fin,R_kt_cir_fin,MP_width_x*Ax_out,label_seg,MAX_KT_t); sprintf(fn,"%dsnp2cnt_%dNch_%dTcond_%dMT_%1.2fhkt_%1.2fwkt_%1.2fGap_%1.2fRktcir_%1.2fXoff_%dCflg_%dtrot.pov",100000+tt,n_ch,segT,NMT[1],hkt,wkt,arc_l_fin,R_kt_cir_fin,MP_width_x*Ax_out,label_seg,MAX_KT_t); fp = fopen(fn,"w"); fprintf(fp,"background{color rgb <0, 0, 0>}\n"); fprintf(fp,"\n"); //centrosome location for(ncnt=1;ncnt<=NCentro_max;ncnt++){ fprintf(fp,"sphere{<%f, %f, %f>, %f\n",Centro_x[ncnt][tt],Centro_y[ncnt][tt],Centro_z[ncnt][tt],0.2); fprintf(fp,"\n"); if(ncnt==1)fprintf(fp,"pigment{color rgbt <0.9, 0.5, 0.5, 0.0>}\n"); if(ncnt==2)fprintf(fp,"pigment{color rgbt <0.5, 0.9, 0.5, 0.0>}\n"); fprintf(fp,"finish{ambient 0.5 specular 0.5}}\n"); }; fprintf(fp,"\n"); //microtubules from centrosome to kinetochore, shrinking, uncaptured for(ncnt=1;ncnt<=NCentro_max;ncnt++){ for ( i=1; i<=NMT[ncnt]; i++ ){ fprintf(fp,"cylinder{<%f,%f,%f>, <%f,%f,%f>, 0.01\n", Centro_x[ncnt][tt],Centro_y[ncnt][tt],Centro_z[ncnt][tt], PMT_x[ncnt][i],PMT_y[ncnt][i],PMT_z[ncnt][i]); if(ncnt==1)fprintf(fp,"pigment{color rgb <1.0, 0.0, 0.0>}\n"); if(ncnt==2)fprintf(fp,"pigment{color rgb <0.0, 1.0, 1.0>}\n"); fprintf(fp,"finish{ambient 0.5 specular 0.5}}\n"); }; }; fprintf(fp,"\n"); //captured MTs from centrosome for(i=1; i<=N_ch; i++ ){ for(j=1;j<=Nkt_cap_l[i];j++){ ncnt = KT_l_pole[i][j]; if(ncnt !=0) fprintf(fp,"cylinder{<%f,%f,%f>, <%f,%f,%f>, 0.05\n", Centro_x[ncnt][tt],Centro_y[ncnt][tt],Centro_z[ncnt][tt], (PKTl0_x[i]+PKTl_x[i])/2.0,(PKTl0_y[i]+PKTl_y[i])/2.0,(PKTl0_z[i]+PKTl_z[i])/2.0); if(ncnt==1)fprintf(fp,"pigment{color rgb <0.0, 1.0, 0.0>}\n"); if(ncnt==2)fprintf(fp,"pigment{color rgb <0.0, 1.0, 0.0>}\n"); fprintf(fp,"finish{ambient 0.5 specular 0.5}}\n"); }; for(j=1;j<=Nkt_cap_r[i];j++){ ncnt = KT_r_pole[i][j]; if(ncnt !=0) fprintf(fp,"cylinder{<%f,%f,%f>, <%f,%f,%f>, 0.05\n", Centro_x[ncnt][tt],Centro_y[ncnt][tt],Centro_z[ncnt][tt], (PKTr0_x[i]+PKTr_x[i])/2.0,(PKTr0_y[i]+PKTr_y[i])/2.0,(PKTr0_z[i]+PKTr_z[i])/2.0); if(ncnt==1)fprintf(fp,"pigment{color rgb <0.0, 1.0, 0.0>}\n"); if(ncnt==2)fprintf(fp,"pigment{color rgb <0.0, 1.0, 0.0>}\n"); fprintf(fp,"finish{ambient 0.5 specular 0.5}}\n"); }; }; //chromosome fprintf(fp,"union\n"); fprintf(fp,"{\n"); for ( i=1; i<=N_ch; i++ ){ fprintf(fp,"cylinder{<%f,%f,%f>, <%f,%f,%f>, %f}\n", PCHu_x[i],PCHu_y[i],PCHu_z[i],PCHd_x[i],PCHd_y[i],PCHd_z[i],R_ch); //temp }; fprintf(fp,"\n"); fprintf(fp,"pigment{color rgb <0.7, 0.7, 0.7>}\n"); fprintf(fp,"finish{ambient 0.5 specular 0.5 }\n"); fprintf(fp,"\n"); fprintf(fp,"}\n"); fprintf(fp,"\n"); //kinetochore for(i=1;i<=N_kt;i++){ if (label_kt_l_cir[i]==1){//i.e. KT has segregated //left KT fprintf(fp,"cylinder{<%f,%f,%f>, <%f,%f,%f>, %f\n", PKTl_x[i],PKTl_y[i],PKTl_z[i], PKTl0_x[i], PKTl0_y[i],PKTl0_z[i],R_kt); if(label_kt_l[i] == 1)fprintf(fp,"pigment{color rgb <0.5, 0.9, 0.1>}\n"); if(label_kt_l[i] == 0)fprintf(fp,"pigment{color rgb <0.1, 0.1, 0.9>}\n"); fprintf(fp,"finish{ambient 0.5 specular 0.5 }}\n"); }else{ fprintf(fp,"cylinder{<%f,%f,%f>, <%f,%f,%f>, %f\n", PKTuL_x[i],PKTuL_y[i],PKTuL_z[i], PKTdL_x[i],PKTdL_y[i],PKTdL_z[i],R_kt_cir); if(label_kt_l[i] == 1)fprintf(fp,"pigment{color rgb <0.5, 0.9, 0.1>}\n"); if(label_kt_l[i] == 0)fprintf(fp,"pigment{color rgb <0.1, 0.1, 0.9>}\n"); fprintf(fp,"finish{ambient 0.5 specular 0.5 }}\n"); }; if (label_kt_r_cir[i]==1){//i.e. KT has segregated //right KT fprintf(fp,"cylinder{<%f,%f,%f>, <%f,%f,%f>, %f\n", PKTr_x[i],PKTr_y[i],PKTr_z[i], PKTr0_x[i], PKTr0_y[i],PKTr0_z[i],R_kt); if(label_kt_r[i] == 1)fprintf(fp,"pigment{color rgb <0.5, 0.9, 0.1>}\n"); if(label_kt_r[i] == 0)fprintf(fp,"pigment{color rgb <0.1, 0.1, 0.9>}\n"); fprintf(fp,"finish{ambient 0.5 specular 0.5 }}\n"); }else{ fprintf(fp,"cylinder{<%f,%f,%f>, <%f,%f,%f>, %f\n", PKTuR_x[i],PKTuR_y[i],PKTuR_z[i], PKTdR_x[i],PKTdR_y[i],PKTdR_z[i],R_kt_cir); if(label_kt_r[i] == 1)fprintf(fp,"pigment{color rgb <0.5, 0.9, 0.1>}\n"); if(label_kt_r[i] == 0)fprintf(fp,"pigment{color rgb <0.1, 0.1, 0.9>}\n"); fprintf(fp,"finish{ambient 0.5 specular 0.5 }}\n"); }; }; //light source fprintf(fp,"\n"); fprintf(fp,"light_source\n"); fprintf(fp,"{\n"); fprintf(fp,"<-60, 70, %f> color rgb <1, 1, 1>\n",4.0*R_nuc); fprintf(fp,"area_light <-150, 0, 0>, <0, 150, 0>, 25, 25\n"); fprintf(fp,"jitter\n"); fprintf(fp,"adaptive 1\n"); fprintf(fp,"rotate <0, -degrees(atan2(sqrt(0.5), 1)), 0>\n"); fprintf(fp,"rotate <0, 0, 45>\n"); fprintf(fp,"rotate <0, 22.5, 0>\n"); fprintf(fp,"rotate <12, 0, 0>\n"); fprintf(fp,"\n"); fprintf(fp,"}\n"); // camera position fprintf(fp,"\n"); fprintf(fp,"camera\n"); fprintf(fp,"{\n"); //Aspect ratio of image fprintf(fp," up <0,1,0>\n"); fprintf(fp," right <1,0,0>\n"); fprintf(fp,"location <%f, %f, %f>\n", PCHm_x[1],PCHm_y[1],1.5*R_nuc); fprintf(fp,"look_at <%f, %f, %f>\n",0.0,0.0,0.0); fprintf(fp,"\n"); fprintf(fp,"}\n"); fclose(fp); } void Output_ch_kt_coordinate(){ int i,j; char fn[70]; FILE *fp; sprintf(fn,"pos_ch.dat"); fp = fopen(fn,"w"); for(i=1;i<=N_ch;i++){ fprintf(fp,"%f %f %f\n",PCHu_x[i],PCHu_x[i],PCHu_x[i]); fprintf(fp,"%f %f %f\n",PCHm_x[i],PCHm_x[i],PCHm_x[i]); fprintf(fp,"%f %f %f\n",PCHd_x[i],PCHd_x[i],PCHd_x[i]); fprintf(fp,"\n"); }; fclose(fp); sprintf(fn,"pos_kt.dat"); fp = fopen(fn,"w"); for(i=1;i<=N_kt;i++){ fprintf(fp,"%f %f %f\n",PKTl_x[i],PKTl_y[i],PKTl_z[i]); fprintf(fp,"%f %f %f\n",PKTr_x[i],PKTr_y[i],PKTr_z[i]); fprintf(fp,"\n"); }; fclose(fp); } void Output_cap_time(int p, int n_ch, int eff_ensem){ int i,j; char fn[100]; FILE *fp; sprintf(fn,"ktCapT_%dNch_%dTcond_%dMT_%1.2fhkt_%1.2fwkt_%1.2fGap_%1.2fRktcir_%1.2fXoff_%dCflg_%dtrot.dat",n_ch,segT,NMT[1],hkt,wkt,arc_l_fin,R_kt_cir_fin,MP_width_x*Ax_out,label_seg,MAX_KT_t); if(n_ch==N_ch_min) { fp = fopen(fn,"w"); }else fp = fopen(fn,"a"); double denom = (double)(eff_ensem)*60.0; if(T_cap[n_ch] !=0)fprintf(fp,"%d %f\n",n_ch,(double)T_cap[n_ch]/denom); fclose(fp); } void Calculate_avg_mt_length(int n_ch){ int count =0; double l_avg = 0.0; if(t_count_l !=0){ l_avg += l_av_mt_l/(double)t_count_l; count += 1; }; if(t_count_r !=0){ l_avg += l_av_mt_r/(double)t_count_r; count += 1; }; L_av[n_ch] += l_avg/(double)count; } void Success_failure(int n_ch){ int i,j; char fn[100]; FILE *fp; sprintf(fn,"suc_%dMT_%1.2fRkt_%1.2flkt_%1.2fRch_%1.2flch_%dtrot.dat",NMT[1],R_kt,l_kt,R_ch,l_ch,MAX_KT_t); if(n_ch==2) { fp = fopen(fn,"w"); }else fp = fopen(fn,"a"); fprintf(fp,"%d %f %f\n",n_ch,(double)num_succ[n_ch]/(double)MAX_ensem,(double)num_fail[n_ch]/(double)MAX_ensem); fclose(fp); } void Distrb_cap_time(int eff_ensem, int n_ch, int t){ char fn[100]; FILE *fp; sprintf(fn,"CapTdstrb_%dTcond_%dN_ch_%dtrot.dat",segT,n_ch,NMT[1],MAX_KT_t); if(eff_ensem==1){ fp = fopen(fn,"w"); }else fp = fopen(fn,"a"); fprintf(fp,"%d %f\n",eff_ensem, (double)t/60.0); fclose(fp); } void DetachBondswithAging(int t){ int i,j,x,y,z1,z2,sum=0; for(i=1;i<=N_kt;i++){ if(KT_off_l[i] == 0){ for(j=1;j<=Nkt_cap_l[i];j++){ //increase bond lifetime & //check for removal as per aging criterion if(KT_l_time[i][j] !=0){ KT_l_time[i][j] += 1; if(KT_l_time[i][j]>=MAX_disolv_t){ KT_l_pole[i][j] = 0; KT_l_pole[i][j] = KT_l_pole[i][Nkt_cap_l[i]]; KT_l_pole[i][Nkt_cap_l[i]] = 0; KT_l_time[i][j] = KT_l_time[i][Nkt_cap_l[i]]; KT_l_time[i][Nkt_cap_l[i]] = 0; Nkt_cap_l[i] -= 1; }; }; if(Nkt_cap_l[i] <=0){ KT_off_l[i] = 0; Nkt_cap_l[i]=0; label_kt_l[i] = 0; KT_t_l[i] = 0; break; }; }; }; if(KT_off_r[i] == 0){ for(j=1;j<=Nkt_cap_r[i];j++){ //increase bond lifetime & //check for removal as per aging criterion if(KT_r_time[i][j] !=0){ KT_r_time[i][j] += 1; if(KT_r_time[i][j]>=MAX_disolv_t){ KT_r_pole[i][j] = 0; KT_r_pole[i][j] = KT_r_pole[i][Nkt_cap_r[i]]; KT_r_pole[i][Nkt_cap_r[i]] = 0; KT_r_time[i][j] = KT_r_time[i][Nkt_cap_r[i]]; KT_r_time[i][Nkt_cap_r[i]] = 0; Nkt_cap_r[i] -= 1; }; }; if(Nkt_cap_r[i] <=0){ KT_off_r[i] = 0; Nkt_cap_r[i]=0; label_kt_r[i] = 0; KT_t_r[i] = 0; break; }; }; }; }; } /**********************************************/ /* MAIN PROGRAM */ /**********************************************/ int main(int argc, char** argv ) { int x, y, p, n_mt,n_ch,i,j,s,nb,x_p,y_p,ensem,run, cap_flag, eff_ensem=0,chk_flg=0; int input_flag=0, t_move,t0, MTperCentro,ncnt,t_ErrFreq; double div, f_cat, Rcent_ini; double factor_R_kt_cir,factor_d_kt_shift,factor_R_kt_cir_fin,factor_d_kt_shift_fin; for( i = 1; i0){ T_cap[n_ch] += i; cap_flag = 1; eff_ensem += 1; num_succ[n_ch] += 1; CalculateFinalErr(n_ch); break; }; }; if(chk_flg==0){ cap_flag = 1; num_fail[n_ch] += 1; }; }; Output_Errors(n_ch,MAX_ensem); Output_cap_time(p,n_ch,eff_ensem); }; }