Modeling approaches have not been employed in biomedical research related to spinal cord injury. However, in other areas of biomedicine, mathematical and computational modeling has become a standard approach over the last few decades.
Various embodiments disclosed herein relate to an algorithm for computational modeling of stem-cell-driven tissue growth of adult spinal cord tissue. Simulations based on this model can be used for parameter testing and making predictions about the growth dynamics of biological spinal cord tissue. Such predictions include alterations in the growth dynamics and the final properties of the tissue induced by experimental or medical intervention. The computational modeling has significant potential for making the development of biomedical therapies for the treatment of spinal cord injury more rational and efficient, thereby streamlining this process.
In accordance with one or more embodiments, a computer-implemented method is disclosed for predicting stem-cell-driven growth of biological spinal cord tissue under given conditions. The method includes the steps of: (a) electronically accessing a computational model from a computer storage for simulating stem-cell driven spinal cord tissue growth, the model based on a cellular automata (CA) framework comprising a plurality of lattice sites, wherein each lattice site is empty or contains a cell having one of four state values: stem cell, progenitor cell, differentiated cell, or dead cell; the model including a plurality of interaction rules for predicting the state value of each lattice site at a next time iteration based on the state values of cells in neighboring lattice sites at a current time iteration; (b) 65238398.1 electronically receiving input data specifying an initial lattice composition and lattice boundary conditions; (c) electronically running a simulation on the computational model based on the input data for a plurality of time iterations to predict development of the biological spinal cord tissue over the time iterations; and (d) electronically outputting data on the development of the biological spinal cord tissue.
Growth of spinal cord tissue is a complex biological phenomenon. To gain a deeper understanding of the theory of this phenomenon and to be able to make predictions about the outcome of experimental and biomedical interventions, we have developed a computational model that is able to simulate, with high precision, the outcome of the underlying biological processes. This tool has significant potential to make the development of therapeutic strategies (including the design of drugs) aimed at curing spinal cord injury more rational and cost effective.
The processes disclosed herein complement existing approaches based on in vivo or in vitro testing of the effect of experimental or pharmacological interventions to improve spinal cord regrowth after injury. Due to the predictability power of our invention, the design of such tests can be made more rational, and thus more cost effective. The modeling and simulation software even allows one to make predictions about the likely outcome of experiments that are difficult to perform due to resource constraints. Modifications of this technology can be applied to drug discovery related to tumor growth.
In the following, we describe the underlying mathematical model and outline the algorithm, which has been implemented in form of a MATLAB script. We also demonstrate the feasibility of this modeling approach by showing the results of a simulation based on the MATLAB script.
Our model employs the cellular automata (CA) framework in which time, space, and state are treated as discrete variables. In this mathematical framework, a finite number of agents equipped with states interact with each other on a uniform lattice. Each agent occupies one lattice space and its state at the current time iteration is determined by interaction rules dependent on state values of agents within its neighborhood prior to the current time iteration. During each time iteration, the state of an agent can take one value out of a finite set. Therefore, not only time and space but also the states of agents are represented by discrete values.
Using the CA framework, we developed a model for stem-cell-driven tissue growth of the spinal cord. Using this model, computer simulations can be run for predictions about development of biological tissues under different conditions. Below, we describe the steps of model construction in detail.
Our CA model is formulated on a two-dimensional (2D) lattice comprised of square lattice sites. Each lattice site is either empty or associated with an agent, which will be referred to as “cell.” This 2D model (
These two assumptions allow the transformation of the 3D problem into a 2D one. Note that the lattice structure assumes a specific regular, axisymmetric topology to be present inside the 3D tissue. The radial size of 3D lattice sites is fixed, and equal to the size of lattice sites in axial direction x. Hence, we have square lattice sites in the 2D grid (
Rules describing the cellular mechanisms of cell activation, division, differentiation, apoptosis, and phagocytosis are adopted from the simple neurosphere model (Sipahi and Zupanc 2018), with the following significant improvements:
For consistency of notation, throughout this description we denote all probabilities with uppercase letter P and all population pressure values with lowercase letter p.
Each lattice site is either empty or contains a cell whose state can take one of the following values:
The rules describing how these state values are updated in one iteration step, i.e., between iteration numbers t and t+1, are summarized in flowchart (
Table 1 shows the biological properties and rules implemented in the cellular automata model for cell types that are used to study spinal cord growth. For each cell type only those properties are listed that are relevant for model implementation.
As mentioned above, the probabilities mediating cell activation, division and apoptosis are dependent on either population pressure p or radial distance y from the central canal. In our model the former dependence is characterized by two-parameter while the latter by single-parameter functions. In particular, probabilities related to population pressure are characterized by either linear, or exponential functions (for specific formulae see Appendix 1). These Pχ(p), χ=M (Mitosis), and D (Death) functions are parameterized using probabilities Pχ,min and Pχ,max at minimum (p=0) and maximum (p=1−) population pressure values, corresponding to cases when the cell's extended neighborhood is completely empty or completely filled, respectively. Furthermore, the probability of mitosis function PM(p) is characterized such that mitosis is not possible when extended neighborhood of the cell is completely filled. Consequently, PM(1+)=0, always holds which can lead to discontinuity in function PM(p). Contrariwise, function PD(p) is always continuous, thus PD(1+)=PD(1−) always holds. For a particular parameter setting, Pχ(p), χ=M, D functions are plotted in
Functions describing the probability of activation PA(y) and symmetric division PS(y) of stem cells are dependent on radial distance y measured from the central canal. In both cases hyperbolic dependence on y is assumed, which results from the following assumptions.
Due to the axisymmetric nature of the model, the volume of lattice sites increases proportionally with respect to radial distance R=Rc+y+½ measured from the center line of central canal. As a result, the final formula for the probability of activation and symmetric division of stem cells is of the form
where χ=A (Activation), S (Symmetric Division) (for more details on derivation see Appendix 2). The radius of central canal Rc can be provided by experimental data, therefore these functions are characterized by a single parameter Pχ,max which is the maximum probability corresponding to y=0, i.e., to lattice sites located along the surface of the central canal. In
As described above, the probability of mitosis and death of daughter cell are mediated by population pressure p. In our model, population pressure is a quantity that is assigned to the center point of each lattice site. It characterizes the occupancy of the extended neighborhood of an associated lattice site taking into account its distance from cells located inside the extended neighborhood. In particular, neighborhood layers (see
where the occupancy rate is computed as ok=No,k/Nk, with No,k being the number of occupied lattice sites, while Nk is the overall number of lattice sites inside the k-th layer. In addition to the population pressure p, quadrant pressures pi, i=1, 2, 3, 4; are also computed for each lattice site. Quadrant pressures characterize the occupancy rates inside the four quadrants (see
where oi,k=No,i,k/Ni,k is the occupancy rate inside the k-th layer of the i-th quadrant, with No,i,k being the number of occupied lattice sites inside the k-th layer of the i-th quadrant and Ni,k being the total number of lattice sites inside the k-th layer of the i-th quadrant. Note that Ni,k is the same for all k, due to symmetry. Lattice sites along the boundary of two quadrants are assigned to both quadrants, since each lattice site is either occupied or empty.
During division, one of the daughter cells is placed into the lattice site occupied by the mother cell, while the other daughter is placed into one of the four immediate neighbor lattice sites. In the following, the immediate neighbor lattice site chosen for the latter daughter is referred to as target site. In contrast to the method for daughter placement in (Sipahi and Zupanc 2018), where the target site was selected randomly from the empty immediate neighbors, here the immediate neighbor with smallest corresponding quadrant pressure is chosen as target site. The process of daughter placement is illustrated in
It is important to note that, in contrast with (Sipahi and Zupanc 2018), the herein detailed rules allow mitosis even when all four immediate neighbor lattice sites are occupied, as long as there is at least one empty lattice site in the extended neighborhood. When all immediate neighbor sites are occupied then, in order to realize division, the cell at target site must be shifted away. This shifting is carried out according to the following rules:
The process of shifting is illustrated in
Because it is possible that in a time iteration step two or more different cells that undergo mitosis claim the same target site for their daughters, an order must be established for the division of cells that undergo mitosis in the same time iteration step. We arbitrarily defined an order in which daughter placement is followed by shifting. As part of our simulations, we sort, during each time iteration step, all cells that undergo mitosis according to the quadrant pressure corresponding to their target site. We assume that cells with smaller quadrant pressure will undergo mitosis faster, thus we carry out daughter placement and shifting in an ascending order of quadrant pressure. Within a single iteration step, all population pressure and quadrant pressure values are updated after each division (daughter placement and shifting).
In order to carry out simulations using the above detailed lattice and rules, the model is equipped with boundary conditions. Given that the 2D lattice corresponds to an axisymmetric 3D grid, at the lower boundary (i.e., at y=0) the cellular composition of the lattice is mirrored around the axis x. With this boundary condition (BC), it is taken into account that no cells can be found inside the central canal. Note that due to mirroring around x, the rules detailed under Section 4.1.2 cannot result in daughter placement or shifting across x. Consequently, these rules cannot result in the violation of the mirroring BC. In addition to the mirroring BC around x, the occupancy of left, top and right boundaries must be provided, for each time iteration, outside the lattice, within a margin of N1 lattice sites. This is necessary because otherwise population pressure values cannot be computed according to Eqs. (2) and (3).
To make the execution of simulations feasible, an initial cellular composition for the entire lattice is provided. This means the precise specification of location and state values of all agents inside the lattice at the initial time iteration t=0. Note that BCs are also provided at t=0.
In the following, we provide an outline of the algorithm that was developed according to the above detailed model.
for t=0: tend−1
end
end
To demonstrate what results our invention can produce, one exemplary simulation outcome at a particular, final time iteration is shown in
The simulation results in
Tissue growth is based on the behavior of a finite number of discrete cells and the interactions among them. These dynamics are determined by molecular and cellular mechanisms, as well as by stochastic processes, which can be translated into a set of rules and associated probabilities. Thus, for modeling of these growth dynamics, CA provides a framework superior to the traditionally used differential-equation approach because its agents, like biological cells, are finite and distinct, and their local interactions depend on sets of rules. By contrast, in approaches based on differential equations space, time, and state are all continuous. The differential equation framework can be particularly prohibitive to the incorporation of rule-based interactions, as they often lead to jump conditions or non-smoothness in differential equations, resulting in systems that are difficult to simulate. Another major advantage of CA for modeling of tissue growth is that the specific rules and the associated probabilities can be readily modified to adjust for differences in the type or developmental stage of tissue. This flexibility makes the designed model a versatile tool for applications in drug discovery and tissue engineering.
The methods, operations, modules, and systems described herein may be implemented in one or more computer programs executing on a programmable computer system.
Each computer program can be a set of instructions or program code in a code module resident in the random access memory 104 of the computer system 100. Until required by the computer system 100, the set of instructions may be stored in the mass storage device 106 or on another computer system and downloaded via the Internet or other network.
Having thus described several illustrative embodiments, it is to be appreciated that various alterations, modifications, and improvements will readily occur to those skilled in the art. Such alterations, modifications, and improvements are intended to form a part of this disclosure, and are intended to be within the spirit and scope of this disclosure. While some examples presented herein involve specific combinations of functions or structural elements, it should be understood that those functions and elements may be combined in other ways according to the present disclosure to accomplish the same or different objectives. In particular, acts, elements, and features discussed in connection with one embodiment are not intended to be excluded from similar or other roles in other embodiments.
Additionally, elements and components described herein may be further divided into additional components or joined together to form fewer components for performing the same functions. For example, the computer system 100 may comprise one or more physical machines, or virtual machines running on one or more physical machines. In addition, the computer system 100 may comprise a cluster of computers or numerous distributed computers that are connected by the Internet or another network.
Accordingly, the foregoing description and attached drawings are by way of example only, and are not intended to be limiting.
Since the diffusive factor is evenly distributed inside the spinal cord, its overall value in a particular lattice site in the three-dimensional grid is proportional to the volume of the lattice site which can be computed as
V(R)=A(R)Δx, (A 2.1)
where Δx is the size of square lattice sites in the associated two-dimensional grid and
is the area of the transversal section of lattice sites, whose center points are located with R distance from the center line of central canal. The number of angular segments in the three-dimensional model is denoted by Nθ. Due to that in each lattice site the probability of activation and symmetric division of stem cell are inversely proportional to the overall value of diffusive factor inside the lattice site's volume
with χ=A, S and c being a scaling factor. In order to parameterize these functions by their maximum values along the central canal, we require
P
χ(Rc+Δy/2)=Pχ,max, (A 2.4)
where Δy=Δx, due to the square grid. From Eq. (A 2.4) c can be expressed and plugged into Eq. (A 2.3) which becomes
Due to that R=Rc+y+Δy/2, furthermore that Δx=1 is employed in our cellular automata model, one ends up with the formula in Eq. (1).
A shifting route is one shortest route chosen randomly with uniform probability from all possible shortest routes between the target site and the empty lattice site closest to it within the quadrant of daughter placement. It is assumed that routes between the center points of these two lattice sites consist of only horizontal and vertical steps, with each step size being equal to the size of square lattice sites. Consequently, the horizontal steps within shortest routes must sum to the horizontal distance dx, between the two endpoints of these routes. The same holds true for vertical steps where vertical distance of the two endpoints of all possible routes is denoted by dy. Since in our model square lattice sites have unit size (Δx=Δy=1), distances dx and dy are equal to the overall number of horizontal and vertical steps in shortest routes, respectively. In order to find all possible shortest routes, one has to solve a combinatorial partitioning problem, where a total of dx+dy number of steps have to be made consisting of exactly dx number of horizontal and dy number of vertical steps. As a result, the overall number of different shortest routes, i.e., different sequences of horizontal and vertical steps, can be computed as
This application claims priority from U.S. Provisional Patent Application No. 62/928,751 filed on Oct. 31, 2019 entitled Cellular Automata Model of Stem-Cell-Driven Growth of Spinal Cord Tissue, which is hereby incorporated by reference.
This invention was made with government support under Grant No. 1538505 awarded by the National Science Foundation. The government has certain rights in the invention.
Number | Date | Country | |
---|---|---|---|
62928751 | Oct 2019 | US |