Time-stage driven pathfinding framework for optimized medical treatments

Abstract

Coming up with the correct patient-specific treatment procedure is a non-trivial task. One of the main challenges is that the problem is combinatorial, and each treatment can be followed by numerous treatments. In real-life scenarios, one cannot experiment with all probable treatment combinations that may provide the necessary positive outcome. In addition, the task of correct drug prescription is challenging because, at different time-stages, the prescription may provide different results, a concept not widely explored in the literature on computational modeling. To address this task, we model the problem as a search problem and propose two algorithms that construct a directed acyclic graph (DAG), a directed cyclic graph, and a tree for each patient to explore various treatment combinations in an iterative and recursive manner. Each patient has three corresponding datasets, one for each stage, representing the features the patient has demonstrated during the recovery process. As a result, we provide a framework for identifying treatment options that may not have been explored previously while incorporating the concept of time-stage-based observations in the search procedure as novel contributions to the existing literature. We provide evaluation and scaling methods and identify current limitations and future research directions.

Keywords:

• networks
• graph search
• drug-treatment optimization
• multi-stage datasets
• medical software systems

1. Introduction

The problem that we try to address in this work is closely related to the task of iterative target class allocation proposed by (Aslanyan et al., 2022), where the goal is to move the object to the marked class with successive classification steps. This task can be addressed using machine learning-based techniques, such as through cascade classifiers; however, we model it using a graph-based approach. For solving this problem for the medical treatment domain from the perspective of graph search, where the target class is the Home discharge outcome, we obtain and use historical treatment and charted patient observations as a result of MIMIC-III clinical database processing.

Although not at a larger scale, searching algorithms have previously been applied to drug discovery tasks. Based on our conducted literature review, such applications imply a one-stage search process and do not extend to multiple stages, which brings additional complexities to the task of finding the optimal drug treatment prescription procedure. In this work, each node in the graph represents a patient who has a corresponding feature vector, and the graph itself describes a patient–feature interaction. There are similar works in the literature that are based on protein–protein interaction networks which search for a target node in the graph. We construct a graph and perform a search not only to identify desirable nodes (features) but also to derive a treatment strategy based on the medications of the patient using feature similarity.

As a result of our experiments, each patient may have zero, one, or many treatment possibilities depending on the algorithmic approach. After constructing the graphs, we search for the optimal solutions using existing pathfinding algorithms, which are needed in case multiple shortest paths exist, propose cost and heuristic functions, and show the optimality for one of the functions in the scope of the current graph search domain. In addition, we also provide an API-driven software pipeline for the user to obtain a graph and a treatment path for the requested features and diagnosis. The paper’s contributions can be summarized in the following points.¹

• We provide two algorithms that, based on stage-based logic, construct a cyclic graph, a DAG, and a tree for each patient based on the admission features by dynamically searching for the optimal treatment paths. One of the algorithms can have both deterministic and non-deterministic behavior. The algorithms are applied to datasets that group patients’ recovery into three stages, namely initial, intermediate, and final, which include charted patient observations (features) and drugs given to the patient for a given stage.

• We provide heuristic and cost functions devised for performing the search on the constructed graphs for the medical treatment-search task.

• We provide a clustering-based evaluation of the DAG construction experiment to examine its effectiveness. Similar to the approach described in (Chen et al., 2016), we also provide betweenness and cosine similarity-based evaluations of both experiments to identify the most important nodes and observe their similarities to the target features.

• We provide a serverless pipeline based on the modern backend and cloud technologies as a scaling platform for algorithmic results.

2. Related work

Network-based approaches are not uncommon in developing drug combination therapies and candidate-drug prediction tasks. Rivas-Barragan et al. (2020) present Drug2ways, a novel methodology that leverages multi-modal causal networks to predict drug candidates. They implement an algorithm that reasons over causal paths in large-scale biological networks to propose drug candidates for a given disease. The algorithm traverses the ensemble of paths between a drug and disease to propose the drugs that are likely to cure the disease based on the information contained in the network. Directionality inferred through the paths can be used to identify drug candidates, and the algorithm can also find ways of using existing drugs and propose combination therapies.

Many knowledge-graph methods have been developed to identify new drug therapies for diseases. Most such methods are based on counting the number of biomedical concepts, such as drugs, diseases, and genes in the indirect paths defined as sequences of two triples between two drugs, two diseases, or between a drug and a disease. The assumption is that a high number of intermediate concepts indicate similarity between drugs (Malas et al., 2019). Successful network-based frameworks have also been developed for accelerating drug discovery and quantifying disease–disease and drug–disease relationships (Cheng et al., 2018; Guney et al., 2016; Menche et al., 2015). These approaches have allowed moving from “one-drug, one-target” paradigm towards “multiple-drugs, multiple-targets” possibilities (Hopkins, 2008; Vidal et al., 2011; Yıldırım et al., 2007).

(Chen et al., 2016) apply a shortest path algorithm to identify candidate genes in the protein–protein network, and the analysis shows that some of them are deemed breast cancer-related genes according to the most recently published literature. In this work, the model is based on the information on protein–protein interactions, and such methods have been previously applied to identify genes of other diseases, such as age-related macular degeneration (Zhang et al., 2013) and hepatocellular carcinoma (Jiang et al., 2013). The weighted network is constructed by using proteins as nodes, and two nodes are adjacent if the corresponding proteins can interact with each other. The score of the interaction is defined as $\mathit{Q}(p_1,p_2)$, which ranges between 150 and 950, where $p_1$ and $p_2$ are corresponding proteins of nodes $v_1$ and $v_2$. The final weight is given by 1000-$\mathit{Q}(p_1,p_2)$. The shortest paths are constructed for any two known breast cancer-related genes using Dijkstra’s algorithm, and for each node gene in the network, the number of times the node/gene appears as an intermediary node is counted, and this value is referred to as betweenness in the study. To avoid false candidate-gene discoveries, random gene-set selection and permutation FDR calculations are performed, and the genes with small FDRs are picked as significant breast-cancer genes. Li et al. (2013) propose a computational method and apply Dijkstra’s algorithm for identifying cell lung cancer (SCLC) and non-small cell lung cancer (NSCLC) related genes with a shortest path approach in the protein–protein interaction (PPI) network. Based on data from the STRING database resource, the authors show that some of the genes identified from shortest paths have been reported to be related to lung cancer, and the candidate genes identified possessed more functional similarity with known cancer genes than those identified from gene expression profiles. The shortest path algorithm has also been applied to identify Gastric Cancer-related genes in PPI networks by Jiang et al. (2014) and thyroid cancer-related genes by Jiang et al. (2015).

The task of finding the optimal drug treatment procedures is also related to the research area of constructing drug recommendation systems. Health recommender systems (HRS) can make early diagnoses, predict disease progression, and make appropriate recommendations according to the health status of the patients. Some HRSs help support medical treatment and prognosis, while others offer medicine based on patient reviews using sentiment analysis and feature engineering (Bhimavarapu et al., 2022). There are also examples of graph-based approaches for providing personalized medical recommendations. Bhoi et al. (2021) propose a two-stage recommender system comprising attention-based RNNs and graph networks to model drug co-occurrences in the EHR. They also adapt GAT to incorporate the varying importance of drug interactions for learning effective drug embeddings and justify a particular medication by providing contribution percentages among diagnoses and procedures. Similar to this work, MIMIC-III is one of the datasets used for the experiment. An existing recommender system, such as LEAP (Zhang et al., 2017), uses patients’ current visits and drug-to-drug interactions to predict a list of medications; however, patients’ previous visits are not considered. Another example is the recommender system DMNC (Le et al., 2018) which uses temporal differences in patients’ past visits when predicting medications. Mao et al. (2022) present a medical intelligent system called MEDGCN that can automatically recommend patients’ medications based on their incomplete lab tests and estimate lab values that have not yet been taken using graph convolutional networks (GCN). The experiments are validated on NMEDW and MIMIC-III datasets. Some other examples of drug recommendations using graph-based methods are GAMENet by (Shang et al., 2019), which integrates drug–drug interaction knowledge for providing safe and personalized medical recommendations using graph neural networks, and G-BERT by (Shang et al., 2019), which integrates GNNS with transformer-based encoders for solving a medical recommendation task, achieving state-of-the-art performance.

Finding an optimal path with various quality measures is often regarded as an NP-hard problem, even for low-dimensional spaces (Canny & Reif, 1987). Most common deterministic and sampling-based algorithms aim to capture topological space connectivity using a discrete graph structure, and the final solution is extracted by searching the graph using Dijkstra’s algorithm or one of its variants. Overall, there exists a challenge in designing high-quality paths using a representative graph structure, as stated in Raveh et al. (2011), who introduce an algorithm for merging an arbitrary number of input motion paths into a hybrid output path, which is based on the observation that the quality of certain sub-paths within each solution may be higher than the quality of the entire path. To optimize the combinations of therapeutic interventions and better manage the clinical outcomes of complex diseases, Calzolari et al. (2008), propose a novel application of searching algorithms. Their findings demonstrate the effectiveness of the algorithms in identifying optimal combinations compared to random search and suggest that modified searching algorithms have the potential to enhance the discovery of novel therapeutic drug combinations. Calzolari et al. (2008) also state that with searching strategies for combined interventions, there is no need for an experimental design where possible combinations of drugs for all selected doses are tested.

There are also approaches aimed at finding new drugs or a combination of drugs that can be used in combination therapies. Combination therapies are thus widely used in the treatment of complex diseases, from hypertension to cancer and infectious diseases. The systematic identification of drug combinations that offer high efficacy and low toxicity is mostly driven by intuition rather than by established principles. In their work, to tackle such challenges (Cheng et al., 2019), come up with a network-based drug combination strategy to quantify the relationship between drug targets and disease proteins.

The availability and scalability of medical software systems that can computationally serve the methodological findings in the medical/healthcare domain are as much important as the findings themselves. Gribova et al. (2019) propose a concept of a cloud environment that can be used for the development and use of medical intelligent systems. Software systems such as BioXcel Therapeutics, BERG, XtalPi’s ID4, Deep Genomics, IBM’s Watson, and Google’s DeepMind Health can reduce the mortality rate in complicated medical diagnoses based on clinical data (Rajan & Paranthaman, 2022).

3. Materials and methods

3.1. Data

There are datasets for three stages for each diagnosis. The first stage is the initial stage, the second is the intermediate stage, and the third stage is the final stage. Each stage has two datasets. The first dataset has data about the historical patient measurements for 10 features. We label them as the Features dataset. The second dataset has data about the historical list of drugs for that stage, and we label them as the Drugs dataset. Such datasets are constructed for four diagnoses, which are PNEUMONIA, CONGESTIVE HEART FAILURE, SEPSIS, and DIABETIC KETOACIDOSIS. The drugs are the procedures or the graph edges, which tell which procedure should be given to the patient who is at a certain state with the given features. The datasets for stages are constructed based on the concept of time. Features datasets contain averaged information about 10 distinct features for $\mathit{t}$ days for each stage. An example can be seen from Table 1. The features are O2 saturation pulse oximetry, Heart Rate, Respiratory Rate, Noninvasive Blood Pressure mean, Noninvasive Blood Pressure systolic, Noninvasive Blood Pressure diastolic, Temperature Fahrenheit, Arterial Blood Pressure systolic, Arterial Blood Pressure diastolic, Arterial Blood Pressure mean. $\mathit{t}$ is 2 for stages 1 and 2 for all patients, and 2 $\le \mathit{t}\le$ 4 for stage 3. This variation in $\mathit{t}$ allows to increase the number of observations for the experiment. Similar to the Features datasets, the Drugs datasets contain drug prescriptions based on the same logical division of $\mathit{t}$ described above. Possible discharge locations are ICF, HOME, DEAD/EXPIRED, SHORT TERM HOSPITAL, HOME HEALTH CARE, etc. If the location is HOME, we encode it as 1; otherwise, the encoding is 0. This indicates that there is only one positive discharge, and the optimal path should lead to a location of this type. All the datasets are derived from pre-processing four database tables from the MIMIC-III clinical database, a freely accessible database comprising de-identified health-related data from over 40,000 patients (Johnson et al., 2016). The tables that have been used are ADMISSIONS, D_Items, and PRESCRIPTIONS datasets and a subset of the CHARTEVETS dataset.

Table 1. Features dataset example for a single patient in a stage before averaging

Timestep	Feature 1	Feature 2	$\cdots$	Feature 10
$t_1$	value 1	value 1	$\cdots$	value 1
$t_2$	value 2	value 2	$\cdots$	value 2
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
$t_{10}$	value 10	value 10	$\cdots$	value 10

3.2. DAG construction experiment

We first model the target class allocation problem as a directed acyclic graph (DAG) construction and search problem using existing pathfinding algorithms adapted to our multi-stage treatment logic. Figure 1(1) shows that the number of time-stages in our datasets corresponds to the depth of the tree and is 3. Each node or state represents the health state of the patient represented by a 1-dimensional feature vector having 10 numerical observations, one for each feature. Actions are represented by drug sequences that move the patient from state $\mathit{a}$ at time-stage $\mathit{i}$ to state $\mathit{a}{ }^{\mathrm{\prime }}$ at time-stage $\mathit{i}+1$. Each state can have at most five child nodes; in other words, the set of actions for each state is limited to five by default but can vary as an input parameter. How the algorithm responsible for constructing the state space works is described in detail in the following steps. Each time-stage has its historical list of patients, features, and treatments that have been prescribed in each stage. The data for these three stages is denoted as the “training” set. For the first stage, in addition, we also possess information about patients and their features who have just been admitted, and regard these observations are part of the “test” set. We construct one directed acyclic graph (DAG) for each patient. More specifically, we iterate over each “test” instance, which possesses the features at the first stage, and also iterate over the available “training” instances for the first stage. Our step cost function $\mathit{C}$ is based on an RMSE measure (1), which takes two nodes and calculates the square root of the sum of squared distance between the feature vectors $\mathit{F}$ of these nodes. The ideal RMSE score will have a score of 0. We add 1 to make the concept of the path more intuitive, making the minimum path cost to be 1. The first five “training” instances with the lowest cost values become the child nodes or the next possible states for the 2nd stage. The paths (edges) that connect the “testing” instance with its successor states are the drug-treatment actions. This procedure is repeated for each of the obtained five child nodes, finding the most similar child nodes in terms of features from 2nd stage “training” features, then is repeated for the resulting child nodes, connecting them to the states at the final 3rd stage. This process is conducted for each “testing” instance. What makes this process unique is that all the states available at the final time-stage are the instances that have been positively discharged (discharged to Home). The $\mathit{target}$ is a state selected from such states, and this approach makes sure the final selected path will lead to a state having the features of a patient with a successful historical discharge. If the final node has a successful discharge, at least conceptually, each child node from each time-stage will at least be as good as its parent node, given they all are part of a successful drug treatment path. Each drug will have a target node, meaning there is at least one possible solution path for each graph. The intuition is to construct a stage-progression-based drug treatment graph based on similar but not identical historical treatment information for that stage. By default, the search happens using Dijkstra’s algorithm, which for a given vertex in the graph finds the path with the lowest cost between the vertex and any other vertex. It is also often used for finding the shortest path from a single vertex to a single destination vertex (Konakalla, 2014), in this case, the target node. The algorithm, however, is uninformed as the target node is not taken into consideration until the node obtained from path expansion is found to be the target. Search can become informed by taking into account the value of some heuristic function, which can say some nodes are more promising and will lead to a solution faster. More formally, a heuristic function $\mathit{h}(\mathit{n})$, where $\mathit{n}$ is a node, is said to be admissible if $\mathit{h}(\mathit{n})$ is always less than or equal to the actual cost of a lowest-cost path from node $\mathit{n}$ to a target (David Poole & Mackworth, 2017). To identify the best path using heuristic information, we offer two heuristic functions and do the search using the A Star algorithm, a variant of Dijkstra’s algorithm, which takes into account heuristic information. The first heuristic function (2) is analogous to a “straight-line” heuristic, which for a node $\mathit{n}$ directly measures the $\mathit{RMSE}$ score between node $\mathit{n}{ }^{\mathrm{\prime }}\mathit{s}$ features $\mathit{F}$ and the target node’s features. Note that for all nodes located $n_{\mathit{d}-1}$ step away from the target, the cost function $\mathit{C}(\mathit{n})$ and heuristic function $\mathit{h}(\mathit{n})$ will be equal, which still satisfies the admissibility condition. The heuristic is admissible in the scope of this experiment but may not be admissible for other datasets. We also provide another heuristic function (3) and a proof of admissibility in (3.2). $\mathit{kstest}$ is the default function name from the SciPy library, performing a Kolmogorov–Smirnov test for goodness of fit (Virtanen et al., 2020) and providing a $\mathit{P}$ value, where if the $\mathit{P}$ value is bigger than 0.05, we accept the null hypothesis that the distributions are identical. The intuition is that the estimated cost from node $\mathit{n}$ to the target node is the node’s depth weighted by the similarity to the target node, which for the given cost function, is always admissible. The algorithm is provided in the Appendix A as Algorithm (1). In the cost and heuristic function formulas, we use the terms RMSE and kstest to avoid extra formula complications and demonstrate that the names correspond to function names from external libraries.

Figure 1. DAG visualizations for the diagnoses.

(1) $\mathit{C}(\mathit{n},\mathit{m})=1+\mathit{RMSE}(F_n,F_m)$(1)

(2) $\mathit{h}(\mathit{n},\mathit{target})=1+\mathit{RMSE}(F_n,F_{\mathit{target}})$(2)

(3) $\mathit{h}(\mathit{n},\mathit{target})=1-\mathit{kstest}(F_n,F_{\mathit{target}})$(3)

(4) $\mathit{C}(n_{\mathit{d}-\mathit{k}},\mathit{target})=\mathit{k}\ast (1+\mathit{RMSE}(F_n,F_{\mathit{target}}))$(4)

(5) $\mathit{h}(n_{\mathit{d}-\mathit{k}},\mathit{target})=\mathit{k}\ast (1-\mathit{kstest}(F_n,F_{\mathit{target}}))$(5)

Definition 3.1. Suppose $\mathit{G}$ is a DAG representing the multi-stage drug treatment state space. $\mathit{start}$ is the initial state, $\mathit{T}(s_i,\mathit{a})$ is the transition model that returns state $\mathit{s}{{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{i}+1}$ belonging to stage $\mathit{i}+1$ after taking action $\mathit{a}$ at state $\mathit{s}$ from stage $\mathit{i}$. $\mathit{d}$ denotes the maximum depth $\mathit{d}$ which is equal to the number of treatment stages, and $\mathit{target}$ denotes the end state located at depth $\mathit{d}$.

Definition 3.2. Let $\mathit{C}(n_{\mathit{d}-\mathit{k}},\mathit{target})$ denote the cost of the shortest path from node $\mathit{n}$ at depth $\mathit{d}-\mathit{k}$ to the target, and $\mathit{h}(n_{\mathit{d}-\mathit{k}},\mathit{target})$ denote the heuristic cost of $\mathit{n}$ to the target. $\mathit{C}(n_{\mathit{d}-\mathit{k}},\mathit{target})$ is given by equation (4)(4) $\mathit{C}(n_{\mathit{d}-\mathit{k}},\mathit{target})=\mathit{k}\ast (1+\mathit{RMSE}(F_n,F_{\mathit{target}}))$(4) , and $\mathit{h}(n_{\mathit{d}-\mathit{k}},\mathit{target})$ is given by equation (5)(5) $\mathit{h}(n_{\mathit{d}-\mathit{k}},\mathit{target})=\mathit{k}\ast (1-\mathit{kstest}(F_n,F_{\mathit{target}}))$(5) , where $F_n$ and $F_{\mathit{target}}$ are the feature vectors of node $\mathit{n}$ and the node $\mathit{target}$. The range of $\mathit{C}(n_{\mathit{d}-\mathit{k}},\mathit{target})$ is $R_{\mathit{C}({\mathit{n}}_{\mathit{d}-\mathit{k}},\mathit{target})}:\mathit{k}\ast [1,1024]$, and the range of $R_{\mathit{h}({\mathit{n}}_{\mathit{d}-\mathit{k}},\mathit{target})}:\mathit{k}\ast [0,1]$.

Proposition 3.3. Given a DAG $\mathit{G}$, $\mathit{h}(n_{\mathit{d}-\mathit{k}},\mathit{target})\le \mathit{C}(n_{\mathit{d}-\mathit{k}},\mathit{target})$ for any node $\mathit{n}$ $\mathit{k}$ steps away from the maximum depth $\mathit{d}$.

In other words, the theorem states that the heuristic is admissible for a given node at a given depth of the graph. The admissibility is visible from the definition, but we also provide a formal proof.

Proof. Base Case: When $\mathit{k}=1$, $\mathit{h}(n_{\mathit{d}-1},\mathit{target})=1-\mathit{kstest}(F_n,F_{\mathit{target}})$ and $\mathit{C}(n_{\mathit{d}-1},\mathit{target})=1+\mathit{RMSE}(F_n,F_{\mathit{target}})$. $\mathit{h}(n_{\mathit{d}-1},\mathit{target})\le \mathit{C}(n_{\mathit{d}-1},\mathit{target})$ as the range of $\mathit{h}(n_{\mathit{d}-1},\mathit{target})$ is always smaller than the range of $\mathit{C}(n_{\mathit{d}-1},\mathit{target})$.

Inductive Case: Suppose $\mathit{h}(n_{\mathit{d}-\mathit{k}},\mathit{target})\le \mathit{C}(n_{\mathit{d}-\mathit{k}},\mathit{target})$, then $\mathit{h}(n_{\mathit{d}-(\mathit{k}+1)},\mathit{target})$=$(\mathit{k}+1)\ast (1-\mathit{kstest}(F_n,F_{\mathit{target}}))$ and $\mathit{C}(n_{\mathit{d}-(\mathit{k}+1)},\mathit{target})=(\mathit{k}+1)\ast (1+\mathit{RMSE}(F_n,F_{\mathit{target}})).$ The expressions reduce to the base case if we divide both sides by $\mathit{k}+1$, which proves the induction.

3.3. Graph and tree construction experiment

Although in this experiment we construct a cyclic graph and a tree, which again are graphs, it differs from the DAG experiment (3.2) in multiple ways. First, the stages are not fixed, a target state may exist at depth 0 or depths deeper than 3, 2 new cost functions, in addition to cost function (1), are provided. Here, every cost function output is compared to a threshold value to check whether the algorithm should terminate or not. The threshold function for cost function (1) is given by (6), which calculates the average RMSE across the rows of the target feature matrix $\mathit{T}$, which is the third-time features of those patients who have been discharged to Home, $\mathit{n}$ denoting the number of rows, $\mathit{T}[\mathit{i},:],\mathit{T}[\mathit{j},:]$, denoting rows of the matrix.

(6) $\mathit{threshold}=1+\frac{1}{n^2}\sum _i^n\sum _j^n\mathit{RMSE}(\mathit{T}[\mathit{i},:],\mathit{T}[\mathit{j},:])$(6)

1 is added to prevent the threshold from being 0. The algorithm can have non-deterministic behavior as part of the imposed complexity, and finally, the algorithm is recursive. The path can have no solution, acceptable solution, or a goal solution, the details of which can be found in the Appendix A as part of Algorithm (2). The algorithm of the experiment can search for the optimal drug path by constructing a graph, where directed edges but not self-loops are allowed, or a tree, where undirected edges are not allowed. Similar to the breadth-first search traversal, the algorithm uses a queue for node expansion where the most similar nodes are added to the queue first. The graph can have multiple shortest paths, while the tree will always have a single shortest path given there exists a target solution node. By default, each node can have three child nodes, but the number can vary as an input parameter to the algorithm. Furthermore, both the graph and the tree can be searched in a deterministic and non-deterministic manner. Deterministic search means that the child nodes will have features that belong to a dataset from one of the three stages. The idea is that the existing features as a result of existing procedures are a good approximation, representing the mean of all the potential features that could be obtained as a result of procedure application. Meanwhile, non-deterministic means we do allow features to vary which implies deviating the features with a certain probability from the observed to include some randomness and bring it closer to the real-life scenario. We experiment with two additional cost functions and also provide threshold functions for our recursive algorithm to terminate. The cost function (7) shows the average combined deviation for the features of node $\mathit{n}$, which in the algorithm is compared with the threshold denoting the sum of standard deviations obtained from target features from the third time-stage (8). We add 1 to the cost and threshold functions for the same reason as described in the case of cost function (1).

(7) $\mathit{C}(\mathit{n})=1+\sum _{\mathit{i}=1}^{10}\sqrt{{(F_n[\mathit{i}]-\mathit{A}[\mathit{i}])}^2}$(7)

where $\mathit{F}$ and $\mathit{A}$ are 1D vectors with 10 observations, representing the feature vector and the target feature vector averaged across rows, respectively. Note that here $\mathit{C}(\mathit{n})$ depends only on a single input, as $\mathit{A}$ is the same for each $\mathit{n}$.

(8) $\mathit{threshold}=1+\sum _{\mathit{i}=1}^{10}\mathit{\sigma }(\mathit{T}[:,\mathit{i}])$(8)

where $\mathit{T}$ is the target feature matrix, $\mathit{T}[:,\mathit{i}]$ is a feature column in the matrix. The intuition is that if the node’s total feature deviations are within the acceptable deviation range, then the algorithm can stop the execution.

The second cost function given by (9) is a vector-based subtraction operation, counting the number of features that are in the acceptable deviation range, and the threshold is a constant of 2, which is a varying input parameter to the algorithm and is not explicitly given. The intuition is that if the patient does not demonstrate acceptable values for the majority of the features, then the search should continue.

(9) $\mathit{C}(\mathit{n})=1+(10-\mathit{count}(\sqrt{F_n-A^2}<\mathit{\sigma }(\mathit{T})))$(9)

4. Results

By looking at Figure 1(1), you can see the DAG visualizations for three of the diagnoses. The concept of time-stage is incorporated in subfigure (3.1) but applies to all Figure 1(1’s) subfigures. The nodes at each stage show the features that the patient can have after being given a drug procedure. Each such DAG has a solution located at depth 3. Meanwhile, by looking at Figure 2 and Table 2, the latter showing sample paths per patient using a given experimental approach, you see a tree constructed with a deterministic approach and what a drug treatment path looks like for a given diagnosis. If a path exists, it either has found an acceptable solution or a goal solution. The tree, described in the table, has no solution, as it reaches a depth of 5, and a target node still has not been found, so a base case is met, and the algorithm terminates without giving a solution. Overall the paths are obtained for different cost functions and graph types, but only some of the combinations are provided in Table 2. From the table, “start:0” means the path is for the patient with id “0”, “start:4” is for the patient with id 4, etc. Figures 3 and 4 are examples of a non-deterministic tree and a deterministic graph, respectively, and the search paths can be found in Table 2. The paths can be explained as follows: if you look at the first observation from Table 2, you see that the path is “start:0” $\to$ “t1:5”$\to$“t2:35”, $\to$ “t3:9”. This means after the initial patient features are observed, the patient should be given the medication of the patient with id 5 from time-stage 1, after that, the medication of patient with id 35 from time-stage 2, then the medication of patient 9 in time-stage 3 for a successful discharge. Note that a non-deterministic path has some extra alpha-numeric letters attached to the path, such as “t2:14:0.88:jKp”. This is done to assign the node a unique id and make sure no two nodes share the same label as a result of a non-deterministic action. You can also see that for all DAGs, the paths include one node from each stage, meaning one procedure is assigned relevant to the stage the patient currently is in. For the tree and the graph, the path search is more complex, and if the patient shares similar features with patients from later time-stages, such medications can be given disregarding the stage progression. An interesting example is the result of the non-deterministic tree obtained for SEPSIS diagnosis for the patient with id “0”, where all medications are chosen from the second time-stage before the algorithm suggests that a solution has been found. The corresponding medications are not listed, but each patient for each time-stage has a drug medication sequence from the Drugs dataset. Although the real effectiveness of the procedure can only be verified with actual hospital-level medication, this method of searching for the best procedure computationally is expected to yield the desired results.

Figure 2. Deterministic tree.

Figure 3. Non-deterministic tree.

Figure 4. Deterministic graph.

Table 2. Shortest path examples

Search	Diagnosis	Graph Type	Cost Method	Path
A Star	PNEUMONIA	DAG	1	‘start:0’ $\to$ ‘t1:5’$\to$
				‘t2:35’ $\to$ ‘t3:9’
A Star	CONGESTIVE HEART FAILURE	DAG	1	‘start:3’ $\to$ ‘t1:0’ $\to$
				‘t2:0’ $\to$ ‘t3:0’
A Star	DIABETIC KETOACIDOSIS	DAG	1	‘start:4’ $\to$‘t1:4’ $\to$
				‘t2:5’ $\to$ ‘t3:2’
Dijkstra	PNEUMONIA	Tree	3	no solution path
		Deterministic
Dijkstra	SEPSIS	Tree	3	‘start:0’ $\to$ ‘t2:2:0.84:cAV’ $\to$
		Non-Deterministic		‘t2:18’ $\to$ ‘t2:7’ $\to$ ‘t2:14:0.88:jKp’
Dijkstra	DIABETIC	Graph
	KETOACIDOSIS	Deterministic	2	‘start:1’ $\to$ ‘t1:5’ $\to$ ‘t2:5’

4.1. Evaluation

4.1.1. K-Means

We provide a clustering-based approach for evaluating the performance of the DAG construction algorithm. K-Means algorithm with three clusters is fit on the features datasets of each stage, assigning each unique instance to a given cluster. The algorithm is fit on whole datasets, meaning we obtain clusters for the “testing” instances of each diagnosis. We then obtain the treatment suggestions from the graph search for each of the three stages. If the suggested path (treatment) of the patient belongs to the same cluster as the patient’s actual treatment cluster, then the algorithm has been successful in performing the search. By looking at Table 3, we see that in all cases, the algorithm with above 60% accuracy assigns the correct medication cluster. We take these results as a baseline for improvement. We also note that the results are subject to small variations for different runs as the clustering algorithm is fit during execution. Changing the number of clusters may also affect the results.

Table 3. Evaluation results

Diagnosis	Results
PNEUMONIA	72.22%
CONGESTIVE HEART FAILURE	61.11%
SEPSIS	61.11%
DIABETIC KETOACIDOSIS	80.0 %

4.1.2. Betweenness and cosine similarity

We calculate how many times each node appears as part of a solution path for each experiment, and this count is referred to as a betweenness value similar to Chen et al.’s (2016) work. For the nodes having the highest betweenness values, we calculate cosine similarity for each node using node’s and target node’s vectors, where values closer to 1 indicate high similarity.² The paths are constructed using the default cost function (1) and heuristic function (3). We can observe from Figures 5-7 that the nodes that appear the most as part of a treatment path have above 75% similarity with the target node/features for all diagnoses, which is an indication of potentially effective candidate treatment procedures.

Figure 5. DAG evaluation.

Figure 6. Tree evaluation.

Figure 7. Graph evaluation.

5. Parallel processing and scaling

To be able to perform graph construction and path generation efficiently, we also construct a serverless API-driven pipeline based on Django (Django Software Foundation, 2022), AWS Lambda Functions (Amazon Web Services, 2023b), and S3 (Amazon Web Services, 2023a). Figure 8 demonstrates the pipeline logic, where the client sends a request with patient features and the diagnosis. Serverless Lambda Functions start and perform the graph construction and optimal treatment recommendation, then return the generated path and image URL; meanwhile, the images get uploaded to the S3 storage bucket. By calling the URL of the image, the backend makes a request to S3 storage and returns the constructed graph. The framework for calling the Lambda Function can itself be deployed on Lambda, making the overall pipeline available on demand, providing efficiency. This approach provides high scalability, meaning for newly admitted patients, the results can be obtained efficiently as the pipeline is serverless and fast, as the process is parallel, with each process running in its container environment. Further instructions on how to test the pipeline on the local Django server are provided in the project repository.

Figure 8. API-based serverless path recommendation pipeline.

5.1. Setting

The experiments are conducted in Python, using the Networkx library (Hagberg et al., 2008). The software pipeline is constructed using the Django backend framework, and the serverless cloud functionality is achieved with AWS Lambda and S3 services.

6. Discussion

In this work, graph construction algorithms, cost, and heuristic functions, along with a software pipeline were proposed with the aim of finding the correct medical treatment for patients across multiple recovery stages using graphs where each node denoted patient features. As a limitation, we acknowledge that the approaches need further medical testing and verification. In addition, due to the stage-based dataset pre-processing and novel graph construction logic for performing a graph-search-based recommendation, which also provides dataset size limitations, the effectiveness of the approach with existing works can be validated only in terms of methodological and experimental similarities. K-Means and Cosine Similarity-based validation are provided for the evaluation of the methods. As part of future work, providing a graph construction pipeline for an increased number of time-stages, which can also include dataset benchmarking, scalability in searching performance for bigger graphs or trees, and further optimization of the algorithms, can be considered.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was partially supported by the RA MESCS Science Committee and BRFFR in the frames of the joint research project SC No21SC-BRFFR-1B029 and by the research project [21T-1B314].

Notes on contributors

Karen Gishyan

Karen Gishyan. Bachelor of Arts in Business, American University of Armenia, Yerevan, Armenia 2015–2019. Master of Science in Computer Science, University of Bath, Bath, UK 2019–2020. PhD student at the Institute for Informatics and Automation Problems of the National Academy of Sciences of the Republic of Armenia. Research interests: Intelligent Control and Cognitive Systems, Software Systems and Design, Artificial Intelligence and Time Series Analysis.

Hasmik Sahakyan

Hasmik Artem Sahakyan. Graduated from Yerevan State University. She received her PhD in 2002 and her doctorate in 2018. She is currently the scientific secretary of the Institute for Informatics and Automation Problems of the National Academy of Sciences of the Republic of Armenia and leading researcher at the Discrete Mathematics Department. Scientific interests: Combinatorics, Discrete Optimization, Discrete Tomography, Artificial Intelligence, Data Science, and Data Mining.

Levon Aslanyan

Levon Hakob Aslanyan. Graduated from Novosibirsk State University in 1968. He received his PhD in 1976, and his doctorate in 1997. Professor since 1997. Corresponding member of National Academy of Sciences of the Republic of Armenia, 2014. Currently heads the Department of Discrete Modeling, Analysis, and Recognition Technologies at the Institute for Informatics and Automation Problems of the NAS RA. Research interests: mathematical logic, discrete mathematics, mathematical theory of pattern recognition and artificial intelligence.

Notes

1. Code source: https://github.com/karen-gishyan/project

2. For Algorithm 1, similarity is calculated with the target node’s feature vector, for Algorithm 2, it is calculated with the averaged vector of 3rd stage’s target feature vectors.

Appendix A.

Algorithms

Algorithm 1 DAG construction pseudocode

1. Obtain “train” patient features for the first and second stages.

2. Obtain those “train” patient features for the third stage, which have been discharged to Home (target node will be selected from one of these features).

3. Obtain first stage “test” features as test nodes.

4. Initialize a DAG with the start node for each “test” node.

5. For each “testing” node, iterate over stage 1 “training” features, calculate closeness based on cost function (1), select top five closest “training” features from stage 1, make DAG nodes, construct edges, and assign closeness values to the edge (If there are no actual data for the given “training” feature for that stage, assign a fixed penalty of 2 * *10).

6. For each closest node from step 5, iterate over stage 2 “training” features, calculate closeness based on cost function (1), select top five closest “training” features from stage 2, make DAG nodes, construct edges, and assign closeness values to the edge (If there are no actual data for the given “training” feature for that stage, assign a fixed penalty of 2 * *10).

7. For each closest node from step 6, iterate over stage 3 “training” features, calculate closeness based on cost function (1), select top five closest “training” features from stage 3, make DAG nodes, construct edges, and assign closeness values to the edge (If there are no actual data for the given “training” feature for that stage, assign a fixed penalty of 2 * *10).

8. Search the graph from step 3 with A Star algorithm, using the start node and the target node.

Algorithm 2 Graph/Tree construction pseudocode (recursive)

1. Obtain “train” patient features for the first and second stages.

2. Obtain those “train” patient features for the third stage, which have been discharged to Home (target node will be selected from one of these features).

3. Obtain first stage “test” features as test nodes.

4. Select target features as the average (over the rows) for the feature matrix from step 2, and choose a cost function and a respective threshold.

5. Initialize a Graph with the start node for each “test” node. Initialize a “frontier” queue, dictionary storing visited nodes, dictionary storing cycle-resulting nodes (for the graph version).

6. Base case 1: Terminate the algorithm if there are no more nodes in the “frontier” queue to explore.

7. Base case 2: Terminate the algorithm if the node is at a depth of 3 or deeper, has features that are better than the “start node” features, mark the solution as “acceptable”.

8. Base case 3: Terminate the algorithm if the node is at a depth of 5, marking the solution as “not found”.

9. Base case 4: Terminate the algorithm if the node’s features are within the threshold, mark the solution as ‘goal’.

10. Inductive case. Select a node from the queue, calculate similarity between target features using cost functions (1), (7), or (9). If one of the base cases is not met, iteratively find the three most similar features from the features from steps 1 and 2, consider them as child nodes.

11. If tree version: For a given child node, if it results in a cycle, remove the edge connecting the parent to the child node. Add the node to the frontier queue with default features if the algorithm is deterministic; otherwise deviate the features with some probability before adding them to the queue. Relabel the node if it has non-deterministic features. After considering all child nodes, do a recursive call.

12. If graph version: Add a given child node to the frontier queue with default features if the algorithm is deterministic and has not already been visited; otherwise deviate the features with some probability before adding to the queue. Relabel the node if it has non-deterministic features. After considering all child nodes, do a recursive call.

13. Search the Graph from step 5 with Dijkstra’s algorithm, using the start and target nodes (the tree version will always have a single solution path if one exists).