 Research
 Open Access
 Published:
Concentration optimization of combinatorial drugs using Markov chainbased models
BMC Bioinformatics volume 22, Article number: 451 (2021)
Abstract
Background
Combinatorial drug therapy for complex diseases, such as HSV infection and cancers, has a more significant efficacy than singledrug treatment. However, one key challenge is how to effectively and efficiently determine the optimal concentrations of combinatorial drugs because the number of drug combinations increases exponentially with the types of drugs.
Results
In this study, a searching method based on Markov chain is presented to optimize the combinatorial drug concentrations. In this method, the searching process of the optimal drug concentrations is converted into a Markov chain process with state variables representing all possible combinations of discretized drug concentrations. The transition probability matrix is updated by comparing the drug responses of the adjacent states in the network of the Markov chain and the drug concentration optimization is turned to seek the state with maximum value in the stationary distribution vector. Its performance is compared with five stochastic optimization algorithms as benchmark methods by simulation and biological experiments. Both simulation results and experimental data demonstrate that the Markov chainbased approach is more reliable and efficient in seeking global optimum than the benchmark algorithms. Furthermore, the Markov chainbased approach allows parallel implementation of all drug testing experiments, and largely reduces the times in the biological experiments.
Conclusion
This article provides a versatile method for combinatorial drug screening, which is of great significance for clinical drug combination therapy.
Background
In the practice of clinical treatment, a single drug often fails to achieve the desired efficacy because the single drug in general aims at a single target of diseased cells and cannot remedy all aberrantly functioning pathways because of the robustness of organisms. The drug may also have poor safety profiles owing to various factors [1], including compensatory changes in cellular networks upon drug stimulation [2], redundancy [3], crosstalk [4], and offtarget activities [5]. In contrast, drug mixtures are generally more effective than single effectors because multiple drugs simultaneously act on different pathways and cell targets, potentially leading to higher efficacy and lower toxicity because of drug synergy [6]. Therefore, in the clinical treatment of complex diseases, such as parasitic nematode infections or herpes simplex virus (HSV), a variety of drugs have been used in combination for treatment improvement [7]. The infection of parasitic nematodes (or roundworms) poses a serious safety hazard to humans and livestock [8], and the anthelmintics (or antinematode drugs) are highly susceptible to drug resistance. It has been proved that a variety of combinations of multiple anthelmintic drugs, rather than a single medicine, can enhance the deworming effect [9]. In the case of the eradication of wildtype Caenorhabditis elegans worms, it is more effective to use four combinatorial drugs (levamisole, pyrantel, tribendimidine, and methyridine) than single drugs [10]. Traditional treatments of HSVI, one of the most common sexually transmitted infections, often include virusspecific drugs, which are effective at the beginning but exhibit limited longterm efficacy as drugresistant strains develop. However, a combination of six drugs (IFNα, acyclovir, IFNγ, ribavirin, IFNβ and TNFα) was demonstrated to be the most promising therapy for the reason that the drugs in the combinatorial treatment can act simultaneously on the multiple pathways and cellular protein complexes, and, therefore, regulate all relevant pathways, potentially blocking HSVI replication [11]. Combined use of multiple drugs is also a common practice in the treatment of cancers to achieve higher efficacy and potency. For example, in the treatment of nonHodgkin’s lymphoma, the drugs, pirarubicin, velet, cytarabine and prednisone, are usually used in combination, which the chemotherapy effect is remarkablely enhanced [12].
However, owing to the inherent complexity of biological systems and internal structure of cells and, particularly, to the huge searching space, it is extremely challenging to effectively and efficiently to determine the optimal drug mixture from all possible drug combinations through trial and error. For example, there are n drugs and each drug has m concentration candidates, it is necessary to find the optimal drug mixture in the space of \(m^{n}\) combinations. Obviously, as the types of drug increase, the number of combinations increases exponentially, and it is impossible to test all cases of drug combinations because it takes a considerable amount of time to perform the testing experiments. Therefore, it is important to explore how to reduce the number of experiments and predict the optimal combinatorial drug concentration accurately and quickly.
For these reasons, the optimization of drug combination has attracted considerable attention in recent years, and several methods for predicting the optimal combinatorial drug concentrations have been proposed [13,14,15,16,17,18,19,20]. A feedback system control (FSC) method was developed to search for optimal synergistic combinatorial drugs for the treatment of diseases [16]. The FSC method starts with a set of initial concentrations of combinatorial drugs with defined drug doses, and the efficacies of the combinatorial drugs on the cells at the given concentrations are evaluated according to the phenotypic output response of the cells. Then, the next predictions of the concentrations of drug mixtures are conducted based on the previous drug testing results with a certain searching algorithm, such as the Gur game (GG) algorithm, modified Gur game (MGG) algorithm, differential evolution (DE) algorithm, streamlinedfeedback system control (sFSC) algorithm, continuous adaptive population reduction (CAPR) method, and the FSC method iteratively approaches a globally optimal combinatorial drug mixture [19]. However, in some cases, these algorithms may degrade the overall performance of FSC owing to the inherent shortcomings of these algorithmic frameworks [17,18,19,20]. FSC with the GG and MGG algorithms often falls in oscillatory curves instead of giving a convergent output. It converges too early to a local extremum with the DE algorithm, thereby forming a premature convergence phenomenon, and it lacks a unified parameter controlling strategy with the CAPR algorithm to satisfy various applications. As the search process of sFSC method is based on an ‘iterative cycle’, searching for the optimal concentration requires more iterations. Furthermore, the FSC iterates its searching process, in which the next iteration requires biological experiments with the predicted combinatorial drug doses for further evaluation and prediction. Thus, the optimization of combinatorial drugs with FSC is quite inefficient because a significant amount of time is spent on the testing experiments.
Markov chain is one of the most important and fundamental algorithms in the field of machine learning, and is being used widely in the analysis of biological data and in bioinformatics area. It is novel to use the Markov chain method to model the combinational drug optimization and search for the optimal concentration combination. In this paper, an optimization method based on Markov chain models is proposed to search for optimal combinatorial drug concentrations with excellent performance. In this method, the searching process of the optimal drug concentration is converted into a Markov chain with \(N = m^{n}\) state variables representing all possible drug combinations, where n refers to the number of drugs, and m is the number of discretized concentrations for each drug. This Markov chain can be depicted by a network of N nodes in the space of \(R^{n}\), where the nodes refer to the state variables. Assuming that all the possible drug combinations have equal probability to be the optimal mixture without having prior knowledge about the efficacy of the drug mixtures, a matrix of transition probability can be initialized so that the stationary distribution vector of the Markov chain has an equal value of 1/N for all its states. Then the searching process for the optimal combinatorial drug concentrations is equivalent to updating the transition probability matrix and seeking the the state with the maximum value in the stationary distribution vector.
The proposed method was validated by both simulation and biological experiments. In the simulation experiments, the proposed Markov chainbased method was compared with the five benchmark algorithms (GG, MGG, DE, CAPR and sFSC) in the FSC framework. In biological experiments, the survival rate of cells under two combinatorial drugs is regarded as the response function, and the Markov chainbased method was compared with GG and MGG in FSC. The results of the simulation and biological experiments prove that the algorithm based on the Markov chain outperforms the selected benchmark algorithms in terms of accuracy and efficiency. In summary, this study provides a versatile, novel method for efficiently optimizing combinatorial drug concentrations, and the work is of great significance for clinical drug combination therapy.
The remainder of this article is organized as follows. First, the preliminary theories of the Markov chain are discussed briefly, and the Markov chainbased method is presented. Then the simulation and biological experiments are described, and the experimental results are discussed to compare the performance of the proposed Markov chainbased method with other benchmark algorithms. In the last, we conclude this article.
Methods
In this section, first, some basic theories of the discretetime Markov chain are briefly reviewed. The optimization problem of combinatorial drug therapy is formulated with assumptions, and the general idea of the Markov chainbased approach to the optimization of drug combinations are described. The detailed algorithms in the cases of one drug and two drugs are given in Additional file 2: Fig. S1, Additional file 3: Fig. S2, Additional file 4: Fig. S3, Additional file 5: Fig. S4, Additional file 6: Fig. S5, Additional file 7: Fig. S6.
Markov chain theory
A Markov chain is a special kind of Markov stochastic process with a set of discrete states. It starts in one of these states and moves successively from one state to another, satisfying the Markov property. Markov chains are a mathematical model to describe a process in which the next state of the system depends only on the present state, and not on the preceding states. In other words, the process loses its memory of the past over time.
Definition of a Markov chain: When \(\left\{ {X_{n} ,n \ge 0} \right\}\) is a random sequence taking values in a finite or countable discrete set, where \(\Phi = \left\{ {1,2, \ldots , N} \right\}\) or \(\Phi = N\) typically, the process \(X\left( n \right) = X_{n}\) for \(n = 1,2, \ldots\) is a Markov chain if
where \(n \ge 0\) and \(s_{n + 1} ,s_{n} ,s_{n  1} , \ldots ,s_{0} \in \Phi\), and the values taken by the random variables \(X_{n}\) are called the states of the chain. Moreover, if the transition probability \(P\left( {X_{n + 1} = s_{n + 1} \left {X_{n} = s_{n} } \right.} \right) = p_{ij}\) for \(s_{n + 1} = j\) and \(s_{n} = i\) is independent of \(n\), then \(P = \left\{ {p_{ij} } \right\}\) is called the transition probability matrix.
Stationary distribution: For a Markov chain with the transition probability matrix \(P = \left\{ {p_{ij} } \right\}\), a probability distribution vector \(\pi\) is called a “stationary distribution” if \(\pi\) has entries \(\left\{ {\pi_{j} \ge 0,j \in \Phi } \right\}\) such that the following conditions hold.
where \(\pi = \pi P\) is called the “balance equation.”
Assumptions
Before introducing the combinatorial drug optimization method based on the Markov chain, it is assumed that there are two assumptions:

With the slow change in the concentration of the combination drugs, the effect of the drug on the experimental subject also changes smoothly.

The number of drug combinations is limited.
The above two assumptions are reasonable for the optimization of combined drug concentrations. The organism does not change dramatically under smooth input from the outside world. In the experiment, the concentration of the drug combination is a few discrete points. Under the above two assumptions, the combinatorial drug optimization problem can be expressed using a finitestate Markov chain.
Method description
First, a general example is depicted to illustrate the main idea of the method. Suppose that there are n kinds of drugs and each drug has m possible concentrations, then the state space \({\Phi } = \left\{ {1,2, \ldots ,m^{n} } \right\}\) represents the set of m^{n} combinatorial drug concentrations in an ascending order. Our goal is to find the optimal concentration from the state space. Here, the drug response function or death rate of cells can be represented by a normalized function \(f\left( x \right) \in \left[ {0,1} \right]\) for \(x \in \Phi\). The higher value of the \(f\left( x \right)\), the better effect of the drug combination at the corresponding concentration, leading to higher cell death rate. \(f\left( x \right)\) = 0 means that the drug combination at the concentration level \(x\) is completely ineffective while \(f\left( x \right)\) = 1 indicates that the drug combination achieves its best treatment efficacy. Our aim is to find the best concentration \(x^{*}\) for drug combination with the maximum value of the objective function \(f\left( x \right)\) as follows:
As shown in Fig. 1, in the case of three drugs, it is necessary to construct a threedimensional network structure of Markov chain. Each drug has m concentration levels and a total of \(m^{3}\) concentration combinations constitutes the state space \(\Phi = \left\{ {1,2, \ldots ,m^{3} } \right\}\). The states in \(\Phi\) represent the drug combination of different concentrations. For any \(i,j \in \Phi\); if \(f\left( i \right)\) > \(f\left( j \right)\), it is implied that the efficacy of the drug combination at the concentration level \(i\) greater than that at concentration level \(j\). Likely, in the general case of n drugs, an ndimensional network of Markov chain can be constructed, and the state space \({\Phi }\) consists of \(m^{n}\) states if each drug has m concentration levels.
In order to search for the optimal drug combination, a key assumption with the Markov chain model is that, for any state \(x\left( t \right)\), the state \(x\left( {t + 1} \right)\) at the next moment always comes from the current state \(x\left( t \right)\), and the states have a larger probability shifting to the direction with a larger objective function. In other words, at the step \(t\), if the objective function for the state \(x\left( t \right)\) is \(f\left( {x\left( t \right)} \right)\), then the state \(x\left( t \right)\) can select to transfer to its adjacent states to obtain the next state by comparing its function value with those at the adjacent states and choosing the state with a relatively higher objective function value for the next step. The benefits of this approach are obvious. As \(t\) approaches infinity, the state transfers to the optimal state x*, which means that the probability at the global maximum of the objective function is the greatest.
Reconsidering the Markov chain model described above, from the state transition diagram shown in Fig. 1, it is obvious that it is a random walk with any two adjacent states. Searching for the optimal drug concentration is equivalent to seeking the state with the largest steadystate probability in the stationary distribution. Therefore, a transition probability matrix P of \(m^{n} \times m^{n}\) is initialized and then updated iteratively for searching for the optimal drug combination with the Markov chain model. Firstly, according to the initial state, the corresponding transition probability matrix P is constructed, and two suitable experimental points are selected from the state space. Secondly, the transition probability matrix is updated and then the balance equation is solved to achieve the stationary distribution. After multiple iterations, the algorithm converges with a predefined criteria, and the maximum value in its stationary distribution is the corresponding optimal state sought, that is, the optimal combination of drug concentration levels.
It is noteworthy that the initialization of the transition probability matrix is not unique. Without having prior knowledge about the efficacy of the drug mixtures, it is reasonable to assume that all the possible drug combinations have equal probability to be the optimal mixture and a matrix of transition probability can be initialized so that the stationary distribution vector of the Markov chain has an equal value of \(1/N\) for all its states. In this study, the transition matrix is initialized such that, on the network, every pair of adjacent states has an equal transition probability to move back and forth between each other, and every state has the same transition probability to move to all its adjacent states. In particular, the state on the edge of the network has a certain probability to go back to itself. Then, the Markov chainbased approach to optimizing the combinatorial drugs turns into a process of repeatedly updating the transition matrix by comparing the efficacies of pairs of adjacent drug combinations and then computing the corresponding stationary distribution vector until a certain convergent criterion is satisfied. The steady state that has the maximal distribution probability is referred to as the optimal drug combination.
The general procedure of the optimization algorithm for combinatorial drugs based on the Markov chain model is described as follows:

Step 1: The Markov chain and the corresponding transition probability matrix are initialized according to the numbers of drugs and concentration levels.

Step 2: Suitable adjacent combinations of experimental points are selected.

Step 3: The transition probability matrix of the Markov chain is updated according to the difference in the drug response functions at the corresponding suitable experimental points.

Step 4: The corresponding stationary distribution is solved according to the updated transition probability matrix using the balance equation.

Step 5: It is determined whether the stationary distribution converges. If it converges, the algorithm stops; otherwise, it returns to the second step, or, when the predetermined number of iterations is reached, the algorithm stops.
A single drug and two kinds of drugs are taken as examples to introduce the searching algorithm we proposed in Additional file 2: Figs. S1–S6.
Simulation experiments and discussion
Simulation experiments were conducted to compare the performances of the Markov chainbased algorithm and five benchmark algorithms: GG algorithm, MGG algorithm, DE algorithm CAPR method and sFSC method. The principle of these algorithms and the relevant control parameters selection for these algorithms are introduced briefly in Algorithm 1. The simulation experiments were implemented with the response functions in the cases of single drug, two drugs and three drugs respectively.
Predicting the optimal concentration of singledrug
Three drug response functions of single drug are used to compare the performances of the Markov chainbased algorithm and the benchmark algorithms, including GG algorithm, MGG algorithm, DE algorithm, CAPR algorithm and sFSC algorithm. And the three response functions of single drug, corresponding to the curves in Fig. 2a–c, are defined as below respectively:
In the simulations, it is expected to find the global maximum points for the drug functions. As shown in Fig. 2, the GG and MGG algorithms can not successfully find the global optimum points for all three drug response functions, and sometimes oscillate around some states (Fig. 2d) or stays in a suboptimal states (Fig. 2e, f). The DE algorithm and the CAPR algorithm not only fail to find the global maximum for some functions, but also take much more time to achieve the solutions (Fig. 2g–i). The sFSC algorithm can find the optimal point but sometimes it need more steps (Fig. 2j–l). In the contrast, the proposed Markov chainbased approach always succeeds in obtaining the global optimal state with less iterations (Fig. 2m–o). A further comparison of the Markov chainbased algorithm and the five benchmark algorithms was made using the success rate and the number of iterations, as shown in Table 1, which demonstrate that the Markov chainbased algorithms is more reliable and efficient than the benchmark methods.
In these simulations, the transition probability matrices of the Markov chain models were updated and the corresponding stationary distributions of the states change accordingly. As shown in Fig. 2p–r, the stationary distributions \(\pi = \left( {\pi_{1} ,\pi_{2} , \ldots ,\pi_{N} } \right)\) change gradually with the update of transition probability matrices and are finally convergent. Finally, the shape of the stationary distribution resembles the shape of the drug response function, which explains why the algorithm we proposed is effective for searching for the optimal combinatorial drugs.
Predicting the optimal combinatorial concentrations of multiple drugs
Simulation experiments were also conducted to validate the Markov chainbased algorithm and to compare its performance with the benchmark algorithms in the case of two drugs and three drugs respectively. The two response functions of two drugs are a Rastriginbased function and a De Jongbased function, which are defined as below respectively:
In order to implement the simulation with \(f_{4}\) and \(f_{5}\), the entire range \([0,5]\), for \(x\) and \(y\) were evenly discretized into 20 distinct values.
The response function of three drugs, a ternary function, is defined as below:
This function has multiple stagnation points and a single maximum point. The process to analytical obtain the extreme value and the maximum value points is placed in the additional file 1. 100 pairs of \((x,y,z)\) including the extreme value and maximum points were randomly selected to represent different combinatorial drug concentrations and the drug response function values at these corresponding points were calculated. Then five searching algorithms are used to search for the maximum point.
The searching processes of the optimal drug concentrations in the case of multiple drugs with the Markov chainbased algorithm and the five benchmark algorithms are displayed in Fig. 3. Figure 3A(a)–(f) in the left column, Fig. 3B(a)–(f) in the middle column, and Fig. 3C (a) ~ (f) in the right column of Fig. 3 represent the search for the optimal solutions of the response functions \(f_{4}\), \(f_{5}\) and \(f_{6}\), respectively with the various algorithms, indicating the relationship between the number of iteration and the optimal solutions at each iteration. Each row of Fig. 3 represent the search for the optimal solutions of the three response functions of multiple drugs with a certain algorithm, and from the top to the bottom are the original GG algorithm, the MGG algorithm, the DE algorithm, the CAPR algorithm, the sFSC algorithm, and the Markov chainbased algorithm respectively. From Fig. 3, the same conclusion can be drawn as the case of single drug that the five benchmark algorithms can easily fall into local optimal values or may take more iterations to find the optimal value, while the proposed Markov chainbased algorithm can find the global optimal value more reliably, but with less iterations.
As stated in the section of Method Description, the essence of the optimization process of the Markov chainbased approach is to update the matrix of transition probability so that the corresponding stationary distribution of the states converges to a pattern, which has a similar shape to the drug responses. As shown in Fig. 4A(b)–(f) and B(b)–(f), the patterns of stationary distribution of the Markov chain models for the response functions of two drugs are displayed at different iterations, and the number of interval iterations between each graph is 10 steps. Figure 4A(a) and B(a) are the patterns of the response functions, the Rastriginbased function and the De Jongbased function,respectively As the number of iterations increases, the smooth surfaces of the stationary distributions gradually converge in shape to the patterns of the corresponding response functions of two drugs.
Based on the results of the simulation experiments, the performances of the Markov chainbased algorithm and the other five algorithms were evaluated and compared with the measurements of success rate and number of iteration. The success rate of a algorithm indicates the effectiveness and reliability of the algorithm to find the optimal solutions, and the number of iterations represents how fast and efficient the algorithm is. The same evaluation measurements have been used to compare algorithm performance in literature [20], and they can effectively evaluate the reliability and efficiency of an algorithm.As shown in Table 1, the performance comparisons of the algorithms were made in the cases of the drug response functions of single drug (\(f_{1}\), \(f_{2}\) \(f_{3}\)), two drugs (\(f_{4}\), \(f_{5}\)), and three drugs (\(f_{6}\)), respectively. The table lists success rate and searching iterations (# of iters) of each algorithm for every response function. For each function, 1000 simulation experiments were conducted. In each experiment, it is recorded as “success” if the maximum point was successfully found. The number of iterations at which the maximum point was found is defined as the iteration number of the experiment. The total number of successes divided by 1000 is the success rate of the algorithm, and the average number of iterations to find the maximum point is the searching iterations. The algorithm is regarded as effective if the optimized output is larger than the threshold (\({\uplambda } = 0.95\)) or the output we predicted is among the top \(P = 5\%\) (even if the results we predicted is far from the real maximum value). From Table 1 we found that the success rate of the proposed Markov chainbased algorithm is much higher than the other benchmark algorithms, and the number of iterations is also less than the others. Furthermore, the success rate of the Markov chainbased algorithm is 1 in all cases of simulations, indicating that the Markov chainbased algorithm can always achieve the optimal solutions in the simulation experiments. Thus, it can be concluded that the reliability and efficiency of the Markov chainbased algorithm we proposed are better than the five benchmark algorithms.
Unlike of the benchmark algorithms of GG, MGG, DE, CAPR and sFSC depicted in Additional file 1, the Markov chainbased algorithm can surmount the deficiencies that the benchmark methods have, and can always successfully predict the optimal concentrations of combinatorial drugs with excellent performance in all simulation experiments (Figs. 2, 3 and Table 1). The the state space of Markov chain models are constructed by discetizing the drug concentrations evenly and the experimental points can be selected randomly. The state with the largest steadystate probability in the stationary distribution is the output of the Markov chainbased method and such output is usually unique. Moreover, in the simulations, with the benchmark algorithms of GG, MGG, DE CAPR, sFSC in the FSC framework, the prediction of combinatorial drug concentrations at a certain iteration is conducted based on the drug response information at the previous iteration, and then the experiments with predicted drug concentrations are carried out for the next iteration of prediction. On the contrary, the Markov chainbased approach allows the experiments at the selected states or drug concentrations to be implemented simultaneously. This makes a significant difference for biologically optimizing the combinatorial drug concentrations since it usually takes a few hours or days to implement one iteration of biological experiments and the computational time is usually much less than the time spent in the biological experiments. Therefore, it make take more than a month to finish a concentration optimization of combinatorial drugs using the benchmark algorithms in the FSC framework. However, the Markov chainbased approach with parallel experiments at the selected set of drug concentrations could save lots of time in the biological experiments. From this perspective, the Markov chainbased algorithm is much more efficient than the benchmark algorithms.
Biological experiments and discussion
Cell culture
The cell lines used in this study were obtained from the School of Medical Device, Shenyang Pharmaceutical University (Shenyang, China). MCF7 cells (human breast cancer cell line) and BXPC3 cells (human pancreatic cancer cell line) were cultured in RPMI1640 (Thermo Scientific HyClone, Logan, UT, USA) containing 10% fetal bovine serum and 1% penicillin–streptomycin solution at 37 °C (5% CO_{2}).
Cell proliferation assay
Cells were plated onto 96well plates (8 × 10^{3} cells/well for MCF7 and 8 × 10^{3} cells/well for BXPC3) and allowed to attach for 24 h. Cells were incubated with free drugs dissolved in an appropriate cell culture medium at serial concentrations for 72 h. For treatments containing DOX and PTX, each contained nine concentrations ranging from 0 to 5000 nM according to DOXequivalent concentration with a total of 81 concentration combinations with six complex holes per concentration. Following incubation, 10 µL of cell counting kit8 (CCK8) (Dojindo) was added to each well in the dark and incubated at 37 °C (5% CO_{2}) for 2 h. After incubation, a microplate reader (Thermo, Multiskan FC) was used to measure the number of viable cells in each well of a 96well plate at a wavelength of 450 nm.
Performance comparison
The response functions of the cells, MCF7 and BXPC3, under the combinatorial action of two drugs, paclitaxel (PTX) and doxorubicin hydrochloride (DOX), were constructed based on the biological responses to compare the performance of the algorithm we proposed and the GG and MGG algorithms. Figure 5A(a) and B(a) are the two drug response functions. (The green circle drawn in the fig is the maximum point of the drug response function, the red square is the maximum point found using the GG algorithm, and the black square is the maximum point found using the MGG algorithm.) Fig. 5A(b) and B(b) are the performance of the original GG algorithm, Fig. 5A(c) and B(c) are the performance of the MGG algorithm, and Fig. 5A(d) and B(d) are the performance of the Markov chainbased algorithm to find the optimal combination of PTX and DOX. Figure 5A(b)–(c) and B(b)–(c) show the nonrobustness of the GG and MGG algorithms. From Fig. 5, we can draw conclusions similar to those in the simulation. The original GG algorithm can easily lead to falling into the local optimal value (Fig. 5A(b) and B(b)). As shown in Fig. 5A(c) and B(c), the MGG algorithm may take many iterations if the starting point and the optimal state are far away. It is obvious that the point found by MGG algorithm is the local optimal value as shown in the black square in Fig. 5A(a) and B(a).
In the Fig. 5A(d) and B(d), when the Markov chainbased algorithm is used, the global optimal combination can be found within only a few iterations. As the experiment and calculation are parallel, the proposed algorithm is much more efficient.
Shown in Table 2 are performance comparisons of the two GGbased algorithms and the Markov chainbased algorithm. Similar to the results of simulation, the efficiency and accuracy of the Markov chainbased algorithm are much better than those of the GG and MGG stochastic algorithms.
As shown in Fig. 6, the stationary distribution of the states (the combinatorial drug concentrations) of the Markov chain model are convergent to the drug response functions of the cells. Figure 6A(a) and B (a) are the two response functions of the cell lines, MCF7 and BXPC3, under the action of the combinatorial drugs respectively (the green circle is the maximum value of the drug response function). Figure 6A(b)–(f) and B(b)–(f) are the stationary distribution at the selected iterations. (The red circle is the maximum point found at the corresponding iteration.) The number of interval iterations between each graph is 10 steps. It can be concluded that the stationary distributions of the combinatorial drug concentrations varied as the transition probability matrices were updated. Finally, the pattern shapes of the stationary distribution is convergent and similar to the corresponding drug response functions. At approximately 20 iterations, the optimal drug concentrations can be found, which explains why the Markov chainbased algorithm performs very well for optimizing the combinatorial drug concentrations.
Conclusion
In this study, a novel Markov chainbased approach was proposed to solve the problem of the concentration optimization of combinatorial drugs. Its basic principle was introduced and the detailed algorithm was depicted in a general case. Both simulation and biological experiments were implemented to validate the proposed approach and to compare its performance with five benchmark algorithms, GG, MGG, DE, CAPR and sFSC, in a FSC framework. The simulation experiments were conducted with the response functions of single drug, two drugs, and three drugs, and the biological experiments were carried out in the case of two drugs with two types of cells, respectively The performances of the Markov chainbased approach and the benchmark algorithms were evaluated using two measurements, the success rate and the number of iteration, indicating the reliability and efficiency of a algorithm to seek the global optimum. The experimental results and the comparisons between the proposed method and the benchmark algorithms demonstrate that the Markov chainbased approach is much more reliable and efficient than the selected benchmark algorithms in the FSC framework. The Markov chainbased algorithms can always succeed in achieving the optimal solutions with much less computational iterations. Moreover, considering that the time spent in the biological experiments is much more than the computational time, the Markov chainbased approach allows parallel experiments at the selected set of drug concentrations and could save lots of time in the biological experiments, and therefore the proposed method will be much more efficient than the benchmark algorithms in the practical application.
Availability of data and materials
The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.
References
 1.
Nishimoto M, Koh H, Tokuwame A, Makuuchi Y, Kuno M, Takakuwa T, Okamura H, Koh S, Yoshimura T, Nanno S, et al. Drug interactions and safety profiles with concomitant use of caspofungin and calcineurin inhibitors in allogeneic haematopoietic cell transplantation. Br J Clin Pharmacol. 2017;83(9):2000–7.
 2.
Forsberg F, Brunet A, Ali TML, Collas P. Interplay of lamin A and lamin B LADs on the radial positioning of chromatin. Nucleus. 2019;10(1):7–20.
 3.
Fang X, Zhong G, Wang Y, Lin Z, Lin R, Yao T. Low GAS5 expression may predict poor survival and cisplatin resistance in cervical cancer. Cell Death Dis. 2020, 11(7).
 4.
Kitamura A, Takata R, Aizawa S, Watanabe H, Wada T. A murine model of atopic dermatitis can be generated by painting the dorsal skin with hapten twice 14 days apart. Sci Rep. 2018;8:9.
 5.
Kwon HJ, Kim W, Jung HY, Kang MS, Kim JW, Hahn KR, Yoo DY, Yoon YS, Hwang IK, Kim DW. Heat shock protein 70 increases cell proliferation, neuroblast differentiation, and the phosphorylation of CREB in the hippocampus. Lab Anim Res. 2019;35:21–21.
 6.
Wages NA, Chiuzan C, Panageas KS: Design considerations for earlyphase clinical trials of immuneoncology agents. J Immunother Cancer 2018, 6.
 7.
Debrah O, AgyemangYeboah F, Asmah RH, TimmyDonkoh E, Seini MM, Fondjo LA, Sight N, OwusuDabo E: SEROprevalence of herpes simplex virus type 1 and type 2 among women attending routine Cervicare clinics in Ghana. Bmc Infect Dis. 2018, 18.
 8.
Coghlan A, Tyagi R, Cotton JA, Holroyd N, Rosa BA, Tsai IJ, Laetsch DR, Beech RN, Day TA, HallsworthPepin K et al: Comparative genomics of the major parasitic worms. Nature Genet. 2019, 51(1):163.
 9.
Secula F, Erismann S, Cerniciuc C, Chater A, Shabab L, Glen F, Curteanu A, Serbulenco A, Silitrari N, Demiscan D, et al. Evidencebased policy making for health promotion to reduce the burden of noncommunicable diseases in Moldova. BMC Proc. 2020;14(Suppl 1):1–1.
 10.
Atakan HB, Xiang R, Cornaglia M, Mouchiroud L, Katsyuba E, Auwerx J, Gijs MAM: Automated Platform for LongTerm Culture and HighContent Phenotyping of Single C. elegans Worms. Sci Rep. 2019, 9:14.
 11.
Ding X, Sanchez DJ, Shahangian A, AlShyoukh I, Cheng G, Ho CM. Cascade search for HSV1 combinatorial drugs with high antiviral efficacy and low toxicity. Int J Nanomed. 2012;7:2281–92.
 12.
Naegelin Y, Naegelin P, von Felten S, Lorscheider J, Sonder J, Uitdehaag BMJ, Scotti B, Zecca C, Gobbi C, Kappos L, et al. Association of rituximab treatment with disability progression among patients with secondary progressive multiple sclerosis. JAMA Neurol. 2019;76(3):274–81.
 13.
Ding X, Liu W, Weiss A, Li Y, Wong I, Griffioen AW, van den Bergh H, Xu H, NowakSliwinska P, Ho CM: Discovery of a low order drugcell response surface for applications in personalized medicine. Phys Biol. 2014, 11(6).
 14.
Ding X, Njus Z, Kong T, Su W, Ho CM, Pandey S: Effective drug combination for Caenorhabditis elegans nematodes discovered by outputdriven feedback system control technique. Sci Adv. 2017, 3(10).
 15.
Wang B, Ding X, Wang FY. Determination of polynomial degree in the regression of drug combinations. IEEECaa J Autom Sin. 2017;4(1):41–7.
 16.
Wei F, Bai B, Ho CM. Rapidly optimizing an aptamer based BoNT sensor by feedback system control (FSC) scheme. Biosens Bioelectron. 2011;30(1):174–9.
 17.
Wong I, Liu W, Ho CM, Ding X. Continuous Adaptive Population Reduction (CAPR) for Differential Evolution Optimization. Slas Technology. 2017;22(3):289–305.
 18.
Wong PK, Yu F, Shahangian A, Cheng G, Sun R, Ho CM. Closedloop control of cellular functions using combinatory drugs guided by a stochastic search algorithm. Proc Natl Acad Sci USA. 2008;105(13):5105–10.
 19.
Yang J, Liu C, Wang B, Ding X. Feedback system control optimized electrospinning for fabrication of an excellent superhydrophobic surface. Nanomaterials 2017, 7(10).
 20.
Yoon BJ. Enhanced stochastic optimization algorithm for finding effective multitarget therapeutics. BMC Bioinform. 2011, 12.
Acknowledgements
None
Funding
This work is supported by the National Key R&D Program of China (Grant No. 2018YFB1304700), the National Natural Science Foundation of China (Grant Nos. U1908215, 61925307, 61903265, 91748212, 91848201, U1813210, 61821005, 61927805), the Instrument Developing Project of the Chinese Academy of Sciences (Grant No. YJKYYQ20180027), and the Key Research Program of Frontier Sciences, CAS (Grant No. QYZDBSSWJSC008).
Author information
Affiliations
Contributions
S.M. performed experiments, coding and manuscript writing. D.D. provided early code for the algorithm. W.W., Y.W., L.L. provided crucial guidance and ideas throughout the project. All authors read and approved the final manuscript.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 1.
Theory of Markov chainbased method and other benchmark algorithms.
Additional file 2
. Figure S1. Statetransition diagram of an Markov chain.
Additional file 3
. Figure S2. Drug response function between two adjacent concentrations.
Additional file 4
. Figure S3. Initializing the Markov chain and updating the corresponding transition probability according to two adjacent experimental points.
Additional file 5
. Figure S4. Twodrug case: a twodimensional network structure with N^{2} states.
Additional file 6
. Figure S5. Drug response function at concentration levels (x,y) and (x,y+1).
Additional file 7
. Figure S6. Initializing the Markov chain and updating the corresponding Markovchainbased transition probability based on the two neighboring states.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Ma, S., Dang, D., Wang, W. et al. Concentration optimization of combinatorial drugs using Markov chainbased models. BMC Bioinformatics 22, 451 (2021). https://0doiorg.brum.beds.ac.uk/10.1186/s12859021043645
Received:
Accepted:
Published:
DOI: https://0doiorg.brum.beds.ac.uk/10.1186/s12859021043645
Keywords
 Combinatorial drug optimization
 Markov chain
 Transition probability
 Stationary balance distribution
 Combinatorial therapy