Systems biology in disease modeling and control, analysis and perspectives

Omics and drug molecules become increasingly inﬂuential in identifying disease mechanisms and drug response. Because diseases and drug responses are co-expressed and regulated in the relevant omics interactions, the traditional way that grabbing molecular data from single isolated layers cannot always obtain valuable inference. Also, adverse effects exist in drugs that impair patients, and launching new medicines for diseases is costly. To resolve the above difﬁculties, systems biology is then applied to predict potential molecular interaction elements by integrating omics data from genomic, proteomic, transcriptional, and metabolic layers. Combined with known drug reactions, the resulting models improve medicines’ therapeutical performance by re-purposing the existing drugs and combining drug molecules without off-target effects. Based on the identiﬁed computational models, drug administration control laws are designed to balance toxicity and efﬁcacy. This review introduces biomedical applications and analyses of interactions among omics and drug molecules for modeling disease mechanism and drug response. The therapeutical performance can be improved by combining the predictive and computational models with drug administration designed by control laws. The challenges are discussed for its clinic uses. The high mortality of many diseases prohibits human longevity. Therapy design is needed to suppress disease progression and aid organisms to recover from abnormal states 1 . However, the cost of launching new drugs is increasing due to the long-term safety procedures in clinical trials 2 . Drug side effects and toxicity should also be considered in treatment design due to its damage to patients. Side effects occur when drugs bind with unexpected targets and disturb normal regulations. The side effects of new medications are unknown and may vary dramatically for different individuals 3 . For example, a new drug Torcetrapib has been designed for cardiovascular disease 4 , but the drug’s uncertainty might cause severe side-effects of hypertension 5 . At the same time, drug toxicity is considered in drug overdose, drug-drug interactions, addition to gene and protein network, the metabolic networks enable researchers to investigate protein functions in the network separately and identify their steady states 7 .

. Drug and target interactions. Drugs may bind with off-target proteins that induce side-effect. Drug molecule interactions may improve or inhibit the efficacy, which offers chances for drug repositioning and drug combination. quantified drugs clinical relevance 56 . Besides, when applying drug combinations, different drug agents may have interactions that induce unexpected side effect to patients or reduce the drugs efficacy 59,60 . It requires researchers to predict potential drug-drug interactions (DDI), which can be explored by expanding DDIs with the shared target components of the drugs. The enriched DDIs can then identify potential targets of the drugs being tested, which reveals new therapeutic use and combination with surprising efficacy or for the diseases that the drugs do not initially aim at 3 . An abridged general view of drug and target interactions is shown in Fig.3.
The text mining-based predictions of molecular interactions rely heavily on the existing reports and references 61 , which indicates that insufficient clues would result in inaccurate predictions for the combination of omics and drug molecules. The predicting performance also depends on the setting of the cutoffs. Instead of the hard cutoffs, soft cutoffs (which may be similar to the soft thresholds 52 ) may be used when considering the generality for different data sets. Some of the frequently used databases for omics interaction analyses are listed in Table.1.

Analysis of static modelling
Importance quantification Although the models of diseases or drugs have been updated and well-studied for decades, the actual biosystems are far more complicated from complete modelling 63 . Though, not all the multi-layer information will be used after quantifying the response measure for efficiency. The importance of the interactions requires measuring and ranking for reducing the complexity of the networks or generating the corresponding weights. The simplified network will include only those most relevant molecules from a biological standpoint and directly associated with questions we concern about 75 . Thus, the initialization of static network analyses can be the importance quantification and ranking of its components, such as genes, interaction importance 63 , treatment or drug efficacy 48 . Z-score quantification have been used to measure the importance of the shortest path lengths between targets of drugs and proteins correlated with the cardiovascular module, which reduces the bias brought by the highly connected nodes 55 . Network proximity was applied to measure the importance of network distance between a drug and a given disease. In PPI networks, the edge weights quantify the confidence of interactions 50 . The importance of cell physiology has been measured by the observable concentration levels, the change rate, and the initial concentration of ligand proteins in a cell signaling models 76 . For example, in a protein interaction network 63 , the top 20 ranked cardiac septal defect candidate genes have been selected to represent most of the functional interactions. The known disease genes were seeds in the network, and the prioritization defined the rank of the importance of all the other nodes towards the seed proteins. Rank of the same proteins would be unique since no overlapping happens within the disease network. All the other proteins contributed to the vectorized interaction measurement toward each of the specific seed proteins. In DTI network, the shared drug-substructures prioritize DTI pairs to build the links between the new chemical entities and known drugs 61 . The 'resource-spreading processes' that utilized adjacency of these interaction pairs were proposed to quantify the interaction probabilities and make predictions. Also, drugs can be repurposed with quantified drug-target binding affinity. For instance, Fabahistin is designed as a symptomatic antihistamine, while it can be used for Alzheimer disease since it has higher off-target affinity to the disease-related ligands 32 . Protein information based on interactions described in published reports. This interaction set is expected to be biased toward known disease genes. 14,31,44 Similarity analysis Similarities describe the closeness between the components in static network. It helps convey the knowledge from the well-studied components to the unexplored ones, and it has the potential for novel clinical findings when two components share the same neighbor 3 . From the view of DTI and PPI, drug target proteins in the significantly differentially expressed signature modules can share the same protein family for the given diseases 31 , which indicates that similar drugs may be suitable to treat a cluster of similar diseases 65 . Drug pairs with more similar structures have significantly higher drug sensitivity correlations by activity area, which suggests that drugs with similar chemical structures show similar inhibitory effects across the cell lines tested 54 .
There are different types of similarities in molecular interactions. Cell lines with similar gene expression profiles exhibit similar drug responses, which can be used to predict the drug response of cancer patients to therapeutic agents 54 . The sequencesimilarity ranking has shown good performance for prediction based on ligand similarity 32 . Similarity in phenotypic side-effect contributes to drug targets identification. This has been validated experimentally that the antiulcer drug rabeprazole binded the serotonin receptor HTR1D and inhibited the dopamine receptor DRD3 through similarity in drug-target connections 71 . The global structure similarity has been used to build connections between the off-targets 23 . In DDI, the similarities of positive DDI pairs were significantly higher than those of random DDI pairs 77 . With similarity inherence and the collective inference approach, the Probabilistic Soft Logic language 65 conveys the similarity of drugs (targets) to new targets (drugs) through the known DTIs.
The similarity in a network-based inference model of drug and target interactions can be derived by the drug-drug structure, the target-target genomic sequence and the drug-target topology network 78 . When identifying off-target effect in drug activity 56 , the two-dimensional structural similarity of a drug to the target's ligand set was quantified as an E value using the similarity ensemble approach. This method computes whether a molecule will bind to a target based on the common chemical features with known ligands. The pathological similarities used for DTI prediction can be quantified by the Partial Spearman correlation calculation 31 with module response scores. The binarized fingerprint of side effects and the Jaccard index can be used to compute similarities of drug-target pairs 79 . Rules of how the shared drugs and targets convey similarity have been compared 65 .
Data fusion Combination of high throughput data contributes to the prediction of biomedical outcome, which means data synthesis of specific module functions and network motifs makes data more informative. Useful information will be maintained after fusion. For omics information, the data fusion is often used by incorporating multiple data types and sources to expand the level of omic and pharmacological data from molecule to organism 3 . Similarity contributes to the fusion of the separated 7/20 networks which consist of data from large scales such as genomics and its transcriptional products 62 . Data fusion has been used to refine essential and common information in cancer type identification 62 . The subnetworks, which were constructed for data of each scale by similarity measurements with non-linear fusion, eventually converged to one single network by fusing iteratively. Higher similarities from various levels were fused together in the resulting network. Lower similarities were discarded since they were regarded as noise. The fusion of different data and knowledge sources has been used to study drug repurposing by applying the tri-factorization algorithm 40 . Data fusion by matrix factorization in drug synergy was applied to identify the proposed method's potential target protein combinations. Besides, kernels or kernel functions that extract data features by matrix multiplication can be powerful data fusion tools. For drug side-effect prediction, the weighted kernels of the drug spaces and the side-effect spaces have been fused for information integration 80 . A multi-kernel based method, which refers to a composite kernel with a weighted average of each single kernel, has been utilized for omics integration in predicting the genomic risk 13 . The robustness of the proposed method for different disease analyses was retained by choosing the appropriate kernels' predictive genomic regions. Additionally, whether the fused information contains enough features from the original data depends partly on the kernels' selection or design. Some recent progresses of omics integration by kernel learning have been summarized 81 . While the integrated information generalizes the scope of predicting model, the weights that are used to overlay multi-layers should be carefully assigned, which can improve the accuracy of the predictions.
Learning-based methods for clustering and classification The evolution in learning-based networks has proven to be a promising architecture for efficient learning from massive datasets 3 . Learning-based approaches efficiently utilize computer power for extracting more complex features from high-throughput data in hyper-dimensions. It can identify global pattern of network from multi-layer omics data by combining them into fewer layers, which retains the maximal covariation from each original level. The learning-based methods can be supervised (with labels for training data) or unsupervised (without labels). A supervised learning approach enables the model to predict unknown part (links) of the network based on the known interacting molecules 58,64,70,82 . Support Vector Machine (SVM) is a classical supervised learning model for classification problem that can be applied to refine topological information from network structure. The multi-class SVM has been used to predict the therapeutic class of FDA-approved compounds and the classifier is trained with integrated drug similarities 72 . Kernels in SVM measure the features between gene pairs to train the classifier 83 , and the classifier makes the binary prediction for the interactions between the existing molecules and the incoming components 64 . Moreover, the convolutional neural network, which outperforms SVM since it does not require manually defined features, has been applied to extract DDIs with text mining on biomedical information 60 . Besides, a supervised bipartite inference can predict unknown target proteins of drugs in drug-target interaction 84 .
Deep learning (DL) methods outperform the previous learning-based algorithms as for the property of structure flexibility 85 . It has recently gained considerable attention as it allows the model to 'learn' to extract relevant features of the task at hand 86 . DL is also capable to extract molecular patterns by mining latent information from the network structures 66,69,87,88 . Given network structure, DL helps predict annotations 38 and bingding affinity 86 of the edges. A convolutional neural network can be used to refine molecular features from arbitrary network frame of different sizes and shapes 88 . This approach transformed 2-dimensional neural graphs into state matrixes with informative nodes and edges. Besides, auto-encoding automatically and dynamically extracts features from molecular structure, so that a minimalistic description can be derived for differentiability among samples 87 . An autoencoder (AE) can transform molecules directly into a numerical representation 69 . AEs can be used to denoise a single-cell RNA sequencing model 89 and to compresses computational dimensionality 3 in detecting the cell type-specific clusters in an ensemble clustering study 90 . The stacked AE has been used to identify the potential DTIs with information from drug molecular structure and protein sequences that generates highly representative features from multiple layers 68 . The feature descriptors were constructed by combining the molecular substructure fingerprint information and fed into the rotation forest for accurate prediction. The experimental cross-validation showed that it achieved superior performance on gold standard data sets with high accuracy. Additionally, the accuracy of the prediction can be assessed by area-under-curve derived from receiver-operating-characteristics curves 3,91 . Consequently, the classifiers generated from the clustered model can be utilized to predict the unknown functions of the nodes in PPI network 73 , unknown interactions between drug and targets 84 , potential interactions between proteins and the regulation between disease and genes 74,91 . A scheme of a convolutional neural network is shown in Fig.4.
The omics interactions identified in static framework reveal the possible connections within organisms. By selecting the disease-related omics from the resulting predictions, one can search the dynamic pathways in the regulation and identify key regulators for disease progression. With potential drug targets, drugs may be repurposed for more relevant disease. Interactions between drug molecules may either reduce or enhance drug efficacy, which may help design novel schemes for drug combination in disease treatment, such as drug cocktail 92 . Additionally, the clinical use of potential drug combinations or drug repurposing requires high confidence in the molecular interactions for patients' safety, so the accuracy of current algorithms can be improved, since it may not have good performance when testing on different dataset. While the commonly acknowledged problem in analyzing interactions between omics and chemical molecules based on database is data scarcity 22 .

Figure 4.
A convolutional neural network. This neural network is comprised of an input layer, convolutional layers (that extract features from hyperplanes of input data by projection/convolution), pooling layers (that reduce the spatial size and mitigate the locational sensitivity), flatten layer (that flattens the features and feeds them into the artificial neural network), fully connected layer (that learns nonlinear function of the extracted features) and an output layer. The input is the raw features such as omics sequence patterns, gene regulation annotations and patterns, omics interaction network motifs, molecular structures and structural associations, drug chemical structures, drug side-effect reports, and so on. The output can be the predicted classification (e.g., omics binding profiles) and regression (e.g., quantified molecular binding affinities) when obtaining new samples.

9/20
More accurate descriptors will be discovered as more data is collected.

Dynamic modeling of regulation in organisms
The term dynamic disease can be explained as a disease whose pathogenesis is mainly caused by the appearance of new dynamics of the organism behavior, independently of the underlying pathogenesis 93 . Disease outcome can be predicted by changes in patients' interactome 19 . In dynamic regulation, the systemic view tells us that genes work as part of complex networks instead of acting in isolation to perform cellular processes 94 . Multiple differential co-expression networks outperform single co-expression networks as for inference accuracy of gene pathways 34 . For instance, drug efficacy and reservoir contribution rates cannot be identified from decay-phase data alone 95 . Besides, the dynamics of diseases and drug reactions are informative for treatment design. For example, the viral reservoir dynamics are significant to understand the natural properties of ongoing viral progression with treatment 95 , which help researchers find a functional cure for diseases. Once the dynamic models have been derived, the performance of a kinetic model can be assessed by its ability to generate testable predictions 76 . Notably, the goal of constructing dynamic networks is to identify the biological patterns by modeling the time-course behaviors of gene expressions 27 .
Interactions in gene regulatory network Dynamic analysis in organisms focuses on the omics variations, which are the results of gene expressions. Genes interact with each other through their RNA and protein expression products, thereby governing the rates at which genes in the network are transcribed into messenger RNA 93 . The causal information in gene expression data, e.g., key drivers of complex traits in phenotype-related gene interaction network 11 , can be identified by variations in DNA and gene expression-related traits 27,39 . Gene regulatory networks (GRN) are thus used to detect potential mechanisms based on dynamic behaviors of epigenomic activity between normal and disease state 94 . Notably, GRN is regarded as the reverse engineering for investigating gene regulatory relations by going backwards with the observed gene expression 33 . The schematic diagram of GRN is shown in Fig.5a.
GRNs have been used to analyze infectious diseases 14,96 to detect gene regulations related to the infectious and viral mechanism. A combination of gene regulatory components including myeloid and lymphoid has been studied to identify cell fate specification 10 . The genes related to the altered magnitude of the plasma hemagglutination-inhibition response, which are known to be regulated by XBP-1 (a transcription protein related to immune system), have been identified in the study of human innate and adaptive responses to vaccination against influenza 70 through the GRN. Interferon-related genes regulation interactions were used to get the molecular signatures for antibody response prediction based on data from the vaccinated patients and the mice trials. A Differential SparsE Regulatory Network model 94 has been proposed to identify topological changes in gene-regulator dependence networks for cancer study. The gene expression levels were modeled by a sparse linear model that enabled the researchers to capture conditional dependencies efficiently from genome-wide expression data. In inference regulatory interaction networks 29 , nodes were genes coded with expression levels, and edges were the genes interactions. To construct GRNs, one can use the time-course microarray profiles with network inference 27 , the combination of regulatory proteins and disease-related genes 33 .
Signaling pathway and transcription Signal transduction pathway contains regulators between molecules in an organism. It attributes to changes in both gene expression and gene connectivity 34 . At the protein level, signaling pathways are made up of protein interactions covering almost all the biological functions in living cells. It captures the inter-and intracellular regulatory mechanisms of gene transcription and protein synthesis 29,33 . For example, the analysis of dynamic signal transduction networks helps explore cancer cell phenotype, tumor cell heterogeneity and drug resistance, since cancer is considered a genetic disease 36 . When analyzing the interferon-stimulatory DNA-sensing pathway, the researchers have measured dsDNA response and built the connections between the predicted genomics with a known signaling network to identify the disease-related proteins and the functional participants 96 . Errors in signal transduction can lead to altered development and incorrect behavioral decisions in organisms, whose dysfunction may result in uncontrolled cell growth or tumorigenesis 33,36,97 . The dynamics of signaling pathways can also be a powerful measurement to model molecular regulatory behaviors of fungal pathogen infections in a fungal signaling network 7 . Similar to DEG, multiple differential modules 34 have been used to predict both unique and shared gene modules across multiple differential co-expression networks. This work explored the systematic profiling of the transcriptome in a heart failure study based on RNA-Seq data. The individual pathways were identified by the specificity inference of multiple differential modules, and the enrichment analyses were also performed for genes in differential modules with known cardiovascular phenotypes. Notably, the pathway-wide association can also be used to identify valuable information from background noise and the context-specific logic of GRN 41 . Transcription factors partly determine the cell fate specification. It regulates the development of innate and adaptive cells of the immune system 10 . The dynamic subnetworks of transcriptional and translational architectures have been explored to construct and model the entire trigger mechanism of the innate response regulated by intercellular and intracellular heterogeneity 14

Analysis of dynamic modelling
Mathematical modelling for dynamic regulation Dynamics of omics data reflects organism's response to the changes of interior milieu or environmental factor. Variations in metabolisms and phenotypes can be huge even if most of the corresponding genes are the same 98 . Modeling of the dynamics exerts mathematical tools that quantify the rate of state change in dynamic regulation 27 in different conditions and time sequences. With knowledge of molecular interactions, the Differential Equations (DEs) can be used to capture how the system reacts to the variations caused by disease or drugs. These equations integrate information from different levels into a unified form, which agrees to the integration of omics information in systems biology 35 . DE models the dynamic variation in transcription 29 , gene expression in regulation 27 , the metabolite concentrations 99 , the reaction rate 35 , and the factors in signal transduction pathway 98 . It also quantifies signal flow in pathways and explores the effect of oncogenic mutations on dynamics of ligands 76 . For example, as for disease modeling, DEs capture the dynamics of viral infection, drug response, and interactions between HIV and the specific immune response in an HIV model 100 . As for drug response, DE has been applied to study drug intervention in an irradiation-induced cellular senescence model 101 . A system of DEs, which consist of drug concentration and populations of tumor cell, effector-immune cell and circulating lymphocyte, were summarized to model the growth, death, and interactions of these populations with chemotherapy 102 . Components in these DEs can be the metabolite concentration (in nonlinear biochemical models) 99 , the signal transduction molecules (in dynamic cell compartment models) 35 , transcription factors and regulatory site 29 .
Additionally, the established dynamic model will have the specified inputs and outputs that enable the prediction of crucial components in biological regulation, and a good dynamic model should be robust to perturbations, which allows investigators to track abnormal variations 24 . With more interactive terms identified to be involved in disease progression and drug reaction, DEs can model organisms' actual regulatory mechanisms more precisely. Note that molecules should not be added into model equations unless it is of significance for the given processes.
Parameter estimation in dynamic models The main task of parameter estimation is to find the best-fit parameters to characterize physical processes and reproduce experiments 103,104 . In systems biology, it is usually part of an iterative process  to develop data-driven models for organisms that should have predictive value 105 . One way to perform parameter estimation is to optimize the fit of a collection of parameters iteratively to a given data set 104 with random disturbance, which explores the parameter space extensively and limits the number of non-convergent solutions 101 . The least-square estimation (LSE) and the maximum likelihood estimation (MLE) are two widely used techniques for this task. LSE estimates parameters for a nonlinear dynamic system by minimizing a quadratic criterion function that comprises of the state estimates and measurements, subject to the constraints of the system 99,106 . MLE is used to search for the value of the parameter vector that maximizes the likelihood function, which is most likely to describe the observed data based on probability distribution 107 . MLE has been used to directly estimate the efficacy of antiviral drugs from the viral load data 95 . As an optimal linear estimator, Kalman Filter (KF) recursively generates the maximum likelihood estimates for a linear dynamic system from a series of noisy measurements 108,109 . It handles the approximate modeling of high-dimensional noisy data and a low number of samples 110 . Based on KF, its expanded versions have been utilized to deal with nonlinear models in biology. For example, the Extended KF, which is the nonlinear version of KF, can be used to estimate states and parameters in a nonlinear system 111,112 . Note that the noise uncertainty of KF is subject to Gaussian distribution. When dealing with non-Gaussian processes, Particle Filters (PFs) can be a good choice. By randomly drawing samples from numerical simulation, the Monte Carlo method can be used to estimate parameter values 95 and quantify their uncertainties 113 . As a sequential Monte Carlo method, PFs obtain weighted samples from the non-Gaussian posterior probability of the state to deal with nonlinear systems 114 . It has been studied that a PF with an orthogonal basis (that was used to approximate the posterior by an orthogonal series expansion) outperformed Extended KF when estimating the Wiener model parameters in an anesthesia delivery model 112 . The process of parameter estimation is shown in Fig.6a.
Sensitivity analyses assess how sensitive the models' outputs are to the fitted models' parameters changes. The integral of absolute sensitivities can be used to analyze the effect of relative changes within a time slot, and the irrelevant sensitivities which have peaks close to zero should be discarded 35 . Dynamic sensitivity analysis has found that treatment effectiveness got gradually reduced and all interventions became largely ineffective later in a cellular model for drug response 101 . Sensitivity can then quantify model uncertainty through finite differencing or variational equations 105 . Additionally, the variability that tells model variance can be assessed by a principal component analysis 101 . Besides, identifiability checks the reliability of the estimates and assesses how well the model explains experimental data 109,115 . The fitted model by parameter estimation may not have the generality for other datasets (e.g., clinical data varies for different patients). This may be solved by identifying the parameters that affect the variances between individuals.
Tools have been integrated for modelling and assessment. Systems Biology Markup Language based Parameter Estimation

12/20
Tool 103 has been developed for parameter estimation under specific customized requirements (e.g., changing signal strength). The evolutionary optimization with stochastic ranking was implemented in this tool. Recently, the 'Data2Dynamics' modeling environment has been widely used for dynamic modeling since it implements different efficient approaches for parameter estimation as well as uncertainty evaluation 116 . The functional package 'Uncertainpy' quantifies the uncertainty and performs sensitivity analysis with visualized results 113 . Time cost is considered in the modeling, which depends on the speed of getting solutions of those DEs. The solutions should be the optima for the estimation problem. The global optima for parameter estimation are always desired, leading the multiple approaches to converge to the same solution ultimately 99 . It is often used with a computational time limit to prevent an endless search 105 . The global optimization problem is stated as the minimization of a weighted distance measure between experimental and predicted values of the state variables in a nonlinear biochemical dynamic model with DEs 99 . When estimating parameters, calling the single clusters recursively may cause the absence of global parameters 35 . Also, we may get stuck in local optimum 105 , which is misleading since it matches the data well in the current data set. However, it will give a bad fit to data from other sources or fails to converge to a satisfying result 117 .
Once the numerical values of model parameters have been estimated, the approximation of dynamic omics regulation and pharmacokinetics from a systemic view is obtained for patients. Next, drug scheduling can be designed for patients' recovery based on the models.
Drug administration with control theory In control system, feedback signal provides error information to adjust the input and track the reference such that the system can be brought to the desired states. The meanings for feedback loops from biological aspects may vary. The negative feedback refers to inhibition and the positive feedback refers to activation from the downstream to the upstream 36 , and the feedback control in signal transduction networks includes positive feedback (signal amplification) and negative feedback (signal attenuation). Notably, feedforward signaling exists in specific transcriptional procedures 14 . Diverse properties of regulation in organisms come from a coherent or incoherent feed-forward loop by combining the inhibition and activation 36,93 . The feedback signal adjusts the drug infusion, which may consider the changing drug concentration, toxicity level, and cell population. Drug concentration in patients is delicate, especially for the fast body function change (e.g., patients do not react to anesthesia delivery during surgeries 118 ). The dangerous therapeutic window of drug concentration in plasma defines whether the concentration level is tolerable and the drug is effective for patients 17,67 .
The components in a control system include the states (e.g., time-evolved behaviors of the system), the inputs (e.g., variables that are manipulated by the control system such as drug doses), targets (that are the reference trajectories such as drug concentration), and the outputs (that are the elements we want to observe). For example, in an HIV-1 epidemic model 119 , state variables can be the density of virus, recombinant virus and the double-infected cells. Its control objectives are to use drugs to prevent new cells from being infected and decrease virus generation. Notably, dynamic models of drug response in real-world will be nonlinear for most of the cases. However, when taking drug infusion rate and plasma concentration as input and output respectively, a pharmacokinetic model can be treated as linear 118 .
As a classical control strategy, the Proportional, Integral and Derivative (PID) controller 120 maintains drug concentration at reference level in organisms. For example, in a neuromuscular blockade model 118 , a PID controller determines the dose of nondepolarizing types of muscle relaxant. The PIDs acted as local controllers in the proposed supervised multimodel adaptive control scheme. It was used to cope with the uncertainty in the system's dynamics caused by interpatient variability and time variations. The concise scheme of PID 120 makes it flexible for functional expansion. Based on a circulation model, a IPD controller was implemented in a chemotherapy control scheme that achieved and maintained drug concentration at the desired level 18 . When implementing PID, the tuning process can improve controllers' performance and robustness 121 . Recently, an auto-tuning algorithm for PID based on system specifications (e.g., response to a given input) has been proposed that almost met equivalent specifications compared to the design based on knowledge of the full process 122 . To achieve an adequate paralysis level during surgery, a time-varying PID controller paired with a tuning procedure assists the closed-loop drug infusion system recover from sustained oscillations in an autonomous way 112 . The instability of equilibrium caused by dynamic alternations in the PID model was addressed by a stable limit cycle, and bifurcation analysis has been performed on the nonlinear dynamics of the designed control scheme to prevent the oscillation at the unstable equilibrium points caused by plant uncertainty.
The optimal control problem is to find the control input on a time interval to minimize control errors (e.g., drive the systems along trajectories) and control efforts (e.g., less energy consumption) 123 subject to system dynamics. In chemotherapy, its goal is to find the balance between cancer cells elimination and drug toxicity. The drug concentration can be maintained within an effective therapeutic window, and side-effect caused by toxicity can be minimized. The general scheme is to design a proper cost function based on time evolution, and then solve the optimization problems subject to certain constraints. The control signals are the injection rates of drugs 124 . Usually, the objective functions are relevant to drug cost and patients tolerance during the treatment. As for chemotherapy model, the control problem can be to minimize cancerous cell growth, the average tumor burden 125 and the summation of the kinetic energies of all the cancerous cells 124 by searching the best distribution of a cumulative drug administration. The objective function in chemotherapy can be composed of tumor cell decrement, the average level of toxicity (tolerable drug concentration) during the treatment 18 . In a two-drug treatment model, the objective is the summation of the tumor volume and the doses of anti-angiogenic therapy and cytostatic drug 126 . In HIV, the goal is to maximize the benefit based on levels of healthy CD4+ T cells and immune response cells by reducing the systemic cost of anti-HIV drugs 100 .
Similarly, to eliminate the spread of HIV-1 by increasing the number of the normal CD4+ T cells, reducing virions, and minimizing the cost of treatment, the objective function consists of beneficiary T Cells and the weighted systemic costs of the therapy 119 . For multi-drug treatment, the goal can be to identify the optimal ratio of co-administration between chemotherapy agents and anti-angiogenic agent 124 . Constraints can be the boundary for toxicity level and the certain therapeutic window for drug concentration 18 . The constraints can also be certain bounds for drug doses and the volume of blood vessels in the neighborhood of the tumor 126 . Algorithms have been proposed for the constrained optimization problems. For instance, in a multi-objective genetic algorithm, the Pareto optimal sets have been used to search for lower values for the objectives simultaneously 18 . Nonlinear optimization can be solved by Bock's direct multiple shooting method with a numerical solution on a fixed control discretization grid 126 . The steepest descent method can be used to search numerical solutions iteratively to minimize cost function 124 . The numerical solutions can be derived with Miser3/Matlab 102 . To realize precision medicine, the objective functions should be evaluated for the individuals, and the solutions should be checked by the rank of the degree of dominance 18 . Model Predictive Control (MPC) can handle the mismatch between nominal and actual processes 127 . Its essence is the ability to handle constrained optimization problems [128][129][130] . MPC solves a finite horizon open-loop optimal control problem to obtain control actions with predicted states from models during each sampling period 129 . Similar to optimal control, objective functions in MPC mainly consist of disease-related molecules and drug infusion 131 . Control signals are formulated to minimize the difference between model values and the predefined reference trajectories (that can be learned based on a comparison of measurement to the desired values 130 ). The local asymptotic stability of the control law is guaranteed when giving sufficiently long time horizon 128 . The weight of each term in the objectives reflects the focus of the optimization problem. The control law calculated by minimizing the cost can be piecewise constants or linear 128 in the finite time horizon. For multiple control inputs (e.g., chemotherapy and immunotherapy), multiple MPC can be applied to design the optimal combination ratio with lower drug toxicity 132 . The loop of model predictive control is shown in Fig.6b. Notably, discrete drug therapy should be considered to allow the normal cells to rebuild in clinical treatment 119 . Besides, controllers can also be "learned" from complex situations by applying learning-based strategies, such that optimal parameterized control signals can be derived 133 . If the equilibrium exists in the chemotherapy model, it should be the case that the tumor is eliminated or cancer cell concentration is at a safe level, so that no more treatment is needed for recovery 124 . With the estimated parameters in dynamic models, various control laws can be designed for drug administration to suppress disease progression. Optimization has been widely used for therapeutic design since it handles patients' physical constraints and the efficacy of the treatment at the same time. Similar to the optimization problem in parameter estimation, the time cost of getting an optimal solution depends on the solvers' convergence. The solvers with higher efficiency are always desired. Also, to design personalized treatments, elements relevant to the variability among patients may be added to the objective functions.
The static network identifies potential elements to update the dynamic models for accurate disease mechanisms and drug responses. With experimental and clinic data, numerical parameters are fitted in dynamic models. Next, the drug infusion rate for better therapeutic performance is designed by control laws.

Outlook
Researchers can identify potential molecules or interactions in the existing biological networks with interactive patterns and make predictions for the incoming molecules. Since these molecules do not work individually, the functional clusters may be identified to account for organisms' particular functions. With a larger map of the biological process in organisms, the potential molecules or connections in omics and drug interactions can be identified for disease mechanisms and drug intervention. These existing models could have novel components to describe the same biological mechanisms in dynamic modeling. This makes the modeling more comprehensive and close to the real phenomena. The resulting models may contain the elements (from omics levels) to differentiate drug responses between individuals. In that case, it may enable personalized medicine that generalizes the models to serve patients with different health conditions. Besides, with the identified disease regulators through the molecular interactions and dynamic pathways, we may only use the drug molecules that bind these particular targets, which indicates that the key regulators should be the only targets for novel therapy design when desiring fewer drugs. On the one hand, with the detected off-targets, the side-effect can be minimized by selecting drug molecules properly according to the potential DTIs. For individuals, the identified off-targets' binding properties can be used to assess the side-effect of the drug molecules for all the available drugs. Then drug molecules can be combined by avoiding these off-targets (only if no DDI will be detected for side-effects within these drugs). Besides, the severity and frequency of the side effects should also be taken into consideration 79 . In this way, the safe drug combinations would be feasible for reducing the side-effect events. On the other hand, side effects should be avoided for safety issues, while this does not mean that all drugs that cause off-target effects will be abandoned. The treatments such as chemotherapy harm normal cells, so we may select cancer drugs that cause less damage 14/20 according to the identified DTI and combine with the advanced control strategies to get higher treatment performance but lower impairment. Besides, data sets in hand are always low-throughput 134 with high dimensionality. So, a good capability of noise resistance is needed for the feature extraction algorithms. The changing nature of omics data 8 may involve the irrelevant molecules in the models, which reduces the robustness of the prediction. An advanced algorithm, e.g., the learning-based algorithm, is desired to solve these constraints. Moreover, multi-constraints in systems biology make optimization an NP-hard problem 117 , which exists in parameter estimation and the constrained problem. A solver that avoids getting stuck in local optima and handles the constraints simultaneously is necessitated.