This paper presents a numerical strategy for the shape change analysis of spine biotensegrity models in multi-directional modes. The formulation of incremental equilibrium equations and optimization problem for shape change analysis via the forced elongation of cables to achieve the target coordinates of the monitored nodes of spine biotensegrity models are presented. The distance between the monitored nodes and the target coordinates is chosen as the objective function which is minimized subject to inequality constraints on member axial forces and cable forced elongation. Three spine biotensegrity models were analyzed to validate the effectiveness of the proposed method. The deformation characteristics of the Class-1 four-stage biotensegrity models mimicking the natural curvature of the human spine were investigated. A highly successful rate in achieving the target coordinates was observed in a total of 258 analysis cases, with percentages of 99.9%, 99.9% and 98.9% for shape change analysis involving uni-, bi- and tri-directional modes, respectively. The results show that the spine biotensegrity models have more flexibility in undergoing bending in comparison with axial deformation. With the established shape change strategy, the flexibility and versatility of the movement of spine biotensegrity models can be further studied for potential application in the shape change control of deployable structures together with the use of IoT.