Form-Finding of Spine Inspired Biotensegrity Model

Featured Application: The advantageous characteristics of ﬂexibility and versatility of movement in spine biotensegrity models can be further studied for potential application in deployable structures and ﬂexible arm in the robotic industry. Abstract: This paper presents a study on form-ﬁnding of four-stage class one self-equilibrated spine biotensegrity models. Advantageous features such as slenderness and natural curvature of the human spine, as well as the stabilizing network that consists of the spinal column and muscles, were modeled and incorporated in the mathematical formulation of the spine biotensegrity models. Form-ﬁnding analysis, which involved determination of independent self-equilibrium stress modes using generalized inverse and their linear combination, was carried out. Form-ﬁnding strategy for searching the self-equilibrated models was studied through two approaches: application of various combinations of (1) twist angles and (2) nodal coordinates. A total of three conﬁgurations of the spine biotensegrity models with di ﬀ erent sizes of triangular cell were successfully established for the ﬁrst time in this study. All members in the spine biotensegrity models satisﬁed the assumption of linear elastic material behavior. With the established spine biotensegrity model, the advantageous characteristics of ﬂexibility and versatility of movement can be further studied for potential application in deployable structures and ﬂexible arm in the robotic industry.


Introduction
The principle of biotensegrity was introduced by Levin [1] and Ingber, et al. [2] in the 1980s, as an idea of applying the concept of tensegrity to represent the interaction of forces in all hierarchical biological systems. Biotensegrity has been found to demonstrate inherent mechanistic properties that are in good agreement with experimental data of in vivo [3,4]. It also possesses the following four excellent characteristics: efficiency (i.e., geodesic form), self-stabilizing, multi-modularity and multi-functional (i.e., mechanotransduction). To date, the overwhelming body of evidence clearly indicates the applicability of the principle of biotensegrity in anatomy and physiology from macro-to nano-scale biology systems [5,6]. Biotensegrity as a model emulated from the forms and functions of sophisticated and sustainable biological systems has been suggested for several potential applications especially in medical technology, health care practices [5] and even in the construction industry in the form of a biotensegrity robot [7]. This study is aimed at form-finding of possible biotensegrity models inspired by the form of the human spine. Advantageous features of a human spine such as slenderness and natural curvature of the geometry, as well as the stabilizing network consisting of spinal column and muscles were incorporated in the configuration of the spine biotensegrity model. This study was carried out as a first step towards the next phase of study on the shape changeability of the spine biotensegrity models that can be beneficial for applications such as automated robotic tools (i.e., for inspection and automation in manufacturing) or deployable structures (Figure 1). Shape change of regular tensegrity models studied by the authors can be found in [8]. The form-finding analysis employed Moore-Penrose generalized inverse in the solution of rank-deficient system of equilibrium equations. An algorithm incorporating iterative process to determine the self-equilibrated configuration of spine biotensegrity model through the modification of nodal coordinates was proposed. The human spine is one of the unique models in the biological system with the ability to sustain high loads and at the same time demonstrate a superior combination of multidirectional movements. Despite being spongy, the tapered spinal vertebrae body efficiently maintains body weight and withstands external impacts. In addition to the special architecture with sagittal curvature of spinal column, the spine can avoid distress (i.e., exhaustion) by re-establishing new equilibrium positions. The dexterity in avoidance of distress in the system is mainly due to the harmonization of the passive (i.e., spinal column), active (i.e., surrounding muscles and tendons) and control (i.e., neural centers) subsystems [9]. It is noted that the spine always maintains the stability of the system even after disturbance by a self-structural reorganization. Based on the above descriptions, a spine has a distinctive characteristic of flexibility and versatility in accommodating different movements. Coupling the advantage of tensegrity (i.e., lightweight and flexibility) and the versatility of a human spine, a novel biotensegrity model can be realized. In addition, the principle of biotensegrity is highly suitable for explaining the mechanism of compressive passive subsystem (i.e, spinal column) and tensional active subsystem (i.e., surrounding muscles and tendons) in the spine. According to Levin, a spinal system re-structures to form many new postures in order to sustain the stability of the system once the tensional forces in the network have been disturbed [10]. Levin also recommended a geodesic, close packing and hierarchical spine model using tensegrity mast. More information about spinal tensegrity mast models can be found in [11].
Configurations of biotensegrity have been proposed through form-finding process by solving complex geometrical inputs and forces [12]. The first physical biotensegrity model that mimics living cells, namely cellular tensegrity model, was built with the spherical central nucleus [13]. Thereafter, mathematical modeling with various cellular models using regular tensegrities were studied [12]. Many have claimed excellent prediction of cellular mechanistic properties through cellular tensegrity models [5,6,14]. However, the major limitation of these models lies in the use of simple configurations to represent the complex structure in nature [11]. Despite the fact that the basic mathematical formulation to express the mechanical models for biotensegrity has achieved a certain maturity, many of the investigations focused extensively on cellular tensegrity models [6,12]. In this study, a form-finding of biotensegrity models to specifically represent the complex living organism at a higher hierarchy is addressed. There are many established form-finding strategies on tensegrity structures, which have been classified as kinematical and statical methods [15]. A force density-based analytical form-finding method was applied for regular tensegrity structures by group elements [16]. Form-finding method combining symmetry-based qualitative analysis with particle swarm optimization for symmetric tensegrities and prestressed cable-strut structures was proposed [17]. A method using strut finite element and unconstrained nonlinear programming for form-finding of hyper-elastic tensegrity structures was studied [18]. A nonlinear programming approach for form-finding of tensegrity structures by using fictitious materials properties was presented [19]. An optimization method to generate nodal coordinates of super-stable tensegrity structures with known connectivity data and force densities was proposed [20]. Form-finding analysis of an asymmetric large deployable antenna using a force density-based optimization algorithm was carried out [21]. A form-finding method whereby a fully automatic labeling of elements through a double-loop genetic algorithm was used [22]. A modified dynamic relaxation algorithm for the form-finding of tensegrity robots in task space was proposed [23]. Many of the investigations mentioned above were aimed at finding new configurations of tensegrity subjected to different constraints or optimization criteria. This paper presents form-finding strategies for spine biotensegrity models.
The remainder of the paper is organized as follows: Section 2 defines the geometry of the spine biotensegrity models that are based on key parameters from human spine vertebrae and natural sagittal curvature. Section 3 presents the search strategies for form-finding of self-equilibrated configuration of the spine biotensegrity models adopted in this study. Section 4 discusses form-finding results of spine biotensegrity models using the algorithm proposed. Lastly, the conclusions are presented in Section 5.

Geometrical Input of Spine
The geometrical input data for the generation of spine biotensegrity models were gathered from the anatomical parameters of human spines by Busscher, et al. [24] where the human spine specimens were collected from six male cadavers with age and body height at death ranging from 55-84 year-old and 175-192 cm, respectively (mean of 72 year-old and 182cm, respectively). Figure 2 shows typical vertebrae in superior and lateral views whereas Figure 3 shows a typical human spine in sagittal and posterior views. Sagittal and posterior views are medical terms that refer to the appearance of the spine when viewed from the side and behind, respectively. Figure 3 shows a total of twenty-two vertebrae from the cervical (vertebrae C3,C7), thoracic (vertebrae T1T12) and lumbar regions (vertebrae L1-L5) that were considered in this study. Anatomical parameters such as the vertebral end-plate width denoted as EPW (mean value of upper and lower vertebral end-plate width denoted as UEPW and LEPW, respectively), vertebral central body height denoted as VBHC and intervertebral disc height denoted as IDH from the source by [24] were considered in the generation of the spine biotensegrity models.

Natural Sagittal Curvature
The human spine ( Figure 3) has a natural curvature in sagittal view with: (i) lordosis (concave curve) at cervical and lumbar regions and (ii) kyphosis (convex curve) at thoracic and sacral regions. According to Busscher, et al. [24], cervical lordosis of 20.1 • , thoracic kyphosis of 34.5 • and lumbar lordosis of 29.2 • were measured as the angles between the upper end-plate of C3 and the lower end-plate of C7, the upper end-plate of T1 and the lower end-plate of T12 and the upper end-plate of L1 and the lower end-plate of L5, respectively. These angles are reasonably accepted and supported by others. For instance, the fact that thoracic kyphosis lies between 20-45 • has been suggested by Mejia, et al. [25]. Been and Kalichman [26] have reported a normal range of 30-80 • for lumbar lordosis. The apex of the curvature is defined as the farthest point of a vertebra from the C7 plump line. The C7 plumb line is established by extending a vertical line from the center of C7 vertebral body to the posterior corner (at the back) of sacrum S1 [28]. The normal weight of the human body generally acts through the C7 plumb line. Nevertheless, there are insufficient anatomical data about sacrum S1 in the study by Busscher, et al. [24]. Hence, an assumption that the C7 plumb line rests on the right end of the vertebrae L5 adjoining to S1 was made in this study. Figure 4 shows the established natural sagittal line for the spine biotensegrity models. It is noted that there is an apex each within cervical lordosis, thoracic kyphosis and lumbar lordosis in a human spine. The decision on the location of these apexes is based on several pieces of literature. The apex of lumbar lordosis, AP1, is assumed at the intervertebral disc between L3-L4 whereas the apex of thoracic kyphosis, AP2, is assumed at vertebrae T7 [29]. The apex of cervical lordosis, AP3, is assumed at C5 according to Knox and Lonner [30], who suggested the lordotic climax in between C4-C6. One of the most common methods, namely the Cobb method, has been adopted in establishing the curvature of the spine biotensegrity model. A review on this method can be found in [31].

Search Strategies for Spine Biotensegrity Models
This section presents the key steps in the search strategies for spine biotensegrity models: (1) the generation of geometrical model, and (2) the formulation of mathematical equations.

Geometry of the Spine Biotensegrity Model
The geometry of the spine biotensegrity model was defined by their nodal information, which is described in the basic cell such as the triangular cell and element connectivity.

Definition of Basic Cell
The triangular cell has been chosen to mimic the human spine, due to its simplicity of application. Figure 5 shows the definition of a typical triangular cell and triangular surfaces of the spine biotensegrity models. A total of three variations of spine biotensegrity models identified as SB1, SB2, and SB3, with different sizes of triangular cells have been proposed in this study. The geometry of a triangular cell is presented in the global Cartesian coordinate system O o -X o Y o Z o . The position vector of nodal coordinates is denoted as X. The origins of the local Cartesian system O i -1 -X i -1 Y i -1 Z i -1 and O i -X i Y i Z i are positioned at the centroid of lower (X 6i-5 X 6i-4 X 6i-3 ) and upper triangular surfaces (X 6i-2 X 6i-1 X 6i ) of ith (i = 1, 2, . . . , n) triangular cell, respectively. The nodal coordinates in plane XY of the spine biotensegrity models are determined based on triangle altitude, d and half of the vertex angle, θ of the lower and upper triangular surfaces in the local Cartesian system. The tapered section properties of a human spine are incorporated into the spine biotensegrity models (the width of vertebra gradually increases from cervical vertebrae to lumbar vertebrae). Reference is made to the EPW of vertebrae to characterize the widths of the spine biotensegrity models and these widths are adopted as the triangle altitudes for the triangular surfaces in cells.
It is noted that all spine biotensegrity models have the same magnitude of triangle altitude (d) (see Figure 5). SB1, SB2, and SB3 are different in terms of half of the triangle vertex angles (θ) with values of 20 • , 30 • , and 40 • , respectively. While accumulation of both VBHC and IDH represents the height of the spine biotensegrity model, accumulation of IDH represents the saddle height between two independent cells.

Definition of Basic Cell
Tensegrity are classified into class one, two and k depending on the numbers of strut at a node [32]. In this study, class one tensegrity mast is chosen for the spine biotensegrity models to mimic the stacking up of discontinuous compressive vertebral bodies of a human spine. Specifically, class one tensegrity mast consists of struts (i.e., vertebrae) that are connected solely by cables (i.e., muscles, intervertebral disc) and offers no contact in-between the struts, such as a human spine (i.e., vertebrae are cushioned by intervertebral disc and surrounded by muscles). Figure 6 shows the developed element connectivity chart inspired by Pugh's illustration of the diamond pattern. The N-stage class one spine biotensegrity models with N triangular cells are assembled according to the chart. The element connectivity at the second and third stage may be repeated in order to stack up to n triangular cells. The lower and upper triangular surfaces of ith (i = 1, 2, . . . , n) triangular cell are denoted as L i B and L i T, respectively. The triangular cell that is stacked up on the current ith triangular cell is denoted as the (i + 1)th cell. All the elements are classified into specific groups. The struts in every stage are represented by bold lines. The saddle cables (red color) connecting the interface between two triangular cells are represented by irregular dash-dotted lines. Horizontal cables in regular dashed lines join the nodes only at the base and the top of the structures. There are three categories of diagonal cables; namely diagonal cables 1 (cyan color) connecting nodes L i B and L i + 1 B, diagonal cables 2 (pink color) connecting nodes L i T and L i + 1 T, and diagonal cables 3 (blue color) connecting nodes L i B and L i T. On top of that, reinforcing cables in zig-zag lines at the first stage are also included.

Assemblage of Spine Biotensegrity Models
As mentioned earlier, a total of three variations of spine biotensegrity model, namely SB1, SB2 and SB3 are developed. These three four-stage spine biotensegrity models are proposed based on the nature of three apexes in a human spine. Figure 7 shows the assemblage process of a spine biotensegrity model with four triangular cells. Firstly, major positions of the eight triangular surfaces for the four triangular cells are decided. The triangular surfaces are placed as follows: one at the base, two at each of the apex and one at the top. The altitudes of the triangular surfaces at each stage are based on the anatomical parameter EPW at the respective height of the human spine specimens. The altitude of triangular surfaces decreases along with the height until the top of the spine biotensegrity model to change accordingly.the taper-form and slenderness of a human spine. These triangular surfaces are then aligned perpendicular to and have their centroids coinciding with the established curvature line as shown in Figure 4. Note that, nodes X 6i-5 and X 6i-2 of ith triangular cell ( Figure 5) are positioned on global Y axis (x = 0), prior to the twisting process (Figure 7a). Next, the triangular surfaces are rotated so that all the triangular cells have an initial twist angle in the range of 10 • to 50 • in clockwise and counterclockwise for odd-and even-number of ith triangular cell, respectively (Figure 7b). Furthermore, the rotational angle between the triangular surfaces forming the interface of two triangular cells is set at 60 • with the same rotational direction as the twisting direction of the lower cell.
Lastly, the triangular surfaces are removed except for their nodes at the established positions. The generation of the spine biotensegrity model is complete after the placement of all the struts and cables according to the developed element connectivity pattern as shown in Figure 6 (Figure 7c). Each spine biotensegrity model consists of sixty-nine elements (i.e., twelve struts and fifty-seven cables).
The following section presents the formulation of the equations for the search of the self-equilibrated state of a spine biotensegrity model.

Equilibrium Equations of Spine Biotensegrity Model
This study is aimed at establishing the feasible self-equilibrium configuration for four-stage class one spine biotensegrity models by the means of form-finding analysis.

Basic Assumptions
Form-finding of four-stage class one spine biotensegrity models that are stabilized by prestresses at the self-equilibrated state is carried out. The following assumptions are made for form-finding analysis: • Self-weight of struts and cables are neglected due to the reason that the stress resulted from self-weight is insignificantly small compared to the level of pre-stress [19].

•
There is no external force acting on the models. It is necessary to ensure the self-equilibrium position as shown under the topology of biotensegrity model (Figure 7) is rational. Against any disturbance, the self-equilibrium and stability of the obtained topology should be retained. Analysis to study the capability of the established spine biotensegrity to undergo change in shape is not the scope of the study.

Equilibrium Equations: Formulation
As the spine biotensegrity model consists of struts and cables that are pin and connected, the set of equilibrium equations adopted in this study is derived based on the basic formulations of pin-jointed structures. Such formulation for tensegrity structures can also be found in [32,33].
The spine biotensegrity model in this study consists of m elements, n o nodes, n c supports and n u = 3(n o -n c ) unconstrained degree of freedoms. Consider an element k of the spine biotensegrity model connecting node i and j (i < j) in a Cartesian coordinate system O-XYZ as shown in Figure 8.  The length of element k is evaluated using the following expression: The vector of directional cosines λ of an element k is defined as follows: Figure 8 also shows the axial force n of element k for the spine biotensegrity model, and the components of the nodal force at i and j, f i and f j , respectively (where f i =f j ), which are in equilibrium with n. Considering equilibriums at all nodes, the following equation can be obtained: where B is a n u × m matrix consisting of the vector of directional cosines λ of all elements in the spine biotensegrity model with respect to x, y, z axes, F is a vector of nodal forces with size n u and n is a vector of axial forces for all elements with size m. By making use of the Moore-Penrose generalized inverse, the necessary and sufficient condition for the existence of solution for Equation (3) can be expressed using the following equation: which can be rewritten in the following form where I n u is an identity matrix of order n u and B + is Moore-Penrose generalized inverse of matrix B. Singular value decomposition method is used in this study to compute the Moore-Penrose generalized inverse of the matrix B with rank deficiency (i.e., rank < n u ).
For the case of self-equilibrated state without external forces, F = 0. Under such condition, the solution for the vector of axial forces n for the spine biotensegrity model is given by the following equation: where β is a vector of arbitrary coefficient of size m.
Denoting matrix [I m − B + B] in Equation (6) as matrix G of size m × m and applying eigenvector basis decomposition method, the real symmetric matrix G can be expressed as follows: where Λ is a diagonal matrix whose diagonal elements Λi (i = 1, 2, . . . , m) correspond to eigenvalues of matrix G; Φ is an orthogonal matrix whose ith column Φ i (i = 1, 2, . . . , m) corresponds to eigenvector associated with eigenvalue Λi. Since matrix G is an orthogonal projection matrix, the eigenvalues Λ i returns the value of either 0 or 1. The eigenvectors Φ i correspond to vectors of independent self-equilibrium stress modes for biotensegrity models considered in this study. It is noted that the number of vectors of independent self-equilibrium stress modes is equal to the number of non-zero eigenvalues of matrix G or rank [G]. Substituting Equation (7) into Equation (6) and then multiplying both sides of the resulting equation with Φ i −1 as follows: the following expression for vector of axial forces n can be obtained: where β * = Φ T β.
Assuming that the spine biotensegrity model has q number of independent self-equilibrium stress modes, Equation (9) can be further rewritten as where β i * is ith component of vector β*, which corresponds to vector of coefficients for the linearly independent self-equilibrium stress modes. Equation (10) shows that the solution of vector of axial force n is dependent on the choice of the combinations of coefficients β i *, which might be more than one. In this study, only one of the possible combinations of coefficients β i * was used to determine the vector of axial forces n.

Equilibrium Equations: Solutions
In this study, the spine biotensegrity model is found to have three independent self-equilibrium stress modes (q = 3). In order to find the possible linear combination of the three self-equilibrium stress modes, an iterative procedure is adopted in this study. Two main steps are involved in the iterative steps. The first one is the determination of trial set of coefficients β i * for an assumed configuration of spine biotensegrity model. The other is the modification of the assumed configuration until one consistent self-equilibrated configuration is found. For the determination of the trial set of β i *, an unconstrained minimization problem as shown in Equation (11) is adopted where the conjugate gradient method is used: Equation (11) is an objective function to prevent excessive axial forces. The set of β i * (i = 1,2,3) determined by solving Equation (11) is used as the trial values to determine the possible solution for n (see Equation (10)), which satisfies the inequality constraints given in Equations (12) and (13) using an iterative process: where n c , E c and A c are axial forces, yield stress and cross sectional area for cable elements, respectively; n s , E s , I s , L s , σ s and A s are axial forces, Young modulus, moment inertia of section, length, yield stress and cross sectional area for strut elements, respectively. In this study, A c = 3.14 mm 2 , A s = 50.30 mm 2 , E s = 200 GPa, and I s = 201 mm 4 . As shown in Equations (12) and (13), in addition to checking against yielding of cable elements, compressive forces of strut elements are checked against Euler's buckling load. The applied inequality constraints ensure that the axial forces for all the cable and strut elements are within their upper bound axial force limits and satisfy the elastic condition of the material properties.

Algorithm for Form-Finding of Spine Biotensegrity Model
The flowchart in Figure 9 summarizes the strategies to search for a possible self-equilibrated state of the spine biotensegrity model. Algorithm 1 begins with the preparation of geometrical input of the spine biotensegrity model in the form of four-stage class one tensegrity mast with triangular cells. A spine biotensegrity model is generated with the application of initial twist angles to each triangular cell. The angles may be considered as the anatomical parameters. Vector of axial forces n satisfying self-equilibrium condition, as shown in Equation (10), is established. The eigenvector basis decomposition method is adopted to determine the independent self-equilibrium stress modes.
Algorithm 2 is the subset of Algorithm 1 and provides an initial estimate for coefficients β i * through the minimization problem as shown in Equation (11). It is noted that the combination of coefficients β i * is searched through an iterative form-finding process to satisfy the Equations (12) and (13). After obtaining the combination of the coefficients β i *, the substitution of the coefficients β i * into Equation (10), gives the axial forces n in the self-equilibrium state.
The form-finding strategy for searching the self-equilibrated models is studied through two approaches. Under approach one (Figure 10a), different combinations of twist angles are applied to a group of basic cells of a trial spine biotensegrity model. Under approach two (Figure 10b), unconstrained nodal coordinates of a trial spine biotensegrity model are modified through a relatively small magnitude (2% of the height of spine biotensegrity model in this study). Design valuables are coefficients β i * and nodal coordinates. As the objective function and the constraint conditions cannot be explicitly expressed by design variables, the procedure of form-finding implies an iteration. Both approaches involve only minor modification of nodal coordinates due to the priority of searching for models with close resemblance to the geometry of the human spine.

Search Results for Spine Biotensegrity Models
The process of finding the possible self-equilibrated state of spine biotensegrity models using the iterative approach (Figure 9) is presented. Specifically, the results of the processes through approaches one and two whereby number of slackened cables in compression (slackened cables) are continuously reduced and eventually eliminated in order to achieve self-equilibrated configuration of spine biotensegrity model are presented. Due to the similar process, only spine biotensegrity model SB1 was chosen to demonstrate the form-finding approaches.

Analysis Results Using Approach One
Under approach one, a random choice of each stage's twist angle ranging from 10 • to 60 • has been selected for the form-finding analysis for the spine biotensegrity models. The symbol α i is defined by the twist angles for the i-th stage's triangular cells. Arrangement of the spine biotensegrity models with α i ≥ 60 • is found to be not feasible due to the problem of closely spaced struts. Figure 11 shows the form-finding analysis with the application of the same twist angle 10 • to 50 • for all the triangular cells in the spine biotensegrity models (i.e., cases 10101010, 20202020, 30303030, 40404040, and 50505050). The case 10101010 denotes the model as α 1 =α 2 = α 3 = α 4 = 10 • . The cables in compression are slacked and should be removed. The slackened cables of SB1 are indicated by the red line ( Figure 11). From the form-finding analysis, it is found that SB1 with α 1 = α 2 = α 3 = α 4 = 50 • yields the least degree of slackness (9 slackened cables). Larger twist angles tend to reduce the number of slacked cables, but at α 1 = α 2 = α 3 = α 4 = 20 • , the number of slacked cables is maximum. For 10 • ≤ α i ≤ 40 • , slackened cables are found to occur throughout the model. For α i = 50 • , more slackened cables are found to occur at the first and third stages. The information about slackened cables at each trial case is an important input in the iterative analysis. Figure 11. Effect of same twist angles on numbers of slackened cables for SB1. Figure 12 shows the analysis with the application of various combinations of twist angles for triangular cells within the model (i.e., cases 10505050, 20505050, 30505050, and 40505050). It is found that cases where the α 1 for the triangular cell at stage 1 is less than 50 • lead to a lesser number of slackened cables in comparison with other combinations of α 2 , α 3 , α 4 (Figure 12a). A similar effect of such combinations of twist angles is found in SB2 and SB3. When α 1 changes while the other α i are held at 50 • , the numbers of slackened cables ranged from three to nine for SB1, SB2, and SB3 are found (Figure 12b). In the model SB1 and SB3, the numbers of slackened cables are maximum at α 1 = 20 • . In the case of α 1 = 10 o in SB1 and α 1 = 40 • in SB3, the number of slackened cables is minimum. On the other hand, SB2 is not sensitive to α 1 . In SB1, as α 2 and α 3 are assigned larger values, the number of slackened cables decreases. Figure 13 illustrates the numbers and position of slackened cables for parameter α 1 = 10 • −40 • and α 2 = α 3 = α 4 = 50 • in SB1. At α 1 = 10 • , number of slackened cables decreases to five and is uniformly located throughout the stages.

Iterative Results Using Approach Two
Self-equilibrated state of spine biotensegrity models cannot be found by approach one. The case 30505050 (selected based on the form-finding results from SB1, SB2 and SB3) is then employed in approach two in the subsequent process of searching for the self-equilibrated configuration. Under approach two, the modified nodes are 18 nodes excluding the base and top nodes (i.e, node 4 to node 21). The controllable range of the nodal coordinates is less than 15 mm in the y direction. Few cases (i.e., cases A1 to H1) that are close to the self-equilibrated state are chosen to demonstrate the feasibility of the form-finding approach. Figure 14 shows the numbers and positions of slackened cables in SB1. It is seen that SB1 still preserves the spine anatomical parameters, especially the natural curvature of the spine. It is also important to note that the position of slackened cables changes even though the modifications of nodal coordinates are small.  Table 1 presents the axial forces of the potential slackened cables of SB1. As can be seen in Table 1, although the slackened cables 40, 50, 59, and 68 may change into tension state easily, the remaining cables (particularly cable 51 and 65) remain slackened in most of the cases. In the case of B1, it is apparent that the case with the modification of node 5 instead of node 6 reduces the degree of slackness in most of the slackened cable (especially in cables 40, 42, 57, 65, and 67). Additionally, the modification of node 19 in case F1 clearly shows the effect of further reducing the degree of slackness in the slackened cables 51, 59, 65, and 66 compared to case E1. However, such modification in case F1 also leads to an increase in the degree of slackness in slackened cables 42 and 67. This situation is improved by excluding node 13 in case G1. Since node 19 showed its significance in changing slackened cables 51, 59, 65, and 66 to tension state, further modification of coordinates of node 19 was carried out. This led to case H1, where a self-equilibrated SB1 is successfully determined. Similar form-finding strategies are applied to SB2 and SB3.

Spine Biotensegrity Models
Three variations of spine biotensegrity models (SB1, SB2, and SB3) are successfully established by using the proposed form-finding strategies. Figure 15 illustrates the configuration and nodal numbers. The larger nodes represent the modification through approach two. In any model, nodes of all cells excluding the nodes at top and base are considered for modification during approach two. Table A1 in Appendix A shows the nodal coordinates for the spine biotensegrity models.   (12) and (13). It is emphasized that the solution obtained for the self-equilibrated configuration for SB1, SB2, and SB3 represents one possible configuration with close resemblance to the geometry of the human spine. The self-equilibrated configuration could be further refined through solution of constrained optimization problem using maximum stiffness, minimum weight or other suitable criteria as objective functions.
In the process of determination of coefficients β i *, the combinations of coefficients β 1 *, β 2 * and β 3 * that lead to tension in struts or compression in cables are rejected. Non-acceptance of the trial values of coefficients β i *(i = 1,2,3) obtained from solution of problem defined by Equation (11) is always due to non-satisfaction of Equation (12) (the existence of slackened cables) within the trial configuration of spine biotensegrity models.
The proposed form-finding strategy using approach one and two yields updated trial configurations of spine biotensegrity models to reduce slackened cables and to maintain closely the geometrical characteristics of the spine. Approach two eliminates inconsistent cables and further refines the configuration. It is found that, although modification of nodal coordinates in some cases can reduce slackened cables, increment in the degree of slackness in other slackened cables and the generation of new slackened cables are observed. Therefore, a good understanding of the relationship between the modified nodes and the change in axial forces in a spine biotensegrity model is essential in searching for a self-equilibrated model.

Conclusions
This paper presents a trial-and-error iterative process for form-finding of three variations of biotensegrity models inspired by the form of human spine. The developed element connectivity chart inspired by Pugh's diamond pattern is able to generate N-stage class one biotensegrity structures. The proposed form-finding strategy consists of two approaches. In the first approach, the twist angles of the cells in each stage, which simulates the geometrical characteristics of the spine, are varied to reduce slackened cables. In the second approach for further refinement, nodal coordinates are varied to reduce and eliminate slackened cables. Through form-finding analysis, the characteristics of the biotensegrity models are found. Twist angle, except for the lowest cell, are appropriate at 50 o . The slacking of cables is very sensitive to twist angles near the thoracic and lumbar region. The sensitivity to twist angle of the lowest cell is dependent on the shape of the cell. This study has successfully established for the first time the possible self-equilibrated configurations of biotensegrity models inspired by the human spine. Future works such as shape change analysis, building and controlling real biotensegrity models and conducting experiments considering external loading and operational challenges as well as formulation of self-collision and external collision of the spine biotensegrity models should be carried out. The outcome of such future works could then lead to the possible application of spine biotensegrity in area of robotic tools and deployable structures.

Nomenclature
A c Cross sectional area of cable A s Cross sectional area of strut B Matrix of directional cosines B + Moore-Penrose generalized inverse of matrix B β* Vector of arbitrary coefficients for linearly independent self-equilibrium stress modes D Triangle altitude E s Young modulus of strut F Vector of nodal forces Φ Eigenvector associated with the eigenvalue Λ G Matrix for existence of solution for Equation (2) (see Equation (7)) I s Moment inertia of section L s Length of strut L T Length of cable or strut element Λ Eigenvalues of matrix G λ Vector of directional cosines n Vector of axial forces n c Axial forces of cable n s Axial forces of strut θ Half of the vertex angle σ c Yield stress of cable σ s Yield stress of strut − x Vector of nodal coordinates Appendix A