Self-organization of developing embryo using scale-invariant approach

Background Self-organization is a fundamental feature of living organisms at all hierarchical levels from molecule to organ. It has also been documented in developing embryos. Methods In this study, a scale-invariant power law (SIPL) method has been used to study self-organization in developing embryos. The SIPL coefficient was calculated using a centro-axial skew symmetrical matrix (CSSM) generated by entering the components of the Cartesian coordinates; for each component, one CSSM was generated. A basic square matrix (BSM) was constructed and the determinant was calculated in order to estimate the SIPL coefficient. This was applied to developing C. elegans during early stages of embryogenesis. The power law property of the method was evaluated using the straight line and Koch curve and the results were consistent with fractal dimensions (fd). Diffusion-limited aggregation (DLA) was used to validate the SIPL method. Results and conclusion The fractal dimensions of both the straight line and Koch curve showed consistency with the SIPL coefficients, which indicated the power law behavior of the SIPL method. The results showed that the ABp sublineage had a higher SIPL coefficient than EMS, indicating that ABp is more organized than EMS. The fd determined using DLA was higher in ABp than in EMS and its value was consistent with type 1 cluster formation, while that in EMS was consistent with type 2.


Background
Self-organization is a property of the biological structure [1] and is reported to be important in protein folding [2][3][4]. It has also been documented that at higher hierarchical levels such as the organelle level, it has a crucial role in the biogenesis of secretory granules in the Golgi apparatus [5,6]. Martin and Russell have shown that self-organization exists in mitochondria, where redox reactions are localized [7]. The most obvious example of self-organization at the organelle level is the cytoskeleton during the mitotic cycle, where mitotic spindle forms dynamically [8] using molecular motors [9]. Misteli concluded that self-organization could govern the mechanistic principles of cellular architecture [10]. The multicellular embryo develops from a zygote, characterized by a dynamic self-organizing process [11].
At an early stage of embryonic development, the forming cells adhere to each other [12] with coordinated cellular movement to form the primary embryonic body axis [13]. These movements are self-regulated and lead to a defined pattern [14]. In vitro matrix. If the elements in the upper and lower triangular matrix of the matrix have equal values(a i,j = a j,i ), then it is a symmetric square matrix. If all the elements in the upper triangular matrix have negative values of the lower triangular matrix or vice versa (a i,j = -a j,i ), then it is a skewed symmetric (anti-symmetric) matrix, and if the diagonal elements values are equal to zero, the matrix is known as "zero centro-axial skewed-symmetric matrix" (CSSM). The matrix resulting from the exchange of the upper and lower triangular matrices is a transpose matrix. If the a square matrix is subtracted from its transpose, followed by division by two, then the resulting matrix is a skew matrix, while the sum of a square matrix and its transpose followed by division by two is a symmetric matrix. The square matrix is used for generating symmetric and anti-symmetric matrices. The square matrix generated in this study by a special algorithm is called the basic square matrix (BSM).

The scale invariant power law
There are two aspects of the scale invariant power law: scale invariance means that the value of the SIPL coefficient does not change as the scale [31], magnification [32], or tissue growth changes [33,34]; and a power law is a relationship between two variables where one quantity varies as a function of the power of the other [33]. For example, Zhang and Sejnowski revealed that the growth of the volume of the white matter increases disproportionately more quickly than the gray matter, where it follows a power law relation [35]. In fact, one of the properties of power laws is scale-invariance [33]. Therefore, the SIPL defines the coefficient obtained by calculating the BSM determinant, which follows a power law rule and is scale invariant.
The reason for using the power law is the nature of the biological matter, over 21 orders of magnitude consistently follows a simple and systematic empirical power law. This includes metabolic rate, time scales and body size [36]. The most commonly used power laws are fractal dimension and allometry [37]. Fractal dimensions have been used to study diverse structures in nature at different levels and from galaxies [38] to subatomic structures [39]. In biomedicine, there are wide ranging applications; for example, at the molecular level, fractals were proposed for evaluating the physical features of ion channel proteins [40]. Vélez et al. reported the possible use of multifractals in the measurement of local variations in DNA sequence in order to define the structure-function relationship in chromosomes [41], and Mathur et al. used fractal analysis of gene expression in studying the hair growth cycle. Moreover [42], fractal genomics modeling has been used to predict new factors in signaling pathways and the networks operating in neurodegenerative disorders [43]. At the cellular level, fractal dimension was used in evaluating the morphological diversity of neurons and discriminating them on the basis of the neuronal extensions [44]; fractals can also explain higher orders of organization in biological materials such as the organization of tissues [45] and branching of tubular systems such as the respiratory and the vascular systems [46][47][48][49]. On the other hand, one of the best known applications of allometry is the metabolic rate scale (Kleiber's law), which is considered universal among different species, within the same species, or in individual animal at different orders including molecular, cellular and body levels [50,51]. A similarly universal allometric law relating time and body weight, including growth rates and animal age, has been documented [52]. This time scale relation is noticed in development biology [53]. Gillooly et al. reported an allometric relationship between metabolic rate and the developmental growth rate during embryogenesis, which has phylogenically and ontogenically invariant values [54]. Allometry was recently used in pharmacokinetics [55], predicting the pharmacokinetics of drugs [56,57]. In addition to allometry and Kleiber's law, other investigators have reported power law relationships in biomedicine; for example, Grandison and Morris reported that kinetic rate parameters showed a scale free relationship with the gene network and protein-protein interactions, which follows Benford's law [58]. Also, Zipf's Law has been used to discriminate the effect of natural selection from random genetic drift [59]; Furusawa and Kaneko (2003) reported that Zipf's Law applies universally to gene expression in yeast, nematodes, mammalian embryonic stem cells and human tissues [60].
The above discussion suggests that not every power law is fractal; on the other hand, in certain situations the behavior of the system shows fractal-like properties but is not truly fractal [61]. In addition, even natural fractal structures such as the triadic Koch curve could have non-fractal properties [62]. The growth of differentiating cells in a developing embryo certainly follows a power law, so we are justified in calling it SIPL to avoid fallacious attribution of fractal properties.

Analytical descriptions of CSSM and BSM CSSM
Suppose we have n points {(x 1 , y 1 , z 1 ), (x 2 , y 2 , z 2 ), ..., (x n , y n , z n )} ⊂ R 3 and relative 1 × n matrices {(x 1 , x 2 , ..., x n ), (y 1 , y 2 , ..., y n ), ..., (z 1 , z 2 , ..., z n )}. By subtracting the first entry x 1 from x i , for each 1 ≤ i ≤ n, we get (0, x 2 -x 1 , ..., x n -x 1 ), with 0 as its first entry. We do this for the other matrices. Now we use each of the resulting matrices as the first row of the CSSM matrix. The other rows of the CSSM matrix are defined using the recursive formula: We only need to prove that the matrix which is constructed from (a 1 , a 2 , ..., a n ), where x i,j = a j -a i , is anti-symmetric and hence it is CSSM. For each i,j we have x n i,j = x n i−1,j − x n i−1,i. Let x i,j = a j -a i . Then x i,j = -(a i -a j ) = -x j,i , so the matrix is anti-symmetric. Now, we need to prove that when x i,j = a j -a i , then we have x i+1,j = a j -a i+1 .
We have Also, by definition, this holds for the first row. So for every i,j, x i,j = a j -a i This shows the matrix is anti-symmetric.

BSM
Now we prove that there is a one to one correspondence between the CSSM matrices and the BSM matrices which is generated from the CSSM matrices. Suppose that (A i,j ) is a CSSM matrix and the BSM matrix defined by

Now we show that
We know that A i,j A j,i < 0 or A i,j A j,i = 0. If A i,j = A j,i = 0, then B i,j = B j,i = 0 and the claim is true. If A i,j > 0, then B i,j = 4A i,j and B i,j = -2A j,i . This implies that The case A i,j < 0 is similar.

Descriptions of the biological data
In the early stages of C. elgans embryogenesis, the zygote (P0) divides into two daughter cells called the anterior blastomere (AB) and the posterior blastomere (P1) forming a 2-cell stage embryo. This is followed by a second round of mitosis, where AB divides into ABa (anterior) and ABp (posterior), while P1 divides into P2 and EMS forming a 4-cell stage embryo (see Figure 1). ABp differentiates into different types of cells including neurons, body muscle, excretory duct cell and hypodermis, while EMS differentiates into 42 body muscles and intestine [63]. During organogenesis of the C. elegans embryo, ABp differentiates into a nervous system and epidermis, while EMS differentiates into muscular tissues, midgut and pharynx [64]. Axis determination is one of the most important events in the early stages of C. elegans embryogenesis; as the pronucleus breaks down in the zygote, asymmetric division follows forming a large daughter cell (AB) and a smaller one (P1) establishing the first antero-posterior axis.  AB starts the next division, which is initially oriented orthogonally to the antero-posterior axis, but as the cell progresses through anaphase, the orientation of the mitotic spindle of the dividing cell skews, resulting in anterior position of ABa to ABp. P1 commences mitosis a few minutes later resulting in a large EMS progeny cell ventrally located, and a smaller P2 posteriorly located; this round of cell division establishes the dorso-ventral axis [65]. The other important event is garstulation, which begins at the 28-cell stage of development, where Ea and Ep move to the center of the developing embryo and gastrulate forming the three germ layers [64]. There are several reasons for comparing Abp-derived cells (ABp-dc) with those from EMS (EMS-dc). At the developmental level, ABp is derived from the AB blastomere while EMS is derived from the P1 blastomere [66], so ABp and EMS are two different lineages and the use of their cells is relevant in developmental biology. At the organogenesis level, ABp differentiates into nervous system and epidermis, while EMS differentiates into muscular tissue, midgut and pharynx [64], thus ABp and EMS form entirely different organs. At the cellular level, the cells derived from AB blastomere (ABa and ABp) enter the mitotic cycle and divide earlier than those of P1 (EMS and P2), while EMS enters the mitotic cycle earlier than P2 [67]. In addition, ABa and P2 align on the long axis (defining the anterior-posterior poles), while ABp and EMS align on the short axis (defining the dorsal-ventral poles) [68]. Therefore, from temporal and geometrical viewpoints, the derivatives of ABp and EMS are closer to each other, so a more powerful quantitative tool is needed to evaluate their development. At the molecular level, P2 is reported to induce polarization in ABp and EMS using MOM-2/Wnt signaling by direct contact between the cells [69].

Experimental setting
The C. elegans data were obtained according to the previous report of Tiraihi and Tiraihi [70]. Briefly, a C. elegans embryo (from the 4-cell to the 80-cell stages) was considered, the Cartesian coordinates (x,y,z) were estimated from ABp and EMS cell lineages using the images obtained from SIMI-Biocell [71] and Angler softwares [72]. The Cartesian coordinates were entered into a computer program to calculate the distances between the cells.
The distances at 30, 55, 82, 109 and 123 minute intervals (fixed intervals) were used at different scales and the data were entered into a computer program used to calculate the zero centro-axial skew-symmetrical (CSSM) and the basic square matrices (BSM).

Straight line
In the zero order straight line, two points represent the beginning and the end of the line. This line was divided into 48 unit lengths representing the steps (48 steps). The first order line was divided into two segments of equal length. At higher orders, each segment was divided into two equal parts, the box number increasing with the increase of order, leading to duplication in the number of the boxes and reduction in the step (see table 1). The SIPL method was applied to the straight line and the calculations were done as for the Koch curve (see below) except that the components of coordinates were taken from the straight line. The numbers of points at each order used in the study are presented in table 2.

Koch curve
The zero order Koch curve is a line comprising two points at the beginning and end, which are named "initiator points". In order to generate a higher order Koch curve, every line segment was divided into three equal segments. We named the first and the third segments "resting segments" and the second a "generating segment". The resting segments stay unmodified, and as the name implies, generation takes place in the generating segment. For each line, if we build an equilateral triangle on the generating segment, the two new lines are called "generated segments". As the final step in generating the next order, we remove the generating segment(s).
In order to generate a CSSM, we need a set of Cartesian coordinates as the input. For every Koch curve, we consider the end points of every line segment. Before the generation of a higher order curve, the current points satisfying this criterion are called "resting points". After the generation, the newly generated points that satisfy this criterion are called "generated points". The ends of each line segment at a certain order are called that order's principal points.
The algorithm for generating a CSSM from a set of points was described earlier. The initiator line has two principal points, while the 1 st and 2 nd order Koch curve have 5 and 17 points, respectively.
A computer program was developed to generate the two-dimensional Cartesian coordinates of the points (as described above) of a Koch curve. In this program, the Table 1 The box counting method data used in estimating the scale-invariant power law coefficient of the straight line.

Number of box (NB) Logarithm (NB)
Step ( The box number (NB) and the steps (h) as well as the logarithm of NB and the logarithm of the inverse of step are presented. These data are plotted in figure 1, and the coefficients of the regression line were estimated in order to calculate the SIPL coefficient. Table 2 The straight line orders and the related parameters used in calculating the scale-invariant power law coefficient using CSSM.
The program calculated the components of the Cartesian coordinates of the principal points for the five orders at different selected scales. These principal points were entered into another algorithm in order to generate the zero centro-axial skew-symmetrical matrix (CSSM) and construct the basic square matrix (BSM) according to the method described in the appendix. The program calculated the two-dimensional Cartesian coordinates (x,y) at the different orders. In the first order (see Figure 2), the program calculated 5 principal points as (1 × n) matrices, hence there were five elements for the x (x 1 , x 2 , x 3 , x 4 , x 5 ) and y (y 1 , y 2 , y 3 , y 4 , y 5 ) components. These (1 × n) matrices (principal row) were used to generate 5 × 5 CSSMs and 5 × 5 BSMs. In the same way, for the other orders, the principal points of Koch curve were calculated and the principal rows and CSSMs were generated and the BSMs were constructed. If there were identical elements in the principal row, then one element would be included in the principal row and the others omitted, otherwise the constructed BSM of the generated CSSM would result in a singular matrix with zero determinant.
For example, for the first order matrix, the program calculated 5 principal points forming two (1 × n) matrices with 5 elements for each coordinate. This was also done for all the scales in the first order.
The If(x a , x b ) are considered as the elements of the (1, n) matrix, then the first row of this matrix is 0 [1]. This was used to generate the CSSM for the x-components(x a ): where a stands for anti-symmetric.  The basic square matrix was constructed according to the algorithm presented in the appendix: The determinant of this matrix is -8. For the y coordinates, y a = 0 and y BSM = 0, the determinant of y BSM is zero. Another CSSM was generated from the combined determinants of x BSM and y BSM ; the input elements for the construction of this CSSM were (-8,0). It was translated to the point of origin to construct the CSSM; the principal row was [0,8]. The resulting matrix was: And the resulting basic square matrix was: The determinant of the basic square matrix (det(xy combined BSM )) was -512, this value represents the determinant of the BSM at Koch curve segment level. Five orders were used in the study (0 th , 1 st , 2 nd , 3 rd and 4 th ), where the 0 th order is the initiator of Koch curve (straight line). Seven scales (10 0 , 10 -1 , 10 -2 , 10 -3 , 10 -4 , 10 -5 and 10 -6 ) were used in the calculations. At each scale, the determinant of (xy combined BSM ) was calculated.
In the third sub-operation, for a given scale, the data from the different orders at each scale were plotted using a log-log plot, where the abscissa was the logarithm of the inverse value of step (h) (log(1/h)), while the ordinate was the logarithm of the absolute value of n det (BSM order ) log n det (BSM order ) . The log-log plot was fitted for linear regression and the b coefficients for the scales (b scale ) were subsequently used in estimating the SIPL of the Koch curve.
For the next operation (scale level), the logarithm of the scale (log(scale)) was plotted against the logarithm of the absolute b scale and another linear regression was calculated. The b coefficient(b SIPL ) was used in order to estimate the SIPL according to this equation: SIPL = 1-(D), where SIPL is the scale-invariant power law coefficient, and D is b SIPL .

Assessment of validity
The diffusion-limited aggregate method was used to estimate the fractal dimension of the growing embryo according to Moatamed et al. [73]. Briefly, 5 concentric circles with 5 μm increments were superimposed on the center of gravity of ABp-derived cells (ABp-dc) and EMS-derived cells (EMS-dc) at the 123 min. stage of development, and the nuclei of the ABp-derived cells were counted within each circle. The log of the nuclear number was plotted against the log of circle radius, and the slope of the regression line was used as the value of the fractal dimension. The same procedure was done on EMS-dc.

Programming languages
A computer program was written in the C++ language and a text file was generated containing the basic square matrix, which was copied into the command of the matrix of MATLAB ® software (http://www.mathworks.com: MathWorks, Inc, Natick, Massachusetts) and its determinant was calculated. Also, at each scale (10 0 , 10 -1 , 10 -2 , 10 -3 , 10 -4 , 10 -5 and 10 -6 ), the determinants of the basic square matrices were calculated for the following Koch curve orders (0, 1, 2, 3 and 4). The step for each order was also estimated and entered into the calculations.

Straight line
The results for the straight line using the box counting method are presented in detail in table 1 while Figure 3 presents Richardson's plot; the slope of the regression line equals zero and the SIPL coefficient is one. Table 2 presents the data and the calculations for the straight line using the CSSM. It shows the 4 orders and the number of points at each order, the b coefficients of the linear regression at each order with different scales, the logarithm of the absolute value of b coefficients and the scales used for calculating the second regression line in order to estimate the SIPL using the b coefficients (D).   10 -3 and 10 -4 , respectively. The abscissa represents the logarithm of the inverse of the steps plotted against log n det (BSM order ) (ordinate). In each plot, the 0 th , 1 st , 2 nd , 3 rd and 4 th orders were entered in order to calculate the linear regression. The logarithms of the absolute values of the b coefficients of the above plots were used for estimating the SIPL of the straight line, which were plotted (ordinate) against the logarithm of the scales (abscissa). The SIPL was calculated as follows: SIPL = 1-(D), where (D) is the b coefficient of the regression line.    Table 3 Koch curve orders and the related parameters used in calculating the scale-invariant power law coefficient using CSSM.

Developing embryo
The SIPL coefficients were 1.34176236 and 1.33912941 in ABp-dc and EMS-dc, respectively; the regression line of the plotted data is presented in Figure 8.

Discussion
The primary data used in the fractal analysis of the morphological studies are the lengths of the boundaries of the structure, which become more irregular with the increase of scale or magnification [74]. In this study, the components of the Cartesian coordinates were used as the primary data in the calculation of the SIPL coefficient. The development of an embryo from a fertilized egg is a dynamic process with space and time domains [75]. The spatio-temporal dynamics in developing embryos have studied by several investigators; for example, Wu et al. evaluated fetal development by plotting the fractal dimension of the brain surface at different stages of development against the time of brain development (in weeks) [76], and a similar approach was used by Schaffner and Ghesquiere in evaluating the complexity of type 1 astrotrocytes using the changes in the fractal dimensions during in vitro differentiation of the cells against the time of cell culture [77]. A different approach was taken to evaluating the changes in the value of the fractal dimension using it as a function of time [78]. While self-similarity is a property of fractals [79], cells in the developing C. elegans embryo migrate in a defined direction, becoming located in specific positions as they move. At each stage, the position of each cells is entirely different from the preceding and succeeding stages [70], so self-similarity does not hold during the development of the ABp and EMS sublineages of the C. elegans embryo. Fractal analysis was used in this study in order to confirm only the power law property of the SIPL, not self-similarity. The straight line and Koch curve The SIPL coefficient of the straight line using this method was equal to the fractal dimension using the box counting method, which is consistent with the analytical value of the fractal dimension. The standard error of this regression line was zero. The SIPL of the Koch curve was 1.25994438 with standard error and correlation coefficient equal to 0.1214845 and 0.9810513, respectively. The deviation from the analytical fractal dimension value was -0.15177%. The box counting method was reported to be the most suitable and appropriate of the methods discussed by Mandelbrot for estimating the fractal dimension [80] and it is the most popular one [81]. The criterion for selecting the method is based on the consistency of the fractal dimensions of the straight line and Koch curve. Mola et al. used the box counting method to estimate the fractal dimension of a Koch curve [82], and the results deviated by 2%. Applying the box counting method in calculating the fractal dimension using the fractal dimension calculator [83] resulted in a deviation from its analytical value of 6%. Cajueiro et al. used an isometric grid and the result deviated by 0.25% from the analytical fractal value [84]. Jiang et al. also used the box counting method and reported that 3, 4 and 5 recursive calculations of the Koch curve resulted in fractal values of 1.212, 1.226 and 1.255 (3.8, 2.6 and 0.4% deviation from the analytical value: 1.26) [85], respectively. Therefore, the consistency of the SIPL coefficient with the fractal dimension suggests the power law property of SIPL. Moreover, a few studies have used the Cartesian coordinates as the primary data set in calculating the fractal dimension [86,87], and reported it to be difficult [88], though the use of Cartesian coordinates in field of morphogenesis has been reported [89]. Also, the application of different methods of fractal dimensions in spatial classification has been considered useful [44], but in the development of the embryo, the spatio-temporal events are more essential for dynamic studies [70].

Diffusion-limited aggregation (DLA)
While diffusion-limited aggregation is used to study aggregation in clusters of particles in the time domain [79], morphogenesis in the developing embryo was studied using the DLA model, which confirmed the fractal property of the developing embryo [90]. DLA was suggested for quantifying the fractal dimension of blood vessel formation in the developing embryo [91] and other branching systems [92]. Therefore, it was used in this study in order to compare it with SIPL in the developing embryo. DLA was also used in the analysis of tumor vasculature [93,94], tumor growth [73,[95][96][97][98][99] and in vitro tumor spheroids [100]. The emergence of complex morphology in developing organisms was reported to be caused by DLA [101]. An in vitro study of embryonic retinal neurons showed a decrease in the number of neurite branches with an increase in viscosity of the medium, which was interpreted as a DLA mechanism [102]. Vilela et al. reported that DLA could be used in describing biological growth processes [103]. Therefore, it was used in this study in order to verify the feasibility of using SIPL in the developing embryo. DLA is a physical technique for describing the process of particle addition to a growing cluster of particles resulting in a power model for their number [79]. The fractal dimension calculated by DLA demonstrates the gradient of the diffusing substance toward the cluster [93]. Ryabov et al. reported that the gradient of these particles was proportional to the growth rate [104]. Therefore, the fractal dimension measured by the DLA method represents the growth rate, where the DLA value indicates that ABp grows more rapidly than EMS; this is also consistent with the SIPL results. Moreover, the value of DLA in ABp and EMS is different from the value estimated by SIPL, and the difference could be attributed to the feature of fractal dimension measured by DLA [105]. On the other hand, the large difference in values between ABp and EMS can be explained on the basis of the property of the DLA model. Peker et al. reported that two types of cluster formation resulted from DLA: the first underwent simultaneous aggregation-fragmentation processes or restructuring during growth; in the other type, growth depended on the immediate environment of the position of the new particles [96]. In the developing embryo, ABp-dc and EMS-dc can be considered as separate clusters, and the differentiating cells as particles. In previous studies, the results showed that EMS-DC had earlier regionalization, as the cells were tethered to neighboring cells by adhesion molecules [70], which is consistent with the second type of cluster formation. On the other hand, the cells (particles) of ABp-dc are separated into four groups, resulting in the spread of the cells in all directions (x, y, z), hence the value of the fractal dimension between 2 and 3 (which may indicate cell behavior similar to the first type of cluster formation) compared with that of EMS-dc (between 0 and 1), as the cells moved in one direction [70]. This analysis of ABp-dc is a consistent property of the fractal dimension of DLA, where the fractal dimension is a function of time, and the rate of aggregation is reduced as the cluster increases, as a result of spontaneous restructuring [106].

The developing embryo
A qualitative study was done on self-organization based on the expression of PAR-3, PAR-6 and PKC-3 at the anterior pole of the C. elegans ovum with PAR-1 serine/ threonine kinase and the PAR-2 proteins resulting in anterio-posterior polarization of the embryo [107]; also, there was self-organization at all biological hierarchy levels [1]. In fact, an in vitro study on gastrulation in isolated P1-descendent cells at the 26-cell stage of the C. elegans embryo (tracked in the cultured isolated cells using a 4D videomicroscope) showed that the cells gastrulate similarly to those of the intact embryo [108]. Also, a quantitative evaluation of the motion of the isolated cells using an in vitro setting showed that the onset of cellular motion was similar to that in the intact embryo (in vivo) and that the direction of the P1-descendent cells was also similar to that of the in vivo tracked cells [108]. This may indicate that patterning is required for cell population dynamics with the tendency of cells to associate with each other during gastrulation resulting in self-organization [49]. For example, in the vertebrate embryo, vascular patterning is essential for gastrulation movement [109]. On the other hand, in invertebrates, a quantitative study on video microscopy tracked cells at the early stages of C. elegans embryogenesis showed that their motion in the intact embryo was non-random in both EMS and ABp sublineages [70]. In this investigation, we showed that the SIPL coefficient of the EMS lineage is lower than that of the ABp, which is consistent with previous findings about the forward migration index (an index for chemotatic bias) [70]. Therefore, the lower value of the SIPL coefficient is not due to a decline in the complexity [26] or decomplexication of the sublineage [25], but shows that EMS-derived cells start to form an organized structure more rapidly as they tend to regionalize earlier than those of ABp [70]. Tabony documented that self-organization was not related to a single element such as the founder cells in the C. elegans embryo, but arose from the non-linear dynamics of all the elements collectively coupled to each other [110]. Schulze and Schierenberg revealed that embryonic cell lineages of low complexity formed a single sublineage or generated a single tissue type [111], as in the case of the E cell (from the EMS sublineage). Also, Quintana et al. reported that the isolated cells should adhere to each other in order to form an organized structure and suggested that a similar mechanism operated during the early stages of embryogenesis [12].
The results of this investigation are consistent with previous calculations of the diffusion coefficient, a physical parameter reported in a previous communication to be lower in the EMS lineage than in ABp during the early stages of embryogenesis in C. elegans [70]. Moskal and Payatakes reported that a decrease in the fractal dimension indicated a reduction in the diffusion coefficient, so the lower scaleinvariant power law coefficient in the EMS lineage may indicate an overall slower motion of EMS-derived than Abp-derived cells. A reduction in the fractal dimension indicates a reduction in the Brownian diffusion coefficient [112], which possesses a random walk property [113], so the results indicate that ABp cells moved more randomly than those of EMS.
The values of the SIPL coefficients in ABp and EMS are 1.34176236 and 1.33912941, respectively (small difference), but these values were taken from logarithmic regression and their anti-logarithms are 21.97 and 21.83, respectively, which are higher than the original values.

Conclusion
The study demonstrates that self-organization takes place during the early stages of embryogenesis, as confirmed by a scale-invariant power law method, calculated by using a centro-axial skew symmetrical matrix. The latter was generated from the components of the Cartesian coordinates. The SIPL coefficient results are consistent with DLA, where the fd calculated by DLA indicated that the cells in ABp-dc behaved similarly to type 1 cluster formation, while in EMS-dc they behaved similarly to that of type 2.

Appendix
The model used for evaluating the methodology has three points; each was assumed to be a center of gravity of a nucleus of a cell with two neighboring cells. The Cartesian coordinates of the cells were taken in two dimensions and the point of origin of the Cartesian coordinates was translated sequentially to the centers of gravity of all the cells. The x-and y-coordinates of the three cells were used for generating a(3 × 3) zero centro-axial skew-symmetric square matrix (CSSM). Another matrix, a basic square matrix (BSM), was constructed by multiplying the negative elements of the skew matrix with the negative operator (-1), while the positive elements were multiplied by the positive operator (+2). Then all the elements were multiplied by a scalar operator (+2) and the determinants of the matrices were calculated using Matlab software http://www.mathworks.com.

CSSM generation
There are two algorithms for generating CSSM. The first algorithm was accomplished according to the model shown in Figure 9.  Table 4 presents the arrangement of the x-components, while table 5 presents the y-components for the cells A, B and C in the three states. The array of numbers in table 4 can be presented in this skew matrix, which represents the x-component of the Cartesian coordinates of the centers of gravity of the nuclei in A, B and C cells: This is a zero centro-axial skew-symmetrical matrix (CSSM). The second algorithm for CSSM generation was done as described in the methods section. Initially, (1 × 3) matrices for the x-and the y-components [(0,3,11) and (0,7,2), respectively] were used for generating the skew matrices. The CSSM were generated: Where i and j are indices for rows and columns, respectively. In Euclidean space, the formula can be generalized into n-dimensional space where n = 1, 2, 3, 4.....N, and N is a finite number.
x n i,j = x n 1−1,j − x n i−1,i Table 4 The array of the x-coordinates of the three states of translation for the three cells. According to Hadley's [114] notation, the skew symmetric matrices can be rewritten as follows: The determinants of the basic square matrices (and det(y BSM )) were 12672 and 3360, respectively.
A 2 × 2 CSSM can be generated from the values of det(y BSM ) by arranging them into a (1 × 2) matrix in the form of [det(x BSM ), det(y BSM )]. This matrix will then be input to the CSSM generation routine and a BSM will be constructed from the CSSM.
For example, for the above values, the input matrix to the CSSM generation process (as described above) will be [12672, 3360] and the CSSM will be (see table 6): xy CSSM combined = 0 −9312 9312 0 and the BSM Matrix will be: xy CSSM combined = 0 18624 37248 0 Table 5 The array of the y-coordinates of the three states of translation for the three cells.  Table 6 The (2 × 2) array generated by translation of two points (det(x BSM ) and (det (y BSM )) on the real axis to the origin of real axis.