Single-cell gene expression time-series simulation and ode-based immune target prediction method

Abstract

Disclosed is a technology for: generating virtual time-series pathway data from RNA sequence data of cells; estimating differentiation pathways of a plurality of sets of cells by using the generated virtual time-series pathway data; completing an ODE model that describes an evolution rule of the cells, by using the estimated differentiation pathways of the plurality of sets of cells; and determining a target drug that induces a state change of the cells, by using the completed ODE model.

Description

Single-cell gene expression time series simulation and ODE-based immune target prediction method

[0001] The present invention relates to a technology for generating virtual time series path data from RNA sequence data of a cell, estimating differentiation paths of a plurality of sets of cells using the generated virtual time series path data, completing an ODE model explaining the evolutionary rules of the cells using the estimated differentiation paths of the plurality of sets of cells, and determining a target drug that induces a state change of the cells using the completed ODE model, and relates to a technology for applying a technology for observing the gene expression level of a cell and a technology for executing a computer simulation using the same to the field of target drugs for cells.

[0002] Existing computational methods that use single-cell RNA-sequencing (scRNA-seq) for cell fate prediction do not model how cells evolve stochastically and in physical time, nor can they predict how differentiation trajectories are altered by proposed interventions.

[0003] Immune evasion has been identified as a new hallmark of cancer, leading to the development of cancer immunotherapies to treat a variety of solid and hematological malignancies. Currently, several major types of immunotherapies are being utilized for cancer treatment, including immune checkpoint inhibitors, T-cell transfer therapy, monoclonal antibodies, therapeutic vaccines, and immune system modulators.

[0004] Immune checkpoint inhibitors, such as CTLA4, PD-1, and PD-L1, have been approved by the FDA to treat several types of cancer, including melanoma and non-small cell lung carcinoma. Checkpoint proteins regulate the balance between pro-inflammatory and anti-inflammatory signals under homeostatic conditions, maintaining immune homeostasis and preventing autoimmunity. Checkpoint inhibitors are a form of immunotherapy designed to prevent T cell dysfunction or exhaustion in the tumor microenvironment. While immunotherapy targeting the tumor microenvironment represents a promising approach to cancer treatment, resistance and treatment refractoriness remain significant challenges.

[0005] Identifying key components that facilitate or hinder the efficacy of checkpoint immunotherapy is essential to comprehensively understand tumor-associated T cells and improve response rates in patients receiving this treatment modality.

[0006] The present invention provides a technology for generating a differentiation trajectory model of a cell in the context of tumor immunity, estimating multiple differentiation trajectories of the cell using the generated model, completing an ODE model using the estimated multiple differentiation trajectories, and determining a target drug that induces a change in the state of the cell using the completed ODE model.

[0007] The concept of the present invention was verified through a demonstration experiment. For the demonstration experiment, the inventors of the present invention retrieved single-cell RNA-sequencing data of melanoma (Accession No. GEO: GSE120575) from the GEO dataset (https: / / www.ncbi.nlm.nih.gov / geo / ) and reanalyzed it. This dataset is a preprocessed dataset. All single-cells classified as CD8+ were isolated to identify distinct CD8+ T-cell clusters. The single-cell RNA-sequencing expression data can be used to analyze the differentiation trajectory of CD8+ T-cells. For downstream analysis, for example, 2,000 genes can be selected from the dataset that showed the greatest changes in expression levels.

[0008] Data augmentation

[0009] Conventional computational methods that utilize single-cell RNA sequencing (scRNA-seq) for cell fate prediction have limitations in that they do not account for the stochastic and physical temporal evolution of cells and cannot predict how proposed interventions will affect differentiation trajectories.

[0010] To overcome these limitations, the present invention integrates time-series scRNA-seq data and modifies the previously proposed PRESCIENT MODEL (Potential expression under Erlying Single Cell gradients MODEL) proposed by Yeo et al. in 2021. The generative modeling framework of the modified model can learn the underlying differentiation landscape, making it particularly useful for scRNA data, which is often scarce and difficult to obtain. The modified PRESCIENT model proposed by the present invention can be used to predict the entire trajectory of a cell based on a few selected time points. The PRESCIENT algorithm requires gene expression data and cell-related metadata.

[0011] A method for modeling cell differentiation according to one aspect of the present invention treats it as a diffusion process, characterized by Brownian motion-based migration of cells from higher to lower potential regions. To understand the underlying dynamic process, the present invention utilizes a generative framework that integrates observations of cell positions and other information obtained from population-level time series data. To simulate the movement of cells from one state to another over time, a stochastic differential equation considering the current state of the cell is used, a drift term describing the direction of movement, and a noise term describing random fluctuations within the process.

[0012] Figure 1 is a visual representation of a generative model for cell differentiation provided according to one aspect of the present invention.

[0013] The future state X(t+1) of a cell is its current state X(t), the drift term μ(X(t))Δt, and the noise term root(2σ 2 Δt)Z(t) is determined by where Δt converges to zero (0). The state X(t+1) of the cell can be expressed as in Equation 1.

[0014] [Formula 1]

[0015] X(t+1) =X(t) + μ(X(t))Δt+ root(2σ 2 Δt)Z(t)

[0016] In a modeling approach provided according to one aspect of the present invention, the drift term, which determines the direction of cell migration, is represented by the negative gradient of the potential function (μ(x) = -∇Φ(x)). The potential function is visualized as a color gradient in the background and is represented numerically in the equation. The gradient acts as a force that drives the cell from a region of high potential to a region of low potential, as indicated by the solid arrows in the model shown in Fig. 1. To account for the stochasticity inherent in the process, the noise term is modeled using a Gaussian distribution, represented by a random variable Z(t) between 0 and 1. The standard deviation of the data set is denoted by σ. The dotted line presented in Fig. 1 corresponds to noise, which contributes to the random fluctuations observed in the process. To learn the trajectories of cells from one state to another, we generate simulations of intermediate states (represented by dotted circles in Fig. 1) corresponding to states between X(t) and X(t+1). By simulating these intermediate states, we can learn the fate of many unobserved states, thereby better understanding the underlying dynamics of the process. The solid circles in Figure 1 represent observed states, and the model learns to predict the trajectories of cells from one state to another by generating simulated values ​​of these intermediate states.

[0017] In Fig. 1, cells having a state presented near the point indicated by reference number 901 may be Naive CD8 T cells, cells having a state presented near the point indicated by reference number 902 may be Cytotoxic CD8 T cells, cells having a state presented near the point indicated by reference number 903 may be Dysfunctional CD8 T cells (Subtype 2), cells having a state presented near the point indicated by reference number 904 may be Pre-dysfunctional CD8 T cells, and cells having a state presented near the point indicated by reference number 905 may be Dysfunctional CD8 T cells (Subtype 1). The above description is intended to help understanding of Fig. 1, and it is to be understood that in some cases it may be different from what is described above.

[0018] Figure 2 illustrates a predictive model fitting process provided according to one aspect of the present invention.

[0019] In one aspect of the present invention, a modeling method is provided that uses a neural network to parameterize a drift function (μ) that determines the direction of cell migration. The neural network uses as input a PCA projection of gene expression data at an observed time point, which is indicated by a solid circle in Figure 2. To simulate a stochastic process, a first-order time discretization is used as shown in Equation 1, and cell states at the next time step (t+Δt) are simulated repeatedly until the next observation time point (t+1). At each iteration, a loss is calculated as the difference between the simulated and observed cell states. The loss is used as feedback to optimize the parameters of the neural network so that the loss is minimized. This iterative optimization process helps the model according to one aspect of the present invention learn the fundamental dynamics of the process and improve its predictive power over time. To achieve this, one aspect of the present invention follows a series of processing steps as presented below.

[0020] Preprocessing

[0021] In the present invention, a gene expression dataset is preprocessed to prepare for training. This preprocessing may include steps such as normalization, feature scaling, or feature engineering, depending on the specific requirements of the dataset and model.

[0022] The predictive model uses multiple inputs, including normalized gene expression data, metadata related to the gene expression data, and information on cell cycle genes and cell apoptosis genes. Growth rates are calculated based on these inputs, which can provide valuable insights into the dynamics of cell differentiation.

[0023] To prepare a data set for training, the present invention performs a series of preprocessing steps. These may include using time points, cell types (e.g., pre-dysfunction, cytotoxicity, dysfunction, dysfunction (heat shock), and cell cycle-related dysfunction cells), growth rates, related gene names, and PCA projections of all gene expression levels. By performing these preprocessing steps, a preprocessed data set suitable for model training can be generated.

[0024] Overall, the inputs and preprocessing steps used in the above prediction model are designed to provide a comprehensive view of the fundamental dynamics of cell differentiation processes, which may help to better understand the factors contributing to cell fate decisions.

[0025] Learning

[0026] The model can be trained using a preprocessed data set for a set number of epochs, which includes running the optimization process for a specified number of iterations. During training, the model learns to predict the trajectories of cells from one state to another based on observed data, and the parameters of the neural network are optimized to minimize a loss function.

[0027] To train a predictive model, a preprocessed dataset is required, as well as a set of hyperparameters that define the learned parameters and architecture of the neural network. Some of the hyperparameters required for training include the number of fully connected layers, the activation function, the number of layers for the latent function, the learning rate, the time step size, the gradient clipping threshold, the batch size, and the number of epochs.

[0028] The above parameter, "Number of fully connected layers," defines the number of fully connected or densely connected layers in a neural network. This can affect the network's ability to model complex relationships in data.

[0029] The above parameter, called 'activation function', specifies the activation function used in each layer of the neural network. Typically, ReLU, sigmoid, and softplus functions can be used.

[0030] The parameter "Number of layers for latent function" above specifies the number of layers used to parameterize the latent function that drives the stochastic process. Deeper potential functions can model more complex dynamics, but may also be more susceptible to overfitting.

[0031] The above parameter called 'learning rate' determines the step size used during optimization and can have a significant impact on the convergence of the learning process.

[0032] The parameter "standard deviation of Gaussian noise" above specifies the amount of noise added during a stochastic process simulation. A larger standard deviation can generate more random trajectories, while a smaller standard deviation can produce more predictable behavior.

[0033] The above parameter, called "time step size," specifies the time interval between simulated states in a stochastic process. While smaller time steps can provide more accurate simulations, they can also increase the computational cost of the learning process.

[0034] The above parameter, called the "gradient clipping threshold," specifies the maximum gradient allowed during training. Clipping the gradient prevents large updates to network weights, which can destabilize the optimization process.

[0035] The above "batch size" parameter specifies the number of training examples used in each batch during the learning process. A larger batch size can lead to faster convergence, but may also increase the memory requirements of the learning process.

[0036] The above parameter, "Number of epochs," determines the number of iterations of the learning process. A higher number of epochs can lead to better convergence and improved performance, but it can also increase the computational cost of learning.

[0037] Simulation

[0038] After the model is trained, it can be used to run various simulations by providing initial conditions to the model and simulating the stochastic process using the parameters learned during training. These simulations can be used to generate predictions about cell fate under various conditions or to explore the dynamics of the process in more detail. In summary, implementing a simulation involves preprocessing the dataset, training the model, and then using the trained model to run the simulation. These steps are essential for accurately modeling and understanding the fundamental dynamics of the cell differentiation process.

[0039] The output of the above training step is a PyTorch model containing weights for a model capable of simulating accurate cell fate trajectories at different time points. This step allows for simulations of specific cell types (e.g., naive cells) or all cell types. The simulation also requires a hyperparameter regarding the number of steps. This step provides sufficient data to train an additional functional unit (based on CellBox) for the perturbation analysis task, allowing for simulating the cell fates of naive cells for multiple time steps. Because naive cells from pre-treatment responders provide better perturbation results than other cell types or patient types (e.g., responders, non-responders, etc.) after treatment, the simulated output of naive cells can be utilized as ground truth for the dynamic system module according to the present invention.

[0040] <Dynamic System Modeling Part>

[0041] The goal of the dynamical system modeling part is to generate a weight matrix of size NxN for N genes. These weights represent interaction parameters between each pair of genes and help predict cell behavior under unseen perturbations.

[0042] In the present invention, the model can be trained at several different points in time and the interaction parameters can be calculated in an unperturbed state.

[0043] Figure 3 illustrates a workflow of a dynamic system modeling part provided according to one aspect of the present invention.

[0044] According to the present invention, an interpretable machine learning model can be constructed using the systematic responses of cells over a certain period of time. The present invention models system behavior in terms of interaction parameters connecting genes (nodes) using the ordinary differential equation (ODE) of Equation 2.

[0045] [Formula 2]

[0046] ∂x i (t) / ∂t=ε i φ(Σ j≠1 W ij x j (t))-α i x i (t)

[0047] x in formula 2 i (t) represents the log2 normalized gene expression of gene i at time t.

[0048] In Equation 2, ∂x i (t) / ∂t represents the log2-normalized time relative change in gene expression of gene i.

[0049] Ordinary differential equations (ODEs) are fundamental to mathematical modeling, clarifying the relationships between biological variables and their evolution over time. This method is particularly important for understanding the dynamic interactions within biological systems, particularly in understanding drug interactions, absorption, and the overall impact on disease progression. Here, "evolution" refers to time evolution. In the field of drug target discovery, ODE models facilitate complex simulations of biological pathways and disease mechanisms. They help identify potential targets by predicting how changes in these pathways may affect disease outcomes. This technology is particularly valuable in chronic diseases such as cancer, where understanding the long-term impact of therapeutic interventions is crucial. The use of ODE models to identify novel immunotherapy targets represents a significant advancement in the field. This approach provides a systematic framework for analyzing the complex web of interactions within biological systems, enabling the identification of potential targets that may otherwise be hidden by traditional analytical methods. This pioneering strategy can significantly enhance the impact of drug discovery and pave the way for the development of targeted immunotherapy strategies.

[0050] In the present invention, the sigmoid function φ(x)=tanh(x) is used to introduce the saturation effect of the interaction term, thereby converting the sigmoid function to ε i The limit on the absolute value of the contribution of the interaction term to the derivative is partly motivated by the fact that the modeling according to the present invention models only a portion of the cellular components. In the course of the present invention, RELU was also tested with the sigmoid function, but the difference between RELU and tanh is not significant. The interaction parameter w ijrepresents the interactions between gene j and gene i, and is assumed to be a constant property of a pair of molecules in a given cellular environment. Since nodes / genes cannot act on their own, self-interactions are not allowed, and therefore w ii =0 is not allowed. Here x j (t) represents the log2 normalized gene expression of gene j at time t. α i It characterizes the decay effect, which refers to the tendency of gene i to return to its original gene expression level before the state change. The model can be iteratively trained by optimizing the interaction parameters so that the numerically simulated system response matches the ground truth.

[0051] Implementation

[0052] As input data, gene names, log-normalized gene expression values ​​at the initial time point, and log-normalized gene expression values ​​at the Nth time point (e.g., N=500) (generated using the data augmentation module according to the present invention) can be used as ground truth. The Wij weight matrix is ​​initialized with all self-interactions set to 0. Learning without prior knowledge leads to unbiased perturbation results, which can yield better results. Conventional techniques can be used as ODE calculators, and first-order ODEs can be used to reduce the complexity of the method. A 70:30 split between training and test data can be used. The training and test data can be randomly split. 20% of the training data can be used as a validation set. The model according to the present invention can be repeatedly trained by optimizing interaction parameters so that the predicted gene expression results match the observed data. As a loss function, an L1 regularization penalty for network density and the Euclidean distance between the numerical simulation results of the ODE model and the experimental data can be selected to avoid overfitting. In the present invention, a model can be trained for, for example, 6 epochs for learning. Here, each epoch has a fixed number of iterations (1000), and the learning rate can decrease for each epoch.

[0053] Verification

[0054] For validation purposes, one can verify how closely the predicted gene expression of the prediction model according to the present invention matches the actual data. Various criteria, such as monitoring the mean squared error during training, can be used for this purpose. As training progresses, the overall error decreases.

[0055] Perturbation

[0056] Figure 4 is a conceptual diagram illustrating a method for executing perturbation according to one embodiment of the present invention.

[0057] In the present invention, all genes are perturbed for 400 time steps with a fixed perturbation intensity (left side of Figure 4), and then gene expression is output after each gene perturbation passes through an SVM model (middle side of Figure 4). From the SVM results, a state dictionary can be created that indicates the number of each cell type (Naive, Pre-Dysfunction, Dysfunction, Cytotoxic, etc.) present after each gene perturbation (right side of Figure 4).

[0058] Dynamic Simulation

[0059] Figure 5 illustrates a dynamic prediction modeling method provided according to one embodiment of the present invention.

[0060] After learning the dynamic system modeling part provided according to one embodiment of the present invention, the parameter (ε) for each gene, the optimized decay value (α) for each gene, and the optimized interaction parameter (w ij ) can be obtained. Using these parameters, the perturbation response of a specific gene can be simulated.

[0061] Support Vector Machines for Cell Classification

[0062] FIG. 6 is a diagram illustrating a support vector machine provided according to one embodiment of the present invention.

[0063] To analyze the perturbation results, it is necessary to determine whether the cell state has changed. Therefore, a Support Vector Machine (SVM) model can be designed to classify the gene expression of cells. The initial gene expression and related metadata of each cell can be used as input for learning. Cells can be classified into the following six cell types based on gene expression.

[0064] Type 1: Naive Cell

[0065] Type 2: Pre-Dysfunctional Cell

[0066] Type 3: Cytotoxic Cell

[0067] Type 4: Cell Cycle Dysfunction Cell

[0068] Type 5: (Heat Shock) Dysfunction Cell

[0069] Type 6: Dysfunction Cell (not related to cell cycle and heat shock)

[0070] Target Score

[0071] Figure 7 visually depicts the goals of the perturbation model of the present invention.

[0072] It is desirable to secure a greater number of cytotoxic T cells and a smaller number of dysfunctional or metastatic (pre-dysfunction) T cells. The balance between the two cell types can be significantly altered by intervention.

[0073] In the present invention, the SVM model can be used to classify and count the number of each type of cell for each genetic perturbation. For the control value, the perturbation condition is not executed to first count the number of cells of each type. The control value serves as a reference for the rate performance of each genetic perturbation. In the melanoma disease model, the ratio of the number of cytotoxic lymphocytes to the total number of dysfunctional cells (C / D) can be calculated.

[0074] Figure 8 is a table showing the results of genetic perturbation according to one embodiment of the present invention.

[0075] In Figure 8, inhibition of the 16 genes from above shifts the trajectory toward activated T cells, and inhibition of the 9 genes from below shifts the trajectory toward dysfunctional T cells.

[0076] If genetic perturbation can convert some dysfunctional cells into cytotoxic cells, the gene is considered a positive marker and potential target. However, if some cytotoxic cells are converted to dysfunctional or pre-dysfunctional cells, the C / D ratio decreases, and the gene is considered a negative marker. Negative markers can also utilize genes that downregulate their activity, making them useful for drug discovery. Inhibiting these genes upregulates the negative markers, increasing cytotoxic activity and the C / D ratio. As shown in the figure below, genes are ranked by their C / D ratio, with higher C / D ratios indicating higher rankings. No perturbation is used as a control. PDCD1 is a known marker. When inhibited / disrupted, the C / D ratio increases, the number of cytotoxic cells increases from 2224 to 2365, and the total number of dysfunctional cells decreases from 529 to 119. When TCF7 is inhibited / disrupted, the C / D ratio decreases, and the number of cytotoxic cells decreases from 2224 to 2187. Since only the terminal states (cytotoxicity and dysfunction) are of interest, and since pre-dysfunction is a transition state and Naive cells are a starting state, we are only interested in the pre-dysfunction state and time, so pre-function toxicity and Naive cells are not considered for ranking.

[0077] In this specification, the models and / or target drug determination methods provided by the present invention are generalized mathematical frameworks designed to study cell fate and differentiation across a wide range of cell types, including immune cells, neural cells, epithelial cells, stem cells, cancer cells, and other cell types. This flexibility allows the present invention to be applied to a variety of biological contexts, particularly including the following cell types:

[0078] First, the present invention can be applied to immune cells.

[0079] Second, the present invention can be applied to neural cells. The model can be applied to neurogenesis research, where understanding how neural stem cells differentiate into various types of neurons or glial cells is crucial. This could aid in developing treatments for neurological diseases.

[0080] Third, the present invention can be applied to epithelial cells. Models according to the present invention can explore how epithelial cells in various organs respond to injury and participate in tissue repair, which is crucial for treating conditions such as fibrosis and chronic wounds.

[0081] Fourth, the present invention can be applied to stem cells. The model according to the present invention can be used to predict the fate of stem cells in various niches critical to the advancement of regenerative medicine and stem cell therapy, allowing researchers to better utilize these cells for therapeutic purposes.

[0082] Fifth, the present invention can be applied to cancer cells. In addition to immune cells, the model according to the present invention can help understand tumor evolution by exploring the progression and heterogeneity of cancer cells within tumors and develop strategies to overcome treatment resistance.

[0083] Each of the examples described above demonstrates the capabilities of the model of the present invention in integrating and analyzing complex biological data. The present invention can provide insights into cellular dynamics that are important for both basic research and clinical applications. The predictive capabilities of the model provided by the present invention, combined with its ability to process diverse data types, make it a valuable tool for the scientific community, potentially transforming approaches to understanding and treating various diseases.

[0084] According to one aspect of the present invention, there is provided a method comprising: obtaining time series data representing states of a plurality of cells by observing gene expression in a plurality of cells composed of single cells; determining a drift term representing a force acting on the evolution of a cell's state according to the state of the single cell using the obtained time series data; independently executing a differentiation trajectory determination process for multiple times, which determines a differentiation trajectory representing the state evolution of the single cell by applying the drift term and a predetermined random noise term together to a state that the single cell can have, wherein the random noise term represents a probabilistic force acting on the evolution of the state of the single cell regardless of the state of the single cell; determining an ODE model by calculating values ​​of parameters constituting an ODE model that describes an evolution rule of the single cell based on the determined plurality of differentiation trajectories; A method for determining a target drug may be provided, including a step of determining a target drug that induces a state change of the single cell using the determined ODE model.

[0085] At this time, the plurality of cells may be T cells.

[0086] At this time, the plurality of cells may be any one type of cell among neural cells, epithelial cells, stem cells, cancer cells, and T cells.

[0087] At this time, the step of acquiring the time series data is to obtain N virtual times (pt) for each of the plurality of cells. n ) by tagging one virtual time among the virtual time, and defining N cell groups based on the virtual time (n=1, ..., N); and determining a representative value of the gene expression level of cells included in each of the cell groups; wherein the time series data may be composed of N representative values ​​of the gene expression levels determined for the N cell groups.

[0088] At this time, the representative value of the gene expression level determined for the nth cell group among the N cell groups is K elements (e k ) is an array including elements (k=1, ..., K), and the value of the kth element among the K elements may be a representative value of the expression levels of the kth gene of cells belonging to the nth cell group.

[0089] At this time, the representative value of the expression levels of the k gene of cells belonging to the n-th cell group may be the average value of the expression levels of the k gene of cells belonging to the n-th cell group.

[0090] At this time, the defining step may include: a step of observing the state of the plurality of cells at one specific time; a step of observing any one cell among the plurality of cells to determine the evolutionary state of the any one cell; a step of selecting one virtual time from the N virtual times based on the determined evolutionary state; and a step of tagging the selected virtual time to the any one cell.

[0091] At this time, the function of the drift term is implemented as a neural network, and the neural network may be configured to output the next state (X(t+Δt)) of the single-cell at a point in time after a predetermined unit time Δt when the current state (X(t)) of the single-cell at the current point in time is input. And the method for learning the neural network includes the steps of: obtaining a first observation state (X(t1)), which is a state of the single-cell observed at a first point in time (t1), and a second observation state (X(t2)), which is a state of the cell observed at a second point in time (t2) after the first point in time (t1); inputting the first observation state (X(t1)) to the neural network, and inputting the output of the neural network back to the neural network, repeating the process at least once, to produce a second predicted state (X`(t2)), which is a state of the cell predicted at the second point in time (t2) (provided that t2-t1>Δt); and a step of updating the parameters of the neural network to minimize the loss between the second observation state (X(t2)) and the second prediction state (X`(t2)).

[0092] At this time, the above parameters may be parameters that explain the interaction of genes of the single cells among the variables of the ODE (ordinary differential equation) model.

[0093] At this time, the step of independently executing the differentiation path determination process multiple times may include a step of applying the drift term to the first state of the single cell to determine the first change state of the single cell; and a step of applying the random noise term to the first change state of the single cell to determine the second change state of the single cell.

[0094] At this time, among the plurality of differentiation paths, the first differentiation path may include the first state of the single-cell and the first probabilistic state obtained by applying the random noise term to the first state, and among the plurality of differentiation paths, the second differentiation path may include the first state of the single-cell and the second probabilistic state obtained by applying the random noise term to the first state.

[0095] At this time, the step of determining the ODE model may include a step of optimizing the values ​​of parameters constituting the ODE model by supervised learning the ODE model to optimize a loss function between a plurality of gene phenotypes simulated by the plurality of differentiation paths and a gene phenotype predicted by the ODE model.

[0096] According to another aspect of the present invention, there is provided a method comprising: obtaining time series data representing states of a plurality of cells, each of which is composed of a single cell; determining a drift network, which determines an amount of change in a state of a single cell according to the state of the single cell, using the obtained time series data; independently executing a differentiation trajectory determination process, which determines a differentiation trajectory representing state evolution of the single cell by applying a state change by the drift network and a predetermined random state change together to a state that the single cell can have, wherein the random state change represents a probabilistic state change that acts on the evolution of the state of the single cell regardless of a specific state of the single cell; determining an ODE model by calculating values ​​of parameters constituting an ODE model that describes an evolution rule of the single cell based on the determined plurality of differentiation trajectories; A method for determining a target drug may be provided, including a step of determining a target drug that induces a state change of the single cell using the determined ODE model.

[0097] According to another aspect of the present invention, a drug administration method may be provided, comprising a step of administering a target drug determined by the aforementioned target drug determination method to a subject. The subject may be a human or an animal.

[0098] The present invention provides a technique for modeling cell differentiation trajectories in the context of tumor immunity. The present invention also provides a technique for modeling CD8+ T-cell differentiation trajectories, enabling analysis of CD8+ T-cells to predict the success or failure of checkpoint immunotherapy.

[0099] According to the present invention, a technology can be provided for generating a differentiation trajectory model of a cell in the context of tumor immunity, estimating multiple differentiation trajectories of the cell using the generated model, completing an ODE model using the estimated multiple differentiation trajectories, and determining a target drug that induces a change in the state of the cell using the completed ODE model.

[0100] Figure 1 is a visual representation of a generative model for cell differentiation provided according to one aspect of the present invention.

[0101] Figure 2 illustrates a predictive model fitting process provided according to one aspect of the present invention.

[0102] Figure 3 illustrates a workflow of a dynamic system modeling part provided according to one aspect of the present invention.

[0103] Figure 4 is a conceptual diagram illustrating a method for executing perturbation according to one embodiment of the present invention.

[0104] Figure 5 illustrates a dynamic prediction modeling method provided according to one embodiment of the present invention.

[0105] FIG. 6 is a diagram illustrating a support vector machine provided according to one embodiment of the present invention.

[0106] Figure 7 visually depicts the goals of the perturbation model of the present invention.

[0107] Figure 8 is a table showing the results of genetic perturbation according to one embodiment of the present invention.

[0108] Figure 9 is a flowchart showing a target drug determination method provided according to one embodiment of the present invention.

[0109] FIG. 10 is a flowchart illustrating a method for obtaining time series data presented in FIG. 9 according to one embodiment of the present invention.

[0110] FIG. 11 illustrates a concept of defining multiple cell groups by tagging a virtual view for each of multiple cells according to one embodiment of the present invention.

[0111] Figure 12 illustrates an implementation form of a representative value of gene expression levels of cells included in each cell group defined according to one embodiment of the present invention.

[0112] FIG. 13 illustrates an example of an implementation method for determining the value of the kth element among an array representing representative values ​​of gene expression levels of cells included in a specific cell group defined according to one embodiment of the present invention.

[0113] FIG. 14 is a flowchart illustrating a method for tagging a plurality of cells with virtual vision and dividing the plurality of cells into N cell groups according to one embodiment of the present invention.

[0114] FIG. 15 is a flowchart illustrating a learning method of a neural network that provides a drift term function according to one embodiment of the present invention.

[0115] FIG. 16 is a block diagram illustrating the main components of a single cell gene expression analysis device provided according to one embodiment of the present invention.

[0116] FIG. 17 is a system schematic diagram showing the data flow between an analysis device and a target drug determination system provided according to one embodiment of the present invention.

[0117] Figure 18 is a diagram illustrating a structure in which single cells are captured in a microfluidic-based cell separation module.

[0118] Figure 19 is a conceptual diagram showing the process of barcoding single-cell derived cDNA and measuring gene expression signals.

[0119] Figure 20 is a schematic diagram illustrating the overall configuration and operation flow of a decision support system provided according to one embodiment of the present invention.

[0120] Hereinafter, embodiments of the present invention will be described with reference to the attached drawings. However, the present invention is not limited to the embodiments described herein and may be implemented in various other forms. The terminology used herein is intended to aid understanding of the embodiments and is not intended to limit the scope of the present invention. Furthermore, the singular forms used below also include the plural forms, unless the context clearly indicates otherwise.

[0121] Figure 9 is a flowchart showing a target drug determination method provided according to one embodiment of the present invention.

[0122] FIG. 10 is a flowchart illustrating a method for obtaining time series data presented in FIG. 9 according to one embodiment of the present invention.

[0123] The following description is provided with reference to FIGS. 9 and 10.

[0124] In step (S110), gene expression within a plurality of cells consisting of single cells can be observed to obtain time series data representing the states of the plurality of cells.

[0125] In one preferred embodiment, the single cell and the plurality of cells may be T cells.

[0126] Step (S110) may include steps (S111) and (S112) of FIG. 10.

[0127] In step (S111), N virtual views (pt) are provided for each of the plurality of cells. n ) can be tagged with one virtual time, and N cell groups can be defined based on the virtual time (n=1, ..., N). That is, cells tagged with the same virtual time belong to the same group.

[0128] In step (S112), a representative value of the gene expression level of cells included in each of the above cell groups can be determined.

[0129] At this time, the time series data may be composed of representative values ​​of the N gene expression levels determined for the N cell groups.

[0130] In step (S120), using the acquired time series data, a drift term representing a force acting on the evolution of a cell's state according to the state of the single cell can be determined.

[0131] In step (S130), a differentiation trajectory determination process that determines a differentiation trajectory representing the state evolution of the single cell by applying the drift term and a predetermined random noise term together to a state that the single cell can have can be independently executed multiple times.

[0132] At this time, the random noise term represents a probabilistic force that acts on the evolution of the state of the single cell regardless of the state of the single cell.

[0133] In one embodiment, step (S130) may include a step (S131) ​​of applying the drift term to the first state of the single-cell to determine a first change state in which the single-cell has changed, and a step (S132) of applying the random noise term to the first change state of the single-cell to determine a second change state in which the single-cell has changed.

[0134] At this time, among the plurality of differentiation paths, the first differentiation path may include the first state of the single-cell and the first probabilistic state obtained by applying the random noise term to the first state, and the second differentiation path may include the first state of the single-cell and the second probabilistic state obtained by applying the random noise term to the first state. That is, even if the single-cell is dominated by one drift term, since it is also influenced by the random noise term, a specific state of the single-cell may evolve into a different state by the random noise term.

[0135] In this specification, the term "drift term" may be replaced with the term "drift network." Accordingly, the drift network may be a network that determines the amount of change in the state of a single cell based on the state of the single cell, using the acquired time-series data. The drift network or drift term is represented by μ in Figure 2.

[0136] In step (S140), based on the determined plurality of differentiation paths, the ODE model can be determined by calculating the values ​​of parameters that constitute the ODE model explaining the evolution rule of the single cell. In other words, the ODE model can be completed.

[0137] At this time, before step (S140) is executed, the number of nodes constituting the ODE model and the structural information of the ODE model, which is information about the protein that each node represents, are already provided. Step (S140) is a process for determining parameters that describe the interactions between the nodes of the ODE model. At this time, the parameters may be parameters that describe the interactions between genes of the single cells among the variables of the ODE (ordinary differential equation) model.

[0138] The step of determining the ODE model may include a step of supervised learning the ODE model to optimize a loss function between the plurality of gene phenotypes simulated by the plurality of differentiation pathways and the gene phenotypes predicted by the ODE model, thereby optimizing the values ​​of parameters constituting the ODE model. A related concept is presented in Fig. 3.

[0139] In step (S150), a target drug that induces a state change of the single cell can be determined using the determined ODE model.

[0140] With the above ODE model determined, a technique for determining a target drug that induces a state change in the single cell can be implemented by adopting a related technique disclosed in the past.

[0141] FIG. 11 illustrates a concept of defining multiple cell groups by tagging a virtual view for each of multiple cells according to one embodiment of the present invention.

[0142] The status of multiple cells, each consisting of a single cell, cultured in a culture dish (100) of Fig. 11 can be individually observed. If a total of N virtual views are defined, the multiple cells can be divided into at most N groups.

[0143] Figure 12 illustrates an implementation form of a representative value of gene expression levels of cells included in each cell group defined according to one embodiment of the present invention.

[0144] Among the existing genes, a total of K genes of interest, i.e., gene1, gene2, gene3, ..., and geneK, can be selected.

[0145] As shown in Fig. 12, the representative value of the gene expression level determined for the nth cell group among the N cell groups is divided into K elements (e k) may be an array including (k=1, ..., K). For example, in Fig. 12, the representative values ​​of gene expression levels of cells included in the first cell group (G_1) may be an array of {0.0, 0.121, 0.0, ..., 0.0}. At this time, the value of the k-th element among the K elements is a representative value of the expression levels of the k-th gene of cells belonging to the n-th cell group.

[0146] In the example presented in Fig. 12, the numbers of cells included in the first cell group, the second cell group, the third cell group, ..., the n-th cell group, ..., and the N-th cell group are exemplified as 400, 1000, 700, ..., 460, ..., and 320, respectively, and these specific numbers are arbitrarily selected to help understanding.

[0147] FIG. 13 illustrates an example of an implementation method for determining the value of the kth element among an array representing representative values ​​of gene expression levels of cells included in a specific cell group defined according to one embodiment of the present invention.

[0148] In Fig. 13, the first cell group (G_1) presented in Fig. 12 is illustrated. As presented in Fig. 12, the first cell group (G_1) includes a total of 400 cells, all of which are assigned a virtual time tp1. The number of genes activated in each of the 400 cells is different. For example, in the cell with cell number (1), gene 2 is activated by 0.121, but in the cell with cell number (2), gene 2 is activated by 0.097. The K numbers presented in the bottommost row of Fig. 13 are representative values ​​of the gene expression levels of cells belonging to the first cell group (G_1) presented in Fig. 12. The representative values ​​of the 400 numbers presented in each column of Fig. 13 are the numbers written in the bottommost row of the corresponding column.

[0149] In one embodiment, the representative value of the expression levels of the kth gene (gene k) of cells belonging to the nth cell group may be the average value of the expression levels of the kth gene of cells belonging to the nth cell group.

[0150] FIG. 14 is a flowchart illustrating a method for tagging a plurality of cells with virtual vision and dividing the plurality of cells into N cell groups according to one embodiment of the present invention.

[0151] The steps included in the method of FIG. 14 may be the steps included in step (S111) of FIG. 10.

[0152] In step (S1111), the status of the plurality of cells can be observed at one specific time.

[0153] In step (S1112), any one of the plurality of cells can be observed to determine the evolutionary state of any one cell.

[0154] In step (S1113), one virtual time point among the N virtual times points can be selected based on the determined evolution state.

[0155] In step (S1114), the selected virtual time can be tagged to any one of the cells.

[0156] The above steps (S1112) to (S1114) can be repeated for all of the plurality of cells, thereby tagging an appropriate virtual time for all of the plurality of cells.

[0157] The above drift term function can be implemented as a neural network.

[0158] The above neural network may be configured to output the next state (X(t+Δt)) of the single cell at a point in time after a predetermined unit time Δt, when the current state (X(t)) of the single cell at the current point in time is input.

[0159] Here, the current state (X(t)) may be a representative value of the gene expression level of cells included in the first cell group corresponding to a virtual time (ex: t is pt1 presented in FIG. 11) representing the current state among the plurality of cell groups, and may be an actually observed value. For example, in the current state (X(t)), t may be pt1 presented in FIG. 11. And the next state (X(t+Δt)) may be a simulated value based on the current state (X(t)), as a representative value of the gene expression level of cells included in the first cell groups, at a time slightly later than the virtual time (ex: t+Δt = pt1+Δt) representing the current state.

[0160] FIG. 15 is a flowchart illustrating a learning method of a neural network that provides a drift term function according to one embodiment of the present invention.

[0161] In step (S210), a first observation state (X(t1)), which is a state of the single cell observed at a first time point (t1), and a second observation state (X(t2)), which is a state of the cell observed at a second time point (t2) after the first time point (t1), can be obtained.

[0162] In step (S220), the process of inputting the first observation state (X(t1)) into the neural network and inputting the output of the neural network back into the neural network is repeated at least once to produce the second predicted state (X`(t2)), which is the state of the cell predicted at the second time point (t2) (provided that t2-t1>Δt).

[0163] In step (S230), the parameters of the neural network can be updated to minimize the loss between the second observation state (X(t2)) and the second prediction state (X`(t2)).

[0164]

[0165] According to one embodiment of the present invention, a single-cell gene expression analysis device specifically designed to precisely measure gene expression at the single-cell level is provided. This device is designed to accurately measure temporal and spatial changes in gene expression levels, thereby effectively reconstructing the dynamic evolutionary trajectory of cellular states.

[0166] The device comprises, first, a microfluidic-based cell separation module configured to individually capture multiple single cells. This module independently moves each cell within a microfluidic channel and positions them in individual reaction zones, enabling analysis without mutual interference. Second, it comprises a biochemical reaction module for extracting RNA from each cell and reverse-transcription (cDNA synthesis) of the RNA in real time. This module includes cell disruption and RNA purification steps, and the reaction is performed at low temperatures using non-heated enzymes to prevent RNA degradation.

[0167] Third, a barcoding module is provided to assign unique identification information for each cell to the reverse-transcribed cDNA. This barcoding is performed in conjunction with the gene amplification (PCR) step, enabling the simultaneous analysis of multiple cell-derived data within a single sequencing flow. Fourth, a measurement module is included to convert gene expression levels into electrical or optical signals for measurement. This measurement module may, for example, include optical intensity measurements of fluorescently dyed nucleic acid probes or electrical resistance-based RNA quantity detection technology.

[0168] To control each of the above modules and provide the measured gene expression data to a subsequent analysis system in real time, the device includes a dedicated control unit, which is implemented as a special-purpose computing unit based on an ASIC (Application Specific Integrated Circuit) or an FPGA (Field Programmable Gate Array) rather than a general-purpose computer. This dedicated control unit is optimized for high-speed parallel processing operations, enabling high-speed time-series data processing, and can be used to improve the accuracy of ODE modeling and drug target prediction in subsequent stages.

[0169] According to the present invention, using the single cell gene expression analysis device as described above, high-resolution gene expression levels at the single cell level can be accurately measured, and time series analysis and differentiation path prediction can be performed based on the measurement data.

[0170] FIG. 16 is a block diagram illustrating the main components of a single cell gene expression analysis device provided according to one embodiment of the present invention.

[0171] FIG. 17 is a system schematic diagram showing the data flow between an analysis device and a target drug determination system provided according to one embodiment of the present invention.

[0172] Figure 18 is a diagram illustrating a structure in which single cells are captured in a microfluidic-based cell separation module.

[0173] Figure 19 is a conceptual diagram showing the process of barcoding single-cell derived cDNA and measuring gene expression signals.

[0174] The following description is provided with reference to FIGS. 16 to 19.

[0175] The single-cell gene expression analysis device (10) according to the present embodiment is configured to include a microfluidic-based cell separation module (110), a biochemical reaction module (120), a barcoding module (130), a measurement module (140), and a dedicated control unit (150), as illustrated in FIG. 16.

[0176] The cell separation module (110) can independently capture cells within individual fluid droplets or chambers by introducing multiple single cells into the microfluidic channel. The fluidic channel is designed to have a structure that enables stable separation without interference between cells depending on the flow rate and viscosity of the fluid (see FIG. 18).

[0177] The biochemical reaction module (120) is configured to automatically extract RNA from each captured cell and reverse transcribe the RNA into cDNA. This module includes a low-temperature maintenance system to minimize RNA degradation and a microreactor structure containing a non-thermal reverse transcriptase enzyme.

[0178] The barcoding module (130) functions to assign a cell-specific molecular barcode to the reverse-transcribed cDNA. As illustrated in FIG. 19, this barcode is used to distinguish the origin data of each cell in the subsequent parallel sequencing step.

[0179] The measurement module (140) includes a fluorescence-based detector or an electrical resistance-based detector, and is equipped with an electrical or optical sensor for measuring the amount of cDNA. For example, the intensity of light generated according to the amount of cDNA reacting with the fluorescent probe is measured using a photodetector, or the change in charge is detected through a microelectrode to quantify the expression amount.

[0180] The dedicated control unit (150) is comprised of an FPGA-based high-speed parallel operation circuit, preprocesses the measurement results in real time, and transmits the data to a target drug determination server (20) capable of ODE model learning and prediction calculation. Unlike general-purpose computers, the dedicated control unit is specially designed to rapidly process large-scale genetic vectors without signal delay.

[0181] The device is configured to complete the entire measurement and data transfer process within 10 minutes, offering greater time efficiency compared to conventional bulk sequencing-based methods. Furthermore, since cell separation and data transfer are automatically performed within a single device in a consistent pipeline, biological contamination or cross-interference can be minimized.

[0182] FIG. 20 is a schematic diagram showing the overall configuration and operation flow of a decision support system (30) provided according to one embodiment of the present invention.

[0183] The above system (30) includes various components to support medical staff's treatment decisions based on gene expression data and target drug prediction results provided from a single cell gene expression analysis device.

[0184] The decision support system (30) first receives analysis data from a gene expression analysis device or a target drug determination system via an analysis result receiving module (310). The analysis data may include, for example, a gene expression vector measured at the single cell level (=gene expression measurement results), ODE-based simulation results, a response index for a specific gene perturbation, or a list of target drug candidates. The received data is format-converted or preprocessed so that it can be interpreted and integrated and analyzed in subsequent modules.

[0185] Next, the patient data linkage module (320) links with the hospital information system (HIS) or electronic medical record (EMR) system to retrieve the patient's basic information (e.g., gender, age, diagnosis), medical history, current medications, allergy information, and other clinical data in real time. This patient information is combined with analysis data and used as a key reference for drug eligibility assessment and risk assessment.

[0186] Thereafter, the clinical efficacy assessment module (330) quantitatively evaluates the clinical validity of each drug candidate based on analysis results and patient information. This module (330) references external drug databases (e.g., FDA approval history, ClinicalTrials.gov, WHO ATC codes, etc.) and clinical literature-based risk indicators to calculate an efficacy score or priority for each drug candidate, comprehensively considering safety, response rate, and potential side effects. These evaluation results are used as the basis for personalized drug recommendations.

[0187] The visualization interface (340) includes a user interface for intuitively presenting the overall results to medical professionals. The visualization screen displays, for example, cell type distribution, a list of expressed genes, the ratio of dysfunctional T cells and cytotoxic T cells, simulation results of changes in the C / D ratio of drug candidates, and clinical efficacy scores in diagrammatic or graphical form. Furthermore, a selection button and a "reason for non-selection" field are included for each drug, allowing medical professionals to select or exclude drugs based on their own judgment and record the rationale for their selection.

[0188] Each module of the above system (30) operates organically and interoperably, supporting medication decision-making by providing real-time analysis-based information based on the experience-based judgment of medical staff. Specifically, by integrating and processing analysis results based on single-cell expression data with patient-specific information, a far more personalized treatment strategy can be established than conventional, standardized treatment protocols.

[0189] Additionally, the decision support system (30) may also include functions to ensure clinical safety, such as automatically linking and storing selected medications and their reasons for selection in the EMR system, or generating warning alerts for high-risk medications. This configuration goes beyond simply conveying analysis results and serves a practical function as a medical technology system that plays a technical role in the actual clinical decision-making process.

[0190] By utilizing the embodiments of the present invention described above, those skilled in the art will be able to easily implement various changes and modifications without departing from the essential characteristics of the present invention. The content of each claim may be combined with other claims that are not in a citation relationship within the scope of this specification, as long as it is understood.

Claims

1. A step of a single-cell gene expression analysis device observing gene expression within a plurality of cells composed of single cells to generate time series data representing the states of the plurality of cells; A step in which the target drug decision server acquires the generated time series data; A step in which the target drug decision server determines a drift term representing a force acting on the evolution of a cell's state according to the state of the single cell using the acquired time series data; A step of independently executing a differentiation trajectory determination process for determining a differentiation trajectory representing the state evolution of the single-cell by applying the drift term and a predetermined random noise term together to a state that the single-cell can have, wherein the random noise term represents a probabilistic force acting on the state evolution of the single-cell regardless of the state of the single-cell; A step of determining the ODE model by the target drug decision server calculating the values ​​of parameters constituting the ODE model that describes the evolution rule of the single cell based on the determined plurality of differentiation paths; A step in which the target drug decision server determines a target drug that induces a state change of the single cell using the determined ODE model; The step of the target drug decision server transmitting information about the determined target drug to a decision support system for treatment decision; and A step in which the decision support system integrates and analyzes information about the target drug with the patient's medical information, evaluates the clinical effectiveness of the target drug, and visually provides the information to medical staff; including, Method for determining target drugs.

2. A method for determining a target drug in the first paragraph, wherein the plurality of cells are T cells.

3. A method for determining a target drug in the first paragraph, wherein the plurality of cells are any one of neural cells, epithelial cells, stem cells, and cancer cells.

4. In paragraph 1, The step of obtaining the above time series data is: For each of the above multiple cells, N virtual visual angles (pt) n ) by tagging one of the virtual times, and defining N cell groups based on the virtual times (n=1, ..., N); and A step of determining a representative value of the gene expression level of cells included in each of the above cell groups; Includes, The above time series data is characterized in that it is composed of representative values ​​of the N gene expression levels determined for the N cell groups. Method for determining target drugs.

5. In paragraph 4, The representative value of the gene expression level determined for the nth cell group among the N cell groups is K elements (e k ) is an array containing (k=1, ..., K), Among the above K elements, the value of the kth element is a representative value of the expression levels of the kth gene of cells belonging to the nth cell group. Method for determining target drugs.

6. A method for determining a target drug, wherein the representative value of the expression levels of the k gene of cells belonging to the n-th cell group is the average value of the expression levels of the k gene of cells belonging to the n-th cell group.

7. In paragraph 4, The steps defined above are: A step of observing the status of the above plurality of cells at one specific time; A step of observing any one cell among the plurality of cells and determining the evolutionary state of the any one cell; A step of selecting one virtual time from among the N virtual times based on the above-determined evolutionary state; and A step of tagging the selected virtual time to any one of the cells; including, Method for determining target drugs.

8. In paragraph 2, The above drift term function is implemented as a neural network, The above neural network is configured to output the next state (X(t+Δt)) of the single cell at a point in time after a predetermined unit time Δt when the current state (X(t)) of the single cell at the current point in time is input. The method of training the above neural network is as follows: A step of acquiring a first observation state (X(t1)), which is a state of the single-cell observed at a first time point (t1), and a second observation state (X(t2)), which is a state of the single-cell observed at a second time point (t2) after the first time point (t1); A step of inputting the first observation state (X(t1)) into the neural network and inputting the output of the neural network back into the neural network at least once to produce the second predicted state (X`(t2)), which is the state of the single cell predicted at the second time point (t2) (provided that t2-t1>Δt); and A step of updating the parameters of the neural network to minimize the loss between the second observation state (X(t2)) and the second prediction state (X`(t2)); including, Method for determining target drugs.

9. A method for determining a target drug, wherein the parameters in the first paragraph are parameters that describe the interaction of genes of the single cells among the variables of the ODE (ordinary differential equation) model.

10. In paragraph 1, The step of independently executing the above differentiation path determination process multiple times is: A step of determining a first change state of the single cell by applying the drift term to the first state of the single cell; A step of determining a second change state of the single cell by applying the random noise term to the first change state of the single cell; including, Method for determining target drugs.

11. In paragraph 1, Among the plurality of differentiation paths, the first differentiation path includes the first state of the single cell and the first probabilistic state obtained by applying the random noise term to the first state, Among the plurality of differentiation paths, the second differentiation path includes the first state of the single cell and the second probabilistic state obtained by applying the random noise term to the first state. Method for determining target drugs.

12. In paragraph 1, The step of determining the ODE model includes a step of supervising learning the ODE model to optimize a loss function between a plurality of gene phenotypes simulated by the plurality of differentiation paths and a gene phenotype predicted by the ODE model, thereby optimizing the values ​​of parameters constituting the ODE model. Method for determining target drugs.

13. A step of generating time series data representing the states of a plurality of cells by observing gene expression in a plurality of cells composed of single cells by a single-cell gene expression analysis device; A step in which the target drug decision server acquires the generated time series data; A step in which the target drug decision server determines a drift network that determines the amount of change in the state of the single cell according to the state of the single cell using the acquired time series data; A step of independently executing a differentiation trajectory determination process for determining a differentiation trajectory representing the state evolution of the single cell by applying a state change by the drift network and a predetermined random state change together to a state that the single cell can have, wherein the random state change represents a probabilistic state change that acts on the state evolution of the single cell regardless of the specific state of the single cell; A step of determining the ODE model by the target drug decision server calculating the values ​​of parameters constituting the ODE model that describes the evolution rule of the single cell based on the determined plurality of differentiation paths; A step in which the target drug decision server determines a target drug that induces a state change of the single cell using the determined ODE model; The step of the target drug decision server transmitting information about the determined target drug to a decision support system for treatment decision; and A step in which the decision support system integrates and analyzes information about the target drug with the patient's medical information, evaluates the clinical effectiveness of the target drug, and visually provides the information to medical staff; including, Method for determining target drugs.

14. A method for determining a target drug in claim 13, wherein the plurality of cells are any one of neural cells, epithelial cells, stem cells, cancer cells, and T cells.

15. A target drug decision system comprising a single cell gene expression analysis device (10), a target drug decision server (20), and a decision support system (30), The above single cell gene expression analysis device (10) is It is configured to observe gene expression within a plurality of cells composed of single cells and generate time series data representing the state of the plurality of cells, The above target drug decision server (20) is Obtain the above generated time series data, Using the above time series data, a drift term that acts on the state evolution according to the state of the single cell is determined, Multiple differentiation paths are generated by applying the above drift term and a predetermined random noise term together, Based on the above differentiation pathways, the parameters of the ordinary differential equation (ODE) model that explains the evolutionary rules of the single cell are derived, Using the above ODE model, a target drug that induces a state change of the single cell is determined, It is configured to transmit information about the above-determined target drug to the decision support system (30). The above decision support system (30) is The information on the above target drug is integrated and analyzed with the patient's medical information, By evaluating the clinical efficacy of the above target drug, Characterized in that it is configured to visually provide the above evaluation results to medical staff, Target drug decision system.

Images