A virtual organoid-oriented drug response digital twin prediction system

Abstract

The application discloses a kind of virtual organ-oriented drug response digital twin prediction system, it is related to drug response prediction technical field.The system includes the following modules: data acquisition module obtains the multi-modal observation sequence and causal pathway graph of organoid, generates trajectory tensor;Initial graph construction module encodes causal pathway graph into causal skeleton and assigns node state value, generates initial state graph;Structure refinement module filters causal effect edge through backdoor adjustment and transfer entropy estimation, generates refined state graph;Counterfactual intervention module executes do operator intervention to generate initial counterfactual trajectory;Target point preliminary screening module calculates marginal contribution in reverse to screen candidate target points;Target point verification module verifies target point effectiveness through virtual knockout;Report generation module generates drug response mechanism interpretation report.The application explicitly encodes biological causal pathway into prediction model, realizes the leap from statistical association prediction to causal effect reasoning, and improves the credibility and interpretability of drug target discovery.

Description

Technical Field

[0001] This invention relates to the field of drug response prediction technology, and in particular to a digital twin prediction system for drug response for virtual organoids. Background Technology

[0002] Organoids, as three-dimensional cell structures cultured in vitro, can effectively simulate the structure and function of organs in vivo, providing more accurate drug response predictions than traditional two-dimensional cell cultures in drug screening. With the development of digital twin and artificial intelligence technologies, researchers have begun to construct virtual organoid models, integrating multimodal longitudinal observation data (such as mass spectrometry signals, image sequences, and electrical impedance time series) to simulate the drug response process of organoids in computational space. Existing technical solutions mainly employ deep learning models to extract features from big data and establish statistical mapping relationships between inputs and outputs to predict drug responses. Chinese invention patent CN121031367A proposes a digital twin drug data analysis system that achieves drug response prediction through a data fusion module and a multi-layer digital twin modeling module.

[0003] However, existing virtual organoid drug response prediction systems suffer from the following main shortcomings. System modeling heavily relies on deep learning to extract statistical features from data. While these models excel at capturing statistical correlations between variables, they cannot distinguish between correlation and causation—when a system predicts the effectiveness of a drug, the identified "key target" may only be a statistically correlated variable rather than a true causal driver. The observational nature of organoid experimental data dictates that we can only observe "the outcome after intervention," not "what would have happened if there had been no intervention" (i.e., counterfactual), making it fundamentally difficult to learn causal relationships from purely observational data. The "interpretability" of existing systems largely remains at the level of feature importance analysis (such as SHAP values), failing to reach the level of causal inference. The "mechanism explanations" provided by the system may mislead subsequent experimental design, resulting in a waste of R&D resources. Due to the lack of causal interpretability, the credibility of such prediction models in drug target discovery and regulatory approval is limited, downgrading them from "scientific discovery tools" to "black box predictors." Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a digital twin prediction system for drug response in virtual organoids. It aims to solve the technical problems of existing virtual organoid drug response prediction systems lacking causal reasoning capabilities and being unable to explain the causes of drug response at the mechanistic level. By constructing a digital twin prediction system that explicitly encodes biological causal pathways, it achieves a leap from statistical association prediction to causal effect reasoning.

[0005] To achieve the above objectives, the present invention provides the following technical solution: a digital twin prediction system for drug response in virtual organoids, comprising the following steps: The data acquisition module is used to acquire multimodal observation sequences and causal pathway maps of organoids, and to align the vectors in the multimodal longitudinal observation sequences to a unified space to obtain the trajectory tensor. The initial graph construction module is used to encode the causal pathway graph into a causal skeleton. It assigns the state values ​​of each node in the trajectory tensor and the drug intervention node labels and endpoint response node labels in the causal pathway graph to the corresponding nodes of the causal skeleton to generate the initial state graph. The structure refinement module is used to perform backdoor adjustments on the initial state graph and estimate the transfer entropy between each node. It retains the nodes with connections greater than a preset threshold as causal effect edges, removes the nodes with connections less than a preset threshold, and generates a refined state graph containing drug intervention nodes and endpoint response nodes. The counterfactual intervention module is used to execute do operator intervention with the drug intervention node of the refined state diagram as the intervention point, and propagate the intervention effect along the causal effect edge node by node to obtain the initial counterfactual trajectory; The target screening module is used to calculate the marginal contribution of each causal effect edge from the end response node of the refined state diagram in reverse order, and select the nodes connected by the causal effect edges that are sorted in descending order to form a set of candidate targets. The target verification module is used to virtually knock out each node in the candidate target set in sequence and lock it to zero, re-execute the do operator intervention to obtain the trajectory after knockout, compare the difference in the state of the endpoint response node before and after knockout, remove the nodes whose difference is less than the preset difference threshold, and output the refined candidate target set. The report generation module is used to generate a drug response mechanism explanation report based on the refined candidate target set and the initial counterfactual trajectory.

[0006] Furthermore, the acquisition of the multimodal observation sequence and causal pathway map of organoids, aligning the vectors in the multimodal longitudinal observation sequence to a unified space to obtain the trajectory tensor, includes: Organoids in culture were collected using liquid chromatography-mass spectrometry, a high-content imaging system, and a microelectrode array sensor within a continuous observation window to obtain raw mass spectrometry signal streams, raw image sequences, and raw electrical impedance time series. Peak detection and retention time alignment were performed on the raw mass spectrometry signal stream to generate a molecular expression abundance matrix indexed by the observation time point; The original image sequence is subjected to three-dimensional structural reconstruction and cell nucleus segmentation to generate an organoid spatial morphology matrix indexed by the observation time point; Baseline drift correction and low-frequency filtering were performed on the original electrical impedance time series to generate an organoid barrier function matrix indexed by the observation time points; The molecular expression abundance matrix, organoid spatial morphology matrix and organoid barrier function matrix are longitudinally spliced ​​according to the observation time points to form a multimodal longitudinal observation sequence. The causal pathway graph is obtained from a public pathway database. The protein, gene and metabolite entities in the causal pathway graph are mapped to gene symbol identifiers to serve as identifiers for subsequent nodes, thereby generating a standardized causal pathway graph. Using the node order of the standardized causal pathway graph as the arrangement basis, the row vectors of the molecular expression abundance matrix, the row vectors of the organoid spatial morphology matrix, and the row vectors of the organoid barrier function matrix corresponding to each observation time point in the multimodal longitudinal observation sequence are concatenated into a unified state vector in the feature dimension. The node sequence of the unified state vector is rearranged according to the node arrangement order of the standardized causal pathway graph to generate the rearranged unified state vector for each observation time point. The unified state vectors after rearranging the observation time points are stacked in chronological order to form a three-dimensional tensor, which is the trajectory tensor.

[0007] Furthermore, the step of encoding the causal pathway graph into a causal skeleton, assigning the state values ​​of each node in the trajectory tensor and the drug intervention node markers and endpoint response node markers in the causal pathway graph to the corresponding nodes of the causal skeleton to generate an initial state graph, includes: Extract all node names from the causal path graph and generate a node name sequence; extract the directed connections between nodes in the causal path graph and generate a set of directed edges. Construct an empty adjacency matrix using the node name sequence as row and column indices, set the main diagonal elements to zero, assign 1 to non-main diagonal elements with connections and 0 to the rest based on the set of directed edges, and use the adjacency matrix as a causal skeleton. Using the row index order of the causal skeleton as the traversal order, the node name corresponding to the row index is matched with the gene symbol identifier in the second dimension of the trajectory tensor; The third-dimensional feature channel values ​​of the trajectory tensor corresponding to the successfully matched gene symbol are expanded into a one-dimensional state sequence along the first-dimensional observation time point. The one-dimensional state sequence is written into the node storage area of ​​the causal skeleton row by row with the node name as the sequence index, generating the node state matrix of the causal skeleton. The row index of the node state matrix is ​​the observation time point, and the column index is the node name. Write the drug intervention node marker and the endpoint response node marker into the attribute field of the corresponding node in the causal skeleton, indexed by the node name; Assign each row vector of the node state matrix to the node at the corresponding row index in the causal skeleton to generate the initial state graph.

[0008] Furthermore, the backdoor adjustment of the initial state graph and estimation of the propagation entropy between nodes, retaining nodes with connections greater than a preset threshold as causal effect edges, removing nodes with connections not greater than a preset threshold, and generating a refined state graph containing drug intervention nodes and endpoint response nodes, includes: Extract the drug intervention node marker and the endpoint response node marker from the initial state graph, and use them as the drug intervention node and the endpoint response node respectively. Starting from the drug intervention node and ending from the endpoint response node, traverse all directed paths from the starting point to the endpoint in the initial state graph, and extract all intermediate nodes passed through all directed paths to generate an intermediate node set. Each intermediate node in the intermediate node set is taken out as the node to be tested. The conditional mutual information value of the drug intervention node state time series and the endpoint response node state time series under the given state time series of the node to be tested is calculated. The node to be tested with a conditional mutual information value less than the preset mutual information threshold is marked as a redundant node. All redundant nodes are removed from the intermediate node set to generate a candidate adjustment variable set. For each pair of directed edges in the initial state graph, obtain the state time series of the starting node and the state time series of the ending node of the edge. Using the state time series of each node in the candidate adjustment variable set as condition variables, the joint probability density of the state time series of the endpoint node, the starting node, and the condition variables is estimated, and the ratio of the joint probability density to the marginal probability density is integrated to obtain the conditional mutual information value, which is then used as the transit entropy value. Directed edges between nodes with a propagation entropy value greater than a preset entropy threshold are retained as causal effect edges. Directed edges between nodes with a propagation entropy value less than or equal to the preset entropy threshold are removed from the initial state graph. The initial state graph after removing redundant edges is used as the refined state graph.

[0009] Furthermore, the step of executing the do operator intervention using the drug intervention node in the refined state diagram as the intervention point, propagating the intervention effect node by node along the causal effect edge, and obtaining the initial counterfactual trajectory includes: Using the drug intervention nodes in the refined state graph as the intervention objects, the state values ​​of the drug intervention nodes at each observation time point are replaced with preset intervention dose values, all directed edges pointing to the drug intervention nodes in the refined state graph are deleted, and a post-intervention causal graph is generated. Starting with the drug intervention node in the post-intervention causal graph, a topological sort is performed on the post-intervention causal graph to generate a node processing sequence. Each node is extracted in sequence according to the node processing sequence. When a node is extracted, the observed state value of the node at each observation time point is read from the refined state graph. The counterfactual state values ​​of all parent nodes of the node are used as independent variables, and the observed state values ​​are used as dependent variables. Linear regression fitting is performed on the independent and dependent variables to obtain the structural equation coefficients and intercepts of the node. The structural equation coefficients and intercepts are combined into a linear equation, which is then used as the structural equation for that node. Substitute the counterfactual state values ​​of all parent nodes at the current observation time into the structural equation of the node, calculate the counterfactual state value of the node at the current observation time, and repeat this step in the order of observation time until the counterfactual state values ​​of the node at all observation time points have been calculated. The counterfactual state values ​​of each node at each observation time point are stored in the counterfactual state matrix according to the observation time point and node identifier, and the counterfactual state matrix is ​​used as the initial counterfactual trajectory.

[0010] Furthermore, the step of calculating the marginal contribution of each causal effect edge from the endpoint response node of the refined state graph backwards, and selecting the nodes connected by the causal effect edges that are sorted in descending order to form a set of candidate target points, includes: Using the endpoint response node in the refined state diagram as the target node, the difference between the state value of the endpoint response node at the last observation time point in the initial counterfactual trajectory and the observed state value of the endpoint response node at the last observation time point in the trajectory tensor is taken as the total effect value. Extract all parent nodes of the endpoint response node from the refined state diagram, remove each parent node one by one, replace the state value of the removed parent node with zero at all observation time points, and then re-execute the do operator intervention to obtain the counterfactual trajectory after the virtual removal of the parent node. Extract the state value of the endpoint response node at the last observation time point in the trajectory, and take the difference between the state value of the endpoint response node at the last observation time point in the initial counterfactual trajectory and the state value after virtual removal as the marginal contribution value of the parent node to the endpoint response node. Sort all parent nodes in descending order of their marginal contribution values, extract the first preset number of parent nodes from the sorted list, and use the extracted parent nodes as the candidate target set.

[0011] Furthermore, the process involves sequentially and virtually knocking out each node in the candidate target set and locking it to zero, then re-executing the do operator to obtain the knockout trajectory. The difference in the state of the endpoint response nodes before and after knockout is compared, and nodes with a difference less than a preset difference threshold are removed. The refined candidate target set is then output, including: Using the candidate target set as the input set, an empty set is established as the refined candidate set, and each node in the candidate target set is taken out as the current node to be processed. The refined state graph is copied, and the state values ​​of the current node to be processed in the copied graph are replaced with zero at all observation time points, while the state values ​​of the other nodes remain unchanged. The copied graph is used as a virtual knockout state graph. Using the drug intervention node in the virtual knockout state graph as the intervention point, delete all directed edges pointing to the drug intervention node in the virtual knockout state graph, execute the do operator intervention, and obtain the knockout trajectory corresponding to the current node to be processed. Extract the initial state value of the endpoint response node at the last observation time from the initial counterfactual trajectory, extract the knockout state value of the endpoint response node at the last observation time from the knockout trajectory, and perform a subtraction operation between the initial state value and the knockout state value to obtain the difference; If the difference is greater than or equal to the preset difference threshold, the node identifier of the current node to be processed is stored in the refined candidate set; if the difference is less than the preset difference threshold, the storage operation is not performed on the refined candidate set. The next node in the candidate target set is taken as the current node to be processed. This process continues until all nodes in the candidate target set have been processed. The refined candidate set is then output as the refined candidate target set.

[0012] Furthermore, the step of extracting the initial state value of the endpoint response node at the last observation time point from the initial counterfactual trajectory, extracting the knockout state value of the endpoint response node at the last observation time point from the knockout trajectory, and performing a subtraction operation between the initial state value and the knockout state value to obtain the difference includes: Using the initial counterfactual trajectory as the first input matrix, the node identifier of the endpoint response node as the first row index, and the sequence number of the last observation time point as the first column index, the matrix element values ​​corresponding to the first row index and the first column index are extracted from the first input matrix as the initial state values. Using the knockout trajectory as the second input matrix, the node identifier of the endpoint response node as the second row index, and the sequence number of the last observation time point as the second column index, the matrix element values ​​corresponding to the second row index and the second column index are extracted from the second input matrix and used as the knockout state values. The difference is obtained by subtracting the initial state value from the state value after knocking it out.

[0013] Furthermore, the generation of a drug response mechanism explanation report based on the refined candidate target set and the initial counterfactual trajectory includes: Using the node identifier of each node in the refined candidate target set as the row index, the counterfactual state values ​​of each node at all observation time points are extracted from the initial counterfactual trajectory to generate the target counterfactual state matrix. The arithmetic mean of each observation time point is calculated column by column for the counterfactual state matrix of the target point to obtain the population average response curve at each observation time point. The standard deviation of each observation time point is calculated column by column to obtain the population response fluctuation band. Using the node identifier of each node in the refined candidate target set as an index, the transfer entropy value of each node is extracted, and the transfer entropy value is used as the causal effect strength of each target. Divide the marginal contribution value of each target by the causal effect strength of each target, and use the quotient obtained by the division operation as the effect confidence score of each target. Arrange the targets in descending order and mark the sorting number. With the observation time point as the horizontal axis and the counterfactual state value as the vertical axis, an individual response curve for each target point is generated, and the population average response curve and the population response fluctuation band are superimposed on the individual response curve. The individual response curves, marginal contribution values, transfer entropy values, effect confidence scores, and ranking numbers of each target are combined one by one to generate drug response mechanism entries for each target. All drug response mechanism entries for all targets are concatenated in sequential order, and the concatenation result is output as a drug response mechanism interpretation report.

[0014] The present invention provides a digital twin prediction system for drug response in virtual organoids, the beneficial effects of which are mainly reflected in the following aspects: 1. This invention achieves a leap from statistical association prediction to causal effect inference by encoding biological causal pathway graphs into a causal skeleton and embedding it into a prediction model. The system's initial graph construction module assigns drug intervention node labels and endpoint response node labels from the causal pathway graph to corresponding nodes in the causal skeleton. The structure refinement module filters causal effect edges through backdoor adjustment and transit entropy estimation, enabling the system to incorporate known biological prior knowledge during the modeling stage. This mechanism allows the system to answer counterfactual questions such as "How will the organoid's response to drug Y change if gene X is knocked out," rather than simply providing a binary judgment of "whether to respond."

[0015] 2. This invention employs a counterfactual intervention module that uses drug intervention nodes as entry points to execute do-operator intervention, propagating the intervention effect node by node along the causal effect edge to generate an initial counterfactual trajectory. A target screening module calculates the marginal contribution of each causal effect edge backward from the endpoint response node, filtering the candidate target set in descending order. A target verification module compares the state difference of the endpoint response node before and after knockout through virtual knockout and locking to zero operations, eliminating nodes with differences less than a preset difference threshold. The collaborative work of these modules ensures that the candidate targets output by the system possess causal credibility, rather than being based solely on statistical correlation. The report generation module, based on the refined candidate target set and the initial counterfactual trajectory, generates a drug response mechanism explanation report containing individual target response curves, marginal contribution values, propagation entropy values, and effect credibility scores. This provides drug developers with verifiable and traceable causal mechanism explanations, helping to reduce the proportion of invalid targets entering the experimental verification stage, reduce research resource waste, and improve the efficiency of drug target discovery. Attached Figure Description

[0016] Figure 1This is a functional block diagram of a digital twin prediction system for drug response oriented to virtual organoids according to the present invention; Figure 2 This is a flowchart of the structural refinement module of the present invention; Figure 3 This is a flowchart of the target screening module of the present invention. Detailed Implementation

[0017] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0018] Please see Figures 1-3 This invention provides a digital twin prediction system for drug response in virtual organoids, comprising the following steps: The data acquisition module is used to acquire multimodal observation sequences and causal pathway maps of organoids, and to align the vectors in the multimodal longitudinal observation sequences to a unified space to obtain the trajectory tensor. In this embodiment, the step of obtaining the multimodal observation sequence and causal pathway map of organoids, and aligning the vectors in the multimodal longitudinal observation sequence to a unified space to obtain the trajectory tensor, includes: Organoids in culture were collected using liquid chromatography-mass spectrometry, a high-content imaging system, and a microelectrode array sensor within a continuous observation window, yielding raw mass spectrometry signal streams, raw image sequences, and raw electrical impedance time series.

[0019] Specifically, the liquid chromatography-mass spectrometry (LC-MS) system introduces the effluent from the liquid chromatography separation into a mass spectrometer detector to acquire raw data on the changes in ion intensity of metabolites and proteins in organoid culture medium and cell lysates over time. Its output is a sequence of total ion current chromatograms. The high-content imaging system is a cell imaging device integrating a wide-field fluorescence microscope and an automated stage. It is used to acquire bright-field and multi-channel fluorescence microscopic images of organoids under multiple fields of view, and its output is a sequence of two-dimensional or three-dimensional image frames arranged according to observation time points.

[0020] A microelectrode array sensor is a detection device that integrates multiple platinum or gold microelectrodes on a culture substrate. It is used to apply an alternating voltage to the organoid barrier structure and measure the response current. Its output is a continuous recording sequence of the impedance amplitude and phase angle of each electrode channel over time. A continuous observation time window refers to the complete time span from the start of the pre-drug baseline observation period to the end of the post-drug effect observation period. Within this window, three types of raw data are collected sequentially at equal time intervals.

[0021] Peak detection and retention time alignment were performed on the raw mass spectrometry signal stream to generate a molecular expression abundance matrix indexed by the observation time points.

[0022] Specifically, peak detection refers to identifying the chromatographic peak position, peak height, and peak area of ​​each molecular ion in the raw mass spectrometry signal stream, converting the continuous ion intensity signal into discrete molecular characteristics and corresponding abundance measurements. Retention time alignment refers to correcting the peak positions acquired multiple times to a common retention time scale, using internal standards or dynamic time warping algorithms, to address the shift in retention time of the same molecule in different injections caused by column aging or fluctuations in mobile phase composition. After peak detection and retention time alignment, each observation time point corresponds to a set of molecular identifiers and their abundance values. The row index represents the observation time point, the column index represents the molecular identifier, and the table content represents the abundance value of the corresponding molecule at that time point. This two-dimensional table is the molecular expression abundance matrix.

[0023] The original image sequence is subjected to three-dimensional structural reconstruction and cell nucleus segmentation to generate an organoid spatial morphology matrix indexed by the observation time point.

[0024] Specifically, 3D structure reconstruction refers to restoring a series of multi-layered 2D image frames acquired at the same observation time point along the optical axis into a 2D image superposition form where each voxel has a grayscale value, i.e., 3D voxel data, through deconvolution or confocal stacking algorithms. Cell nucleus segmentation refers to using DAPI or Hoechst fluorescent markers in the 3D voxel data, and employing watershed or level set algorithms to individually label the voxel region occupied by each cell nucleus as an independent object.

[0025] After three-dimensional structural reconstruction and cell nucleus segmentation, each observation time point corresponds to multiple morphological feature values ​​such as organoid volume, surface area, total number of cell nuclei, average cell nucleus volume, and centroid coordinates of each cell nucleus. The row index is the observation time point, the column index is the identifier of each type of morphological feature, and the table content is the value of the corresponding feature at that time point. This two-dimensional table is the organoid spatial morphology matrix.

[0026] Baseline drift correction and low-frequency filtering were performed on the original electrical impedance time series to generate an organoid barrier function matrix indexed by the observation time points.

[0027] Specifically, baseline drift correction refers to using a polynomial fitting method to estimate the slowly changing trend components in the original impedance time series and subtracting them, restoring the impedance signal to the zero-mean level. Low-frequency filtering refers to using a Butterworth low-pass filter with a set cutoff frequency to remove noise components above the cutoff frequency caused by culture medium flow or mechanical vibration, while retaining low-frequency impedance change signals that reflect cell barrier integrity and the state of tight cell junctions.

[0028] After correction and filtering, the impedance amplitude and phase angle values ​​of each electrode channel at each observation time point are obtained. The row index is the observation time point, the column index is the electrode channel identifier and the corresponding impedance characteristic type, and the table content is the value of the impedance characteristic of the channel at that time point. This two-dimensional table is the organoid barrier function matrix.

[0029] The molecular expression abundance matrix, organoid spatial morphology matrix, and organoid barrier function matrix are longitudinally spliced ​​according to the observation time points to form a multimodal longitudinal observation sequence.

[0030] Specifically, vertical concatenation refers to using the observation time point as a common reference for row indices, and sequentially connecting all feature column vectors belonging to the same time point in the three matrices along the column direction to form a new row vector containing all features of the three modalities. The new row vectors of all observation time points are stacked in chronological order to form a two-dimensional data table with the number of rows equal to the total number of observation time points and the number of columns equal to the sum of the number of features of the three types. This data table is the multimodal vertical observation sequence.

[0031] The causal pathway graph is obtained from a public pathway database. The protein, gene and metabolite entities in the causal pathway graph are mapped to gene symbol identifiers to serve as identifiers for subsequent nodes, thereby generating a standardized causal pathway graph.

[0032] Specifically, public pathway databases refer to open-access databases such as KEGG or Reactome that store information on species-specific signal transduction pathways and their compositional relationships. Causal pathway graphs refer to the network structure composed of nodes and directed edges extracted from these databases; nodes represent molecules involved in the signaling pathway, and directed edges represent activation or inhibition relationships between molecules. Gene symbol identifiers refer to the standard nomenclature symbols designated by the Human Genome Nomenclature Committee to uniquely identify genes and their encoded proteins, such as TP53 representing the gene encoding the p53 protein.

[0033] The mapping process involves precisely matching and synonymizing the molecular names, aliases, and metabolite names stored in various databases within the causal pathway graph with a lexicon of gene symbol identifiers, uniformly replacing them with the corresponding gene symbol identifiers. After mapping, all nodes in the causal pathway graph are named using gene symbol identifiers, and this graph is the standardized causal pathway graph.

[0034] Using the node order of the standardized causal pathway graph as the arrangement basis, the row vectors of the molecular expression abundance matrix, the row vectors of the organoid spatial morphology matrix, and the row vectors of the organoid barrier function matrix corresponding to each observation time point in the multimodal longitudinal observation sequence are concatenated into a unified state vector in the feature dimension. The node sequence of the unified state vector is then rearranged according to the node arrangement order of the standardized causal pathway graph to generate the rearranged unified state vector for each observation time point.

[0035] Specifically, the node order in a standardized causal pathway graph refers to the sequence of node names taken sequentially from the entire node list of the graph according to a fixed sorting rule. This sequence serves as a unified reference for all subsequent vector arrangement operations. For a given observation time point, the row vectors of that time point in the molecular expression abundance matrix, the row vectors in the organoid spatial morphology matrix, and the row vectors in the organoid barrier function matrix are concatenated end-to-end to form an initial row vector containing all modal feature data, i.e., a unified state vector.

[0036] The rearrangement operation uses the node order of the standardized causal pathway graph as an index to extract the feature sub-vectors corresponding to each gene symbol identifier from the unified state vector, so that the node arrangement order of the final output vector is completely consistent with the node order in the standardized causal pathway graph. The above splicing and rearrangement operation is performed independently at each observation time point, and a rearranged unified state vector is generated at each observation time point.

[0037] The unified state vectors after rearranging the observation time points are stacked in chronological order to form a three-dimensional tensor, which is the trajectory tensor.

[0038] Specifically, the stacking operation refers to rearranging all observation time points and then arranging the unified state vectors sequentially along the first time dimension to form a three-dimensional numerical array with the first dimension being the total number of observation time points, the second dimension being the total number of nodes in the standardized causal path graph, and the third dimension being the total number of feature channels of each node in multimodal observation. This three-dimensional array is the trajectory tensor.

[0039] The initial graph construction module is used to encode the causal pathway graph into a causal skeleton. It assigns the state values ​​of each node in the trajectory tensor and the drug intervention node labels and endpoint response node labels in the causal pathway graph to the corresponding nodes of the causal skeleton to generate the initial state graph. In this embodiment, encoding the causal pathway graph into a causal skeleton, and assigning the state values ​​of each node in the trajectory tensor and the drug intervention node markers and endpoint response node markers in the causal pathway graph to the corresponding nodes of the causal skeleton to generate an initial state graph, includes: Extract all node names from the causal path graph and generate a node name sequence. Extract the directed connections between nodes in the causal path graph and generate a set of directed edges.

[0040] Specifically, a node name refers to the gene symbol identifier corresponding to each signal transduction molecule in the standardized causal pathway graph. A node name sequence is an ordered list obtained by lexicographically sorting all node names in ascending order, ensuring the sequence order remains constant during each execution. A directed connection refers to a directional regulatory edge in the standardized causal pathway graph pointing from an upstream molecule node to a downstream molecule node. Each regulatory edge includes a start node identifier and a termination node identifier. The set of directed edges refers to the data set that summarizes all directed connections in the form of ordered binary pairs consisting of start node identifiers and termination node identifiers.

[0041] Construct an empty adjacency matrix using the node name sequence as row and column indices, set the main diagonal elements to zero, assign 1 to non-main diagonal elements with connections based on the set of directed edges, and assign 0 to the rest, and use the adjacency matrix as a causal skeleton.

[0042] Specifically, an empty adjacency matrix is ​​a two-dimensional square matrix with both rows and columns equal to the total number of nodes, and all elements are initially set to zero. The row and column indices are the node name sequences arranged lexicographically as described above. Setting the main diagonal elements to zero means fixing the elements with the same row and column indices to zero, indicating that no node has a control edge pointing to itself. When assigning values ​​based on the directed edge set, for each ordered binary pair in the set, the element at the intersection of the column index of the row corresponding to the starting node and the column corresponding to the ending node is set to 1, indicating that a connection exists in that direction. All other elements in the matrix that do not correspond to any directed edges remain zero. This assigned adjacency matrix serves as the causal framework.

[0043] Using the row index order of the causal skeleton as the traversal order, the node name corresponding to the row index is matched with the gene symbol identifier in the second dimension of the trajectory tensor.

[0044] Specifically, the row index order refers to the order in which the row names in the causal skeleton adjacency matrix are arranged, i.e., the lexicographical order mentioned above. The traversal order refers to the order in which each row name is processed sequentially from the first row name to the last row name. The second dimension of the trajectory tensor stores a list of gene symbol identifiers for each node. This list records the molecular identity corresponding to each feature channel according to the node arrangement order when the trajectory tensor was generated. String matching refers to using the node name corresponding to the current row index as the comparison object, and sequentially performing case-sensitive exact equality string comparison operations in the gene symbol identifier list of the second dimension of the trajectory tensor until an equal entry is found or the comparison of all entries in the list is completed.

[0045] The third-dimensional feature channel values ​​of the trajectory tensor corresponding to the successfully matched gene symbol are expanded into a one-dimensional state sequence along the first-dimensional observation time point. The one-dimensional state sequence is written into the node storage area of ​​the causal skeleton row by row with the node name as the sequence index, generating the node state matrix of the causal skeleton. The row index of the node state matrix is ​​the observation time point, and the column index is the node name.

[0046] Specifically, a successful match means that the node name corresponding to the current row index is exactly equal to a gene symbol identifier in the second dimension of the trajectory tensor, confirming that both point to the same molecular entity. The third dimension feature channel value of the trajectory tensor refers to all modal observation values ​​recorded by this node at each observation time point, including measurement data of three channels: molecular expression abundance, organoid spatial morphology, and barrier function. Expanding along the first dimension observation time points means taking out the third dimension channel values ​​of all observation time points belonging to the same node in chronological order and rearranging them into a one-dimensional numerical sequence.

[0047] The node storage area refers to the contiguous storage location reserved for each node in the causal skeleton data structure to store dynamic state data. Row-by-row writing means filling the storage area of ​​the corresponding node in the causal skeleton with the one-dimensional state sequence of each node in the order of the node name sequence. The state sequences of all nodes form a two-dimensional table with the observation time point as the row index and the node name as the column index. This table is the node state matrix.

[0048] Write the drug intervention node tag and the endpoint response node tag into the attribute field of the corresponding node in the causal skeleton, indexed by the node name.

[0049] Specifically, drug intervention node markers refer to pre-labeled node identifiers in the standardized causal pathway diagram that carry information about the direct target of drug action. Endpoint response node markers refer to pre-labeled identifiers in the standardized causal pathway diagram that represent the terminal biological function nodes used to characterize drug effects. Using node names as an index involves performing string matching between the drug intervention node markers and endpoint response node markers and the list of node names in the causal framework, respectively, to determine the row and column positions of their corresponding nodes within the causal framework.

[0050] After writing the drug intervention node tag into the attribute field of the matching node, the type attribute of the node is set to intervention type. After writing the endpoint response node tag into the attribute field of the matching node, the type attribute of the node is set to response type.

[0051] Assign each row vector of the node state matrix to the node at the corresponding row index in the causal skeleton to generate the initial state graph.

[0052] Specifically, each row vector in the node state matrix refers to all the column values ​​contained in each row of the node state matrix. The row index of each row corresponds to an observation time point, and the values ​​corresponding to each column index in that row are the state values ​​of each node at that observation time point. Assigning values ​​to the corresponding row index in the causal skeleton means writing the values ​​in each row vector into the state storage location of the node with the same name in the causal skeleton according to the order of the observation time points, so that each node obtains its dynamic state sequence at all observation time points. After all row vectors are assigned values, each node in the causal skeleton carries its own state time sequence and type attribute within the complete observation time window. This data structure is the initial state diagram.

[0053] The structure refinement module is used to perform backdoor adjustments on the initial state graph and estimate the transfer entropy between each node. It retains the nodes with connections greater than a preset threshold as causal effect edges, removes the nodes with connections less than a preset threshold, and generates a refined state graph containing drug intervention nodes and endpoint response nodes. In this embodiment, the process of performing backdoor adjustment on the initial state graph and estimating the transfer entropy between nodes, retaining nodes with connections greater than a preset threshold as causal effect edges, removing nodes with connections less than or equal to a preset threshold, and generating a refined state graph containing drug intervention nodes and endpoint response nodes includes: Extract the drug intervention node marker and the endpoint response node marker from the initial state graph, and use them as the drug intervention node and the endpoint response node respectively. Starting from the drug intervention node and ending from the endpoint response node, traverse all directed paths from the starting point to the endpoint in the initial state graph, and extract all intermediate nodes passed through all directed paths to generate an intermediate node set.

[0054] Specifically, the drug intervention node marker and the endpoint response node marker are recorded in the attribute fields of each node in the initial state graph. The attribute fields refer to the data areas that store the node type labels. Starting from the drug intervention node and ending at the endpoint response node, a depth-first traversal algorithm is executed in the directed topology of the initial state graph to enumerate all paths that can reach the endpoint by traveling along the directed edges from the starting point. Each directed path consists of the starting point, the endpoint, and one or more intermediate nodes connected sequentially. All intermediate node names appearing in all paths are extracted, and duplicate names are removed to form a set of intermediate nodes.

[0055] Each intermediate node in the intermediate node set is taken out as a node to be tested. The conditional mutual information value of the drug intervention node state time series and the endpoint response node state time series under the given state time series of the node to be tested is calculated. The node to be tested whose conditional mutual information value is less than the preset mutual information threshold is marked as a redundant node. All redundant nodes are removed from the intermediate node set to generate a candidate adjustment variable set.

[0056] Specifically, the conditional mutual information value measures the amount of additional information that the drug intervention node's state time series can provide about the endpoint response node's state time series, given the known values ​​of the state time series of the node to be tested. The preset mutual information threshold is obtained by performing a permutation test on the state time series. This test involves randomly shuffling the observation time points of the endpoint response node's state time series one thousand times, calculating the conditional mutual information value each time, and taking the 95th percentile of all calculated values ​​as the threshold. When the conditional mutual information value is less than this threshold, it indicates that, given the node to be tested, the drug intervention node and the endpoint response node are conditionally independent, and the node to be tested does not transmit causal information. Nodes to be tested with conditional mutual information values ​​less than the preset mutual information threshold are marked as redundant nodes and removed from the intermediate node set. The remaining intermediate nodes constitute the candidate adjustment variable set.

[0057] For each pair of directed edges in the initial state graph, obtain the state time series of the starting node and the state time series of the ending node of that edge.

[0058] Specifically, a directed edge refers to a directed edge in the initial state graph that points from the starting node to the ending node. Each edge represents the direct control effect of the starting node on the ending node. For each directed edge, the state values ​​of the starting node and the ending node at all observation time points are read from the state storage area of ​​each node in the initial state graph, forming a state time sequence arranged in chronological order of observation time.

[0059] Using the state time series of each node in the candidate adjustment variable set as condition variables, the joint probability density of the state time series of the endpoint node, the starting node, and the condition variables is estimated. The ratio of the joint probability density to the marginal probability density is integrated to obtain the conditional mutual information value, which is then used as the transit entropy value.

[0060] Specifically, the transitive entropy value is used to quantify the predictive contribution of the starting node's state time series to the ending node's state time series, excluding the influence of the candidate adjustment variable set. Joint probability density estimation refers to fitting the probability density function of the high-dimensional joint distribution of the state time series of the ending node, the starting node, and all nodes in the candidate adjustment variable set using the Gaussian kernel density estimation method. The bandwidth parameter of the Gaussian kernel density estimation is determined using the Scott rule. The formula for calculating the conditional mutual information value is as follows: In the formula, Indicates that, given a set of candidate adjustment variables Given a state-time series, the starting node State time series for endpoint node The propagation entropy value of the state-time series; Indicates the starting node A continuous variable consisting of the state values ​​at all observation points; Indicates the endpoint node A continuous variable consisting of the state values ​​at all observation points; This represents a joint continuous variable consisting of the state values ​​of all nodes in the candidate adjustment variable set at all observation time points; express The joint probability density function of the three factors is obtained by fitting the kernel density estimate. for The marginal probability density function, through the analysis of... exist and Obtained by integration over the dimension; express and The joint marginal probability density function, through the analysis of... exist Obtained by integration over the dimension; express and The joint marginal probability density function, through the analysis of... exist Obtained by integration over the dimension; Represented by natural constant Logarithmic operations with base 0.

[0061] Directed edges between nodes with a propagation entropy value greater than a preset entropy threshold are retained as causal effect edges. Directed edges between nodes with a propagation entropy value less than or equal to the preset entropy threshold are removed from the initial state graph. The initial state graph after removing redundant edges is used as the refined state graph.

[0062] Specifically, the preset entropy threshold is obtained by self-sampling the transit entropy values ​​of all directed edges in the initial state graph. Self-sampling involves drawing a subsequence of equal length to the original sequence with replacement from the observation time point one thousand times, recalculating the transit entropy value of each directed edge each time, and taking the 95th percentile of all resampled values ​​as the threshold. This threshold is used to determine whether a directed edge has a statistically significant causal transitive effect. A transit entropy value greater than the preset entropy threshold indicates that the starting node's state time series has predictive information for the ending node's state time series that exceeds the level of chance, and this directed edge is retained as a causal effect edge.

[0063] If the propagation entropy value is less than or equal to a preset entropy threshold, it indicates that the directed edge between the two nodes does not have significant causal information propagation, and the directed edge is deleted from the set of directed edges in the initial state graph. After filtering all directed edges, the nodes and directed edges retained in the initial state graph constitute the refined state graph, which contains the drug intervention node, the endpoint response node, and the directed edge between them that has a significant causal effect.

[0064] The counterfactual intervention module is used to execute do operator intervention with the drug intervention node of the refined state diagram as the intervention point, and propagate the intervention effect along the causal effect edge node by node to obtain the initial counterfactual trajectory; In this embodiment, the step of executing the do operator intervention with the drug intervention node of the refined state diagram as the intervention point, propagating the intervention effect node by node along the causal effect edge to obtain the initial counterfactual trajectory includes: Using the drug intervention nodes in the refined state graph as the intervention objects, the state values ​​of the drug intervention nodes at each observation time point are replaced with preset intervention dose values. All directed edges pointing to the drug intervention nodes in the refined state graph are deleted, and a post-intervention causal graph is generated.

[0065] Specifically, the `do` operator intervention refers to actively fixing the state value of the drug intervention node in the refined state graph to a preset intervention dose value and severing the causal connections of that node controlled by all its parent nodes. The preset intervention dose value refers to the final concentration of the drug in organoid culture medium, which is derived from the half-maximal inhibitory concentration (MCC) or half-maximal effective concentration (MCD) experimental data determined in in vitro pharmacodynamic experiments. Deleting all directed edges pointing to the drug intervention node means removing all directed edges in the edge set of the refined state graph that terminate at that drug intervention node. The graph structure obtained after completing the state value replacement and directed edge deletion is the post-intervention causal graph.

[0066] Starting with the drug intervention node in the post-intervention causal graph, a topological sort is performed on the post-intervention causal graph to generate a node processing sequence. Each node is extracted sequentially according to the node processing sequence. When a node is extracted, the observed state value of the node at each observation time point is read from the refined state graph. The counterfactual state values ​​of all parent nodes of the node are used as independent variables, and the observed state values ​​are used as dependent variables. Linear regression fitting is performed on the independent and dependent variables to obtain the structural equation coefficients and intercepts of the node.

[0067] Specifically, topological sorting refers to using the Kahn algorithm to linearly sort all nodes in the directed acyclic graph (DAG) after intervention, ensuring that the starting point of each directed edge precedes its ending point, resulting in a node processing sequence. A parent node refers to the direct upstream node in the refined state graph where a directed edge points to the current node.

[0068] The counterfactual state value refers to the hypothetical state value calculated in the previous stage for each node at each observation time point under intervention conditions. For drug intervention nodes, the counterfactual state value is the preset intervention dose value. The observed state value refers to the original measurement value of the node at each observation time point stored in the refined state diagram. Linear regression fitting refers to using the least squares method to solve the regression coefficients and constant terms in the linear equation to minimize the sum of squared residuals between the regression predicted value and the observed state value.

[0069] The structural equation coefficients and intercepts are combined into a linear equation, which is then used as the structural equation for that node.

[0070] Specifically, a structural equation is a linear expression describing the causal relationship between the counterfactual state value of the current node and the counterfactual state values ​​of all its parent nodes. Combining these into a linear equation involves using regression coefficients as multiplier coefficients for each parent node, adding intercepts as constants, and organizing them into a complete algebraic equation according to a standard linear format. This linear equation serves as the structural equation for the current node, used in subsequent calculations to derive the counterfactual state value of the current node from the counterfactual state values ​​of its parent nodes.

[0071] Substitute the counterfactual state values ​​of all parent nodes at the current observation time into the structural equation of the node, calculate the counterfactual state value of the node at the current observation time, and repeat this step in the order of observation time until the counterfactual state values ​​of the node at all observation time points have been calculated.

[0072] Specifically, the current observation time point refers to a specific time position taken from all observation time points in chronological order. Substitution involves multiplying the counterfactual state values ​​of each parent node at the current observation time point by the corresponding regression coefficients in the structural equation, summing the results, and then adding the intercept. The formula for calculating the counterfactual state values ​​is as follows: In the formula, Indicates the first The node at the th Counterfactual state values ​​at each observation point in time; Indicates the first The intercept of the structural equation for each node is obtained by linear regression fitting; Indicates the first The set of indices of all parent nodes of a given node; Represents the set of parent nodes The index of a parent node in the array; Indicates the first Each node relates to its parent node. The structural equation coefficients, i.e., the parent node When the counterfactual state value changes by one unit, the first The change in the counterfactual state value of each node is obtained by linear regression fitting; Indicates the first The parent node at the _th _ ... The counterfactual state value at each observation time point. Starting from the first observation time point, the counterfactual state values ​​of each parent node are substituted sequentially to calculate the counterfactual state value of the current node. After completion, the process moves to the next observation time point until all calculations are completed.

[0073] The counterfactual state values ​​of each node at each observation time point are stored in the counterfactual state matrix according to the observation time point and node identifier, and the counterfactual state matrix is ​​used as the initial counterfactual trajectory.

[0074] Specifically, the counterfactual state matrix is ​​a two-dimensional numerical storage structure with observation time points as row indices and node identifiers as column indices. Each element in the matrix stores the counterfactual state value of the corresponding node at the corresponding observation time point. The process of storing the counterfactual state value calculated for each node at each observation time point involves filling the corresponding cell with the counterfactual state value calculated for that node at each observation time point, based on the row position of its corresponding time point and the column position of its node identifier. The complete two-dimensional matrix formed after filling in the counterfactual state values ​​for all nodes at all observation time points is the initial counterfactual trajectory.

[0075] The target screening module is used to calculate the marginal contribution of each causal effect edge from the end response node of the refined state diagram in reverse order, and select the nodes connected by the causal effect edges that are sorted in descending order to form a set of candidate targets. In this embodiment, the step of calculating the marginal contribution of each causal effect edge from the endpoint response node of the refined state graph backwards, and selecting the nodes connected by the causal effect edges that are sorted in descending order to form a set of candidate target points, includes: Using the endpoint response node in the refined state diagram as the target node, the difference between the state value of the endpoint response node at the last observation time point in the initial counterfactual trajectory and the observed state value of the endpoint response node at the last observation time point in the trajectory tensor is taken as the total effect value.

[0076] Specifically, the endpoint response node refers to the terminal biological function node marked in the refined state diagram to characterize the drug effect. The state value at the last observation time point in the initial counterfactual trajectory refers to the counterfactual state value calculated for the endpoint response node at the last observation time point under the do operator intervention condition.

[0077] The observed state value at the last observation time point in the trajectory tensor refers to the original measurement value of the endpoint response node at the last observation time point without any intervention. The total effect value is the difference between the state value at the last observation time point in the initial counterfactual trajectory and the observed state value at the last observation time point in the trajectory tensor, used to quantify the total causal effect of drug intervention on the endpoint response node.

[0078] Extract all parent nodes of the endpoint response node from the refined state diagram, remove each parent node one by one, replace the state value of the removed parent node with zero at all observation time points, and then re-execute the do operator intervention to obtain the counterfactual trajectory after the virtual removal of the parent node.

[0079] Specifically, a parent node refers to a directly upstream regulating node in the refined state graph whose directed edges point directly to the endpoint response node. These nodes exert a direct influence on the endpoint response node through causal effect edges. "Selecting each parent node one by one" means selecting one target node sequentially from the entire list of parent nodes, while keeping the original state values ​​of the remaining parent nodes unchanged. "Replacing the state value with zero" means setting the counterfactual state value of the target parent node to zero at all observation time points, simulating the biological state of complete knockout of the molecule's function. "Re-executing the do operator intervention" means, based on the replaced refined state graph, executing the do operator intervention according to the drug intervention nodes and propagating the intervention effect node by node along the causal effect edges, calculating the counterfactual state value sequence of each node at each observation time point under the condition of virtual removal of the parent node.

[0080] Extract the state value of the endpoint response node at the last observation time point in the trajectory. Use the difference between the state value of the endpoint response node at the last observation time point in the initial counterfactual trajectory and the state value after virtual removal as the marginal contribution value of the parent node to the endpoint response node.

[0081] Specifically, the trajectory refers to the counterfactual trajectory obtained by re-executing the do operator intervention after the current target parent node is virtually removed. The state value of the endpoint response node at the last observation time point refers to the counterfactual state value of the endpoint response node at the last observation time point under the virtual removal condition.

[0082] The virtual removal state value refers to the above-mentioned values, reflecting the residual effect that drug intervention can still produce after the function of the parent node is removed. The marginal contribution value refers to the difference between the state value of the endpoint response node at the last observation point in the initial counterfactual trajectory and the state value of the endpoint response node after the virtual removal of the parent node. This difference quantifies the magnitude of the independent causal effect transmitted by the parent node to the endpoint response node through the causal effect edge.

[0083] Sort all parent nodes in descending order of their marginal contribution values, extract the first preset number of parent nodes from the sorted list, and use the extracted parent nodes as the candidate target set.

[0084] Specifically, sorting by marginal contribution value from largest to smallest means ranking all parent nodes of the endpoint response node in descending order. The larger the marginal contribution value, the stronger the causal regulation contribution of the parent node to the endpoint response node. The preset number is the upper limit of the number of candidate targets to be screened based on the downstream experimental verification capability. It is usually set to five to twenty based on the throughput of the well plate that can simultaneously perform CRISPR knockout verification. Taking out the top preset number of parent nodes in the sorting means selecting the top-ranked parent nodes from the sorting results as candidate targets and summarizing these parent nodes into a candidate target set.

[0085] The target verification module is used to virtually knock out each node in the candidate target set in sequence and lock it to zero, re-execute the do operator intervention to obtain the trajectory after knockout, compare the difference in the state of the endpoint response node before and after knockout, remove the nodes whose difference is less than the preset difference threshold, and output the refined candidate target set. In this embodiment, the step of sequentially virtually knocking out and locking each node in the candidate target set to zero, re-executing the do operator intervention to obtain the knockout trajectory, comparing the difference in the state of the endpoint response nodes before and after knockout, removing nodes with a difference less than a preset difference threshold, and outputting a refined candidate target set includes: Using the candidate target set as the input set, an empty set is established as the refined candidate set, and each node in the candidate target set is taken out one by one as the current node to be processed.

[0086] Specifically, the candidate target set refers to the list of node names to be verified, sorted by marginal contribution value. The refined candidate set is an empty set used to store the node names that remain after virtual knockout verification; its initial state is empty. Retrieving each node from the candidate target set sequentially means selecting one node at a time as the current processing target according to the order in the candidate target set.

[0087] The refined state graph is copied, and the state values ​​of the current node to be processed in the copied graph are replaced with zero at all observation time points, while the state values ​​of the other nodes remain unchanged. The copied graph is then used as a virtual knockout state graph.

[0088] Specifically, copying the refined state graph refers to generating a new graph instance that is completely identical to the refined state graph in terms of node set, directed edge set, and node state values. The copying operation is used to avoid modifying the original graph structure. Replacing state values ​​with zero means resetting the counterfactual state values ​​of the currently processed node in the copied graph to zero at all observation time points, used to simulate the loss of molecular function after the gene is completely knocked out. Keeping the state values ​​of other nodes unchanged means that the counterfactual state values ​​of all other nodes except the currently processed node remain unchanged. The graph structure obtained after the replacement is the virtual knockout state graph.

[0089] Using the drug intervention node in the virtual knockout state graph as the intervention point, delete all directed edges pointing to the drug intervention node in the virtual knockout state graph, execute the do operator intervention, and obtain the knockout trajectory corresponding to the current node to be processed.

[0090] Specifically, using the drug intervention node as the intervention point means taking the drug intervention node as the target of the do operator intervention. Deleting all directed edges pointing to the drug intervention node means removing all edges in the set of directed edges in the virtual knockout state graph that terminate at the drug intervention node, making the drug intervention node a root node without a parent node. Executing the do operator intervention means, after deleting the directed edges, fixing the state value of the drug intervention node to a preset intervention dose value and calculating the counterfactual state value sequence of each node at each observation time point along the remaining causal effect edges. This sequence is the knockout trajectory corresponding to the currently processed node.

[0091] Extract the initial state value of the endpoint response node at the last observation time from the initial counterfactual trajectory, extract the knockout state value of the endpoint response node at the last observation time from the knockout trajectory, and perform a subtraction operation between the initial state value and the knockout state value to obtain the difference.

[0092] Specifically, the initial counterfactual trajectory refers to the complete sequence of counterfactual state values ​​generated by the do operator before virtual knockout. The initial state value refers to the counterfactual state value of the endpoint response node in the initial counterfactual trajectory at the last observation time point. The post-knockout state value refers to the counterfactual state value of the endpoint response node in the post-knockout trajectory corresponding to the currently pending node at the last observation time point. The subtraction operation refers to subtracting the post-knockout state value from the initial state value; the result is the decrease in the endpoint response node effect size after virtual knockout of that node.

[0093] If the difference is greater than or equal to the preset difference threshold, the node identifier of the current node to be processed is stored in the refined candidate set; if the difference is less than the preset difference threshold, the storage operation is not performed on the refined candidate set.

[0094] Specifically, the preset difference threshold is a cutoff value obtained by multiplying the standard deviation of repeated experiments between organoid culture batches by three. This threshold is used to determine whether the virtual knockout of a node causes a substantial change in the node effect size of the endpoint response. A difference value greater than or equal to the preset difference threshold indicates that the knockout of the current node leads to a decrease in the node effect size of the endpoint response that exceeds the experimental fluctuation range, and the node is retained as a valid causal target. A difference value less than the preset difference threshold indicates that the knockout of the current node does not cause a decrease in effect size that exceeds the experimental fluctuation range, and the node should not be retained.

[0095] The next node in the candidate target set is taken as the current node to be processed. This process continues until all nodes in the candidate target set have been processed. The refined candidate set is then output as the refined candidate target set.

[0096] Specifically, the next node in the candidate target set refers to the name of the next node after the current node to be processed. Processing all nodes in the candidate target set means performing all operations—copying, replacing, intervening, extracting, subtracting, and threshold comparison—on each node in the candidate target set in sequence. Outputting the refined candidate set means summarizing all node names that remain after virtual knockout verification as the final screening result.

[0097] In this embodiment, the steps of extracting the initial state value of the endpoint response node at the last observation time from the initial counterfactual trajectory, extracting the knockout state value of the endpoint response node at the last observation time from the knockout trajectory, and performing a subtraction operation between the initial state value and the knockout state value to obtain the difference include: Using the initial counterfactual trajectory as the first input matrix, the node identifier of the endpoint response node as the first row index, and the sequence number of the last observation time point as the first column index, the matrix element values ​​corresponding to the first row index and the first column index are extracted from the first input matrix and used as the initial state values.

[0098] Specifically, the first input matrix refers to a two-dimensional array of initial counterfactual state values ​​organized with observation time point numbers as rows and node identifiers as columns. The number of rows in this matrix equals the total number of observation time points, and the number of columns equals the total number of nodes. The node identifier of the endpoint response node refers to the gene symbol name of that node in the standardized causal pathway graph, which is used as the row lookup key of the matrix.

[0099] The sequence number of the last observation time point refers to the last sorting number in the entire sequence of observation time points. This number is equal to the total number of observation time points. Extracting from the first input matrix refers to reading the stored value from the cell located by both the row index and column index of the matrix. This value is the initial state value.

[0100] Using the knockout trajectory as the second input matrix, the node identifier of the endpoint response node as the second row index, and the sequence number of the last observation time point as the second column index, the matrix element values ​​corresponding to the second row index and the second column index are extracted from the second input matrix and used as the knockout state values.

[0101] Specifically, the second input matrix refers to a two-dimensional array of virtual knockout counterfactual state values ​​organized with observation time point numbers as rows and node identifiers as columns, with the same number of rows and columns as the first input matrix. The node identifier of the endpoint response node serves as the second row index, using the gene symbol name to match and locate the corresponding column position in the column header of the second input matrix.

[0102] The sequence number of the last observation time point is used as the index of the second column. This sorting number is used to match and locate the corresponding row position in the row header of the second input matrix. The extraction method in the second input matrix reads the cell value located by both the row index and the column index. This value is the state value after knocking.

[0103] The difference is obtained by subtracting the initial state value from the state value after knocking it out.

[0104] Specifically, the subtraction operation refers to taking the initial state value as the minuend and the knocked-out state value as the subtrahend, and calculating the arithmetic difference between the two. This difference is used to quantify the decrease in the endpoint response node effect after the virtual knockout of the current node to be processed.

[0105] The report generation module is used to generate a drug response mechanism explanation report based on the refined candidate target set and the initial counterfactual trajectory.

[0106] In this embodiment, generating a drug response mechanism explanation report based on the refined candidate target set and the initial counterfactual trajectory includes: Using the node identifier of each node in the refined candidate target set as the row index, the counterfactual state values ​​of each node at all observation time points are extracted from the initial counterfactual trajectory to generate the target counterfactual state matrix.

[0107] Specifically, the refined candidate target set refers to the list of node names retained after virtual knockout verification, and the node identifier refers to the gene symbol name of each node in the standardized causal pathway graph. The row index refers to the row label of the matrix using the node identifier, and the initial counterfactual trajectory refers to the two-dimensional numerical matrix storing the counterfactual state values ​​of each node at each observation time point. Extraction refers to taking the counterfactual state values ​​of each node in the refined candidate set in chronological order and writing them row by row into a new matrix, generating a target counterfactual state matrix with the refined candidate nodes as rows and the observation time points as columns.

[0108] The arithmetic mean of each observation time point is calculated column by column for the counterfactual state matrix of the target point to obtain the population average response curve for each observation time point. The standard deviation of each observation time point is calculated column by column to obtain the population response fluctuation band.

[0109] Specifically, column-wise calculation refers to performing statistical operations on the arithmetic mean and standard deviation of the column vector corresponding to each observation time point in the matrix. The population average response curve is the curve formed by connecting the arithmetic means of all observation time points in chronological order. The standard deviation is a measure of the dispersion of each element in the column relative to the arithmetic mean. The population response fluctuation band refers to the band-shaped region enclosed by adding or subtracting the standard deviation of each observation time point as the upper and lower boundaries, with the population average response curve as the center line.

[0110] Using the node identifier of each node in the refined candidate target set as an index, the transfer entropy value of each node is extracted, and this transfer entropy value is used as the causal effect strength of each target.

[0111] Specifically, the transfer entropy value refers to the directed information metric calculated during the structural refinement stage and stored in the attribute fields of each node in the refined state diagram. It quantifies the amount of directional predictive information transmitted by each node along the causal effect edge to the endpoint response node. Extraction refers to retrieving the corresponding transfer entropy value from the node attribute fields of the refined state diagram using the gene symbol name of each node in the refined candidate set as the lookup key. Causal effect strength refers to using the retrieved transfer entropy value as a quantitative indicator of the causal effect of each target point on the endpoint response node.

[0112] Divide the marginal contribution value of each target by the causal effect strength of each target, and use the quotient obtained by the division operation as the effect confidence score of each target. Arrange the targets in descending order and label the sorting number.

[0113] Specifically, the marginal contribution value refers to the magnitude of the independent causal effect of each parent node on the endpoint response node, calculated during the initial target screening stage, while the causal effect strength refers to the transmission entropy value corresponding to that target. The division operation involves using the marginal contribution value as the dividend and the causal effect strength as the divisor, calculating their quotient, which quantifies the magnitude of effect change corresponding to a unit of transmitted information. The descending order sorting refers to rearranging all targets in descending order of effect confidence score, with the sorting number being a unique number assigned to each target, starting from 1 and incrementing accordingly.

[0114] Using the observation time point as the horizontal axis and the counterfactual state value as the vertical axis, individual response curves for each target point are generated, and the population average response curve and the population response fluctuation band are superimposed on the individual response curves.

[0115] Specifically, the observation time points are plotted on the horizontal axis, with each observation time marked chronologically. The counterfactual state values ​​are plotted on the vertical axis, with the counterfactual state values ​​of each target point at each observation time point marked. The individual response curve is a two-dimensional line graph plotted for a single target point, connecting the counterfactual state values ​​of that target point at all observation time points to form a curve. In the same coordinate system, the population average response curve and the population response fluctuation band are plotted together with the individual response curve of the current target point. The population average response curve is plotted as a solid line, and the population response fluctuation band is plotted as a semi-transparent colored area.

[0116] The individual response curves, marginal contribution values, transfer entropy values, effect confidence scores, and ranking numbers of each target are combined one by one to generate drug response mechanism entries for each target.

[0117] Specifically, "target-by-target combination" refers to compiling the individual response curve, corresponding marginal contribution value, transfer entropy value, effect confidence score, and ranking number for each refined candidate target. The "drug response mechanism" entry refers to a comprehensive information unit centered on a single target, containing a visual response curve and four quantitative assessment indicators.

[0118] All drug response mechanism entries for all targets are concatenated in sequential order, and the concatenation result is output as a drug response mechanism interpretation report.

[0119] Specifically, the sorting order refers to the arrangement of the drug response mechanism entries according to the effect confidence score from high to low, and the splicing refers to the sequential connection of the drug response mechanism entries for each target in ascending order from 1 to the highest number to form continuous document content. The drug response mechanism interpretation report is a complete document generated by sequentially connecting all target entries, providing individual response characteristics and comprehensive assessment information for each candidate target.

[0120] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented in software, the above embodiments can be implemented, in whole or in part, as a computer program product. Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution.

[0121] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0122] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A digital twin prediction system for drug response in virtual organoids, characterized in that, include: The data acquisition module is used to acquire multimodal observation sequences and causal pathway maps of organoids, and to align the vectors in the multimodal longitudinal observation sequences to a unified space to obtain the trajectory tensor. The initial graph construction module is used to encode the causal pathway graph into a causal skeleton. It assigns the state values ​​of each node in the trajectory tensor and the drug intervention node labels and endpoint response node labels in the causal pathway graph to the corresponding nodes of the causal skeleton to generate the initial state graph. The structure refinement module is used to perform backdoor adjustments on the initial state graph and estimate the transfer entropy between each node. It retains the nodes with connections greater than a preset threshold as causal effect edges, removes the nodes with connections less than a preset threshold, and generates a refined state graph containing drug intervention nodes and endpoint response nodes. The counterfactual intervention module is used to execute do operator intervention with the drug intervention node of the refined state diagram as the intervention point, and propagate the intervention effect along the causal effect edge node by node to obtain the initial counterfactual trajectory; The target screening module is used to calculate the marginal contribution of each causal effect edge from the end response node of the refined state diagram in reverse order, and select the nodes connected by the causal effect edges that are sorted in descending order to form a set of candidate targets. The target verification module is used to virtually knock out each node in the candidate target set in sequence and lock it to zero, re-execute the do operator intervention to obtain the trajectory after knockout, compare the difference in the state of the endpoint response node before and after knockout, remove the nodes whose difference is less than the preset difference threshold, and output the refined candidate target set. The report generation module is used to generate a drug response mechanism explanation report based on the refined candidate target set and the initial counterfactual trajectory.

2. The digital twin prediction system for drug response based on virtual organoids according to claim 1, characterized in that, The process of obtaining the multimodal observation sequence and causal pathway map of organoids involves aligning the vectors in the multimodal longitudinal observation sequence to a unified space to obtain the trajectory tensor, including: Organoids in culture were collected using liquid chromatography-mass spectrometry, a high-content imaging system, and a microelectrode array sensor within a continuous observation window to obtain raw mass spectrometry signal streams, raw image sequences, and raw electrical impedance time series. Peak detection and retention time alignment were performed on the raw mass spectrometry signal stream to generate a molecular expression abundance matrix indexed by the observation time point; The original image sequence is subjected to three-dimensional structural reconstruction and cell nucleus segmentation to generate an organoid spatial morphology matrix indexed by the observation time point; Baseline drift correction and low-frequency filtering were performed on the original electrical impedance time series to generate an organoid barrier function matrix indexed by the observation time points; The molecular expression abundance matrix, organoid spatial morphology matrix and organoid barrier function matrix are longitudinally spliced ​​according to the observation time points to form a multimodal longitudinal observation sequence. The causal pathway graph is obtained from a public pathway database. The protein, gene and metabolite entities in the causal pathway graph are mapped to gene symbol identifiers to serve as identifiers for subsequent nodes, thereby generating a standardized causal pathway graph. Using the node order of the standardized causal pathway graph as the arrangement basis, the row vectors of the molecular expression abundance matrix, the row vectors of the organoid spatial morphology matrix, and the row vectors of the organoid barrier function matrix corresponding to each observation time point in the multimodal longitudinal observation sequence are concatenated into a unified state vector in the feature dimension. The node sequence of the unified state vector is rearranged according to the node arrangement order of the standardized causal pathway graph to generate the rearranged unified state vector for each observation time point. The unified state vectors after rearranging the observation time points are stacked in chronological order to form a three-dimensional tensor, which is the trajectory tensor.

3. The digital twin prediction system for drug response based on virtual organoids according to claim 1, characterized in that, The step of encoding the causal pathway graph into a causal skeleton, assigning the state values ​​of each node in the trajectory tensor and the drug intervention node markers and endpoint response node markers in the causal pathway graph to the corresponding nodes of the causal skeleton to generate an initial state graph includes: Extract all node names from the causal path graph and generate a node name sequence; extract the directed connections between nodes in the causal path graph and generate a set of directed edges. Construct an empty adjacency matrix using the node name sequence as row and column indices, set the main diagonal elements to zero, assign 1 to non-main diagonal elements with connections and 0 to the rest based on the set of directed edges, and use the adjacency matrix as a causal skeleton. Using the row index order of the causal skeleton as the traversal order, the node name corresponding to the row index is matched with the gene symbol identifier in the second dimension of the trajectory tensor; The third-dimensional feature channel values ​​of the trajectory tensor corresponding to the successfully matched gene symbol are expanded into a one-dimensional state sequence along the first-dimensional observation time point. The one-dimensional state sequence is written into the node storage area of ​​the causal skeleton row by row with the node name as the sequence index, generating the node state matrix of the causal skeleton. The row index of the node state matrix is ​​the observation time point, and the column index is the node name. Write the drug intervention node marker and the endpoint response node marker into the attribute field of the corresponding node in the causal skeleton, indexed by the node name; Assign each row vector of the node state matrix to the node at the corresponding row index in the causal skeleton to generate the initial state graph.

4. The digital twin prediction system for drug response based on virtual organoids according to claim 1, characterized in that, The process of performing backdoor adjustment on the initial state graph and estimating the propagation entropy between nodes, retaining nodes with connections greater than a preset threshold as causal effect edges, and removing nodes with connections less than or equal to a preset threshold, generates a refined state graph containing drug intervention nodes and endpoint response nodes, including: Extract the drug intervention node marker and the endpoint response node marker from the initial state graph, and use them as the drug intervention node and the endpoint response node respectively. Starting from the drug intervention node and ending from the endpoint response node, traverse all directed paths from the starting point to the endpoint in the initial state graph, and extract all intermediate nodes passed through all directed paths to generate an intermediate node set. Each intermediate node in the intermediate node set is taken out as the node to be tested. The conditional mutual information value of the drug intervention node state time series and the endpoint response node state time series under the given state time series of the node to be tested is calculated. The node to be tested with a conditional mutual information value less than the preset mutual information threshold is marked as a redundant node. All redundant nodes are removed from the intermediate node set to generate a candidate adjustment variable set. For each pair of directed edges in the initial state graph, obtain the state time series of the starting node and the state time series of the ending node of the edge. Using the state time series of each node in the candidate adjustment variable set as condition variables, the joint probability density of the state time series of the endpoint node, the starting node, and the condition variables is estimated, and the ratio of the joint probability density to the marginal probability density is integrated to obtain the conditional mutual information value, which is then used as the transit entropy value. Directed edges between nodes with a propagation entropy value greater than a preset entropy threshold are retained as causal effect edges. Directed edges between nodes with a propagation entropy value less than or equal to the preset entropy threshold are removed from the initial state graph. The initial state graph after removing redundant edges is used as the refined state graph.

5. A digital twin prediction system for drug response based on virtual organoids according to claim 1, characterized in that, The process involves executing the do operator intervention using the drug intervention nodes in the refined state diagram as intervention points, propagating the intervention effect node by node along the causal effect edge, and obtaining the initial counterfactual trajectory, including: Using the drug intervention nodes in the refined state graph as the intervention objects, the state values ​​of the drug intervention nodes at each observation time point are replaced with preset intervention dose values, all directed edges pointing to the drug intervention nodes in the refined state graph are deleted, and a post-intervention causal graph is generated. Starting with the drug intervention node in the post-intervention causal graph, a topological sort is performed on the post-intervention causal graph to generate a node processing sequence. Each node is extracted in sequence according to the node processing sequence. When a node is extracted, the observed state value of the node at each observation time point is read from the refined state graph. The counterfactual state values ​​of all parent nodes of the node are used as independent variables, and the observed state values ​​are used as dependent variables. Linear regression fitting is performed on the independent and dependent variables to obtain the structural equation coefficients and intercepts of the node. The structural equation coefficients and intercepts are combined into a linear equation, which is then used as the structural equation for that node. Substitute the counterfactual state values ​​of all parent nodes at the current observation time into the structural equation of the node, calculate the counterfactual state value of the node at the current observation time, and repeat this step in the order of observation time until the counterfactual state values ​​of the node at all observation time points have been calculated. The counterfactual state values ​​of each node at each observation time point are stored in the counterfactual state matrix according to the observation time point and node identifier, and the counterfactual state matrix is ​​used as the initial counterfactual trajectory.

6. A digital twin prediction system for drug response based on virtual organoids according to claim 1, characterized in that, The step of calculating the marginal contribution of each causal effect edge from the endpoint response node of the refined state graph backwards, and selecting the nodes connected by a preset number of causal effect edges in descending order as the candidate target set, includes: Using the endpoint response node in the refined state diagram as the target node, the difference between the state value of the endpoint response node at the last observation time point in the initial counterfactual trajectory and the observed state value of the endpoint response node at the last observation time point in the trajectory tensor is taken as the total effect value. Extract all parent nodes of the endpoint response node from the refined state diagram, remove each parent node one by one, replace the state value of the removed parent node with zero at all observation time points, and then re-execute the do operator intervention to obtain the counterfactual trajectory after the virtual removal of the parent node. Extract the state value of the endpoint response node at the last observation time point in the trajectory, and take the difference between the state value of the endpoint response node at the last observation time point in the initial counterfactual trajectory and the state value after virtual removal as the marginal contribution value of the parent node to the endpoint response node. Sort all parent nodes in descending order of their marginal contribution values, extract the first preset number of parent nodes from the sorted list, and use the extracted parent nodes as the candidate target set.

7. A digital twin prediction system for drug response based on virtual organoids according to claim 1, characterized in that, The process involves sequentially and virtually knocking out each node in the candidate target set and locking it to zero, then re-executing the do operator to obtain the trajectory after knockout. The difference in the state of the endpoint response nodes before and after knockout is compared, and nodes with a difference less than a preset difference threshold are removed. The refined candidate target set is then output, including: Using the candidate target set as the input set, an empty set is established as the refined candidate set, and each node in the candidate target set is taken out as the current node to be processed. The refined state graph is copied, and the state values ​​of the current node to be processed in the copied graph are replaced with zero at all observation time points, while the state values ​​of the other nodes remain unchanged. The copied graph is used as a virtual knockout state graph. Using the drug intervention node in the virtual knockout state graph as the intervention point, delete all directed edges pointing to the drug intervention node in the virtual knockout state graph, execute the do operator intervention, and obtain the knockout trajectory corresponding to the current node to be processed. Extract the initial state value of the endpoint response node at the last observation time from the initial counterfactual trajectory, extract the knockout state value of the endpoint response node at the last observation time from the knockout trajectory, and perform a subtraction operation between the initial state value and the knockout state value to obtain the difference; If the difference is greater than or equal to the preset difference threshold, the node identifier of the current node to be processed is stored in the refined candidate set; if the difference is less than the preset difference threshold, the storage operation is not performed on the refined candidate set. The next node in the candidate target set is taken as the current node to be processed. This process continues until all nodes in the candidate target set have been processed. The refined candidate set is then output as the refined candidate target set.

8. A digital twin prediction system for drug response in virtual organoids according to claim 7, characterized in that, The process of extracting the initial state value of the endpoint response node at the last observation time from the initial counterfactual trajectory, extracting the knockout state value of the endpoint response node at the last observation time from the knockout trajectory, and performing a subtraction operation between the initial state value and the knockout state value to obtain the difference includes: Using the initial counterfactual trajectory as the first input matrix, the node identifier of the endpoint response node as the first row index, and the sequence number of the last observation time point as the first column index, the matrix element values ​​corresponding to the first row index and the first column index are extracted from the first input matrix as the initial state values. Using the knockout trajectory as the second input matrix, the node identifier of the endpoint response node as the second row index, and the sequence number of the last observation time point as the second column index, the matrix element values ​​corresponding to the second row index and the second column index are extracted from the second input matrix and used as the knockout state values. The difference is obtained by subtracting the initial state value from the state value after knocking it out.

9. A digital twin prediction system for drug response in virtual organoids according to claim 8, characterized in that, The drug response mechanism explanation report generated based on the refined candidate target set and the initial counterfactual trajectory includes: Using the node identifier of each node in the refined candidate target set as the row index, the counterfactual state values ​​of each node at all observation time points are extracted from the initial counterfactual trajectory to generate the target counterfactual state matrix. The arithmetic mean of each observation time point is calculated column by column for the counterfactual state matrix of the target point to obtain the population average response curve at each observation time point. The standard deviation of each observation time point is calculated column by column to obtain the population response fluctuation band. Using the node identifier of each node in the refined candidate target set as an index, the transfer entropy value of each node is extracted, and the transfer entropy value is used as the causal effect strength of each target. Divide the marginal contribution value of each target by the causal effect strength of each target, and use the quotient obtained by the division operation as the effect confidence score of each target. Arrange the targets in descending order and mark the sorting number. With the observation time point as the horizontal axis and the counterfactual state value as the vertical axis, an individual response curve for each target point is generated, and the population average response curve and the population response fluctuation band are superimposed on the individual response curve. The individual response curves, marginal contribution values, transfer entropy values, effect confidence scores, and ranking numbers of each target are combined one by one to generate drug response mechanism entries for each target. All drug response mechanism entries for all targets are concatenated in sequential order, and the concatenation result is output as a drug response mechanism interpretation report.

Images