# Optimal Control of Gene Mutation in DNA Replication

Authors: Juanyi Yu; Jr-Shin Li; Tzyh-Jong Tarn
Journal: Journal of Biomedicine and Biotechnology
 (2012)

---

## Abstract

We propose a molecular-level control system view of the gene mutations in
DNA replication from the finite field concept. By treating DNA sequences as state variables,
chemical mutagens and radiation as control inputs, one cell cycle as a step increment, and the
measurements of the resulting DNA sequence as outputs, we derive system equations for both
deterministic and stochastic discrete-time, finite-state systems of different scales. Defining the
cost function as a summation of the costs of applying mutagens and the off-trajectory penalty,
we solve the deterministic and stochastic optimal control problems by dynamic programming
algorithm. In addition, given that the system is completely controllable, we find that the
global optimum of both base-to-base and codon-to-codon deterministic mutations can always
be achieved within a finite number of steps.


---

## Body

## 1. Introduction

Systems biology is an emerging academic field aiming at system-level understanding of biological systems. The early development of systems biology started in the late 1940s [1]. Recent progress in molecular biology has enabled us to gain information on the interactions among the underlying molecules from comprehensive experimental data sets. In general, a system-level understanding of a biological system can be derived from insight into four key properties: (1) the system’s structure, (2) the system dynamics, (3) the control method, and (4) the design method [2]. Equivalently, identifying related components and their interactions, gathering qualitative and quantitative information about the system’s evolution under different circumstances, achieving the desired outputs by controlling the input with appropriate definitions of inputs and outputs of the system, and reconstructing analogous systems by eliminating the undesired properties are four essential steps in systems biology done by collaboration among engineers, biologists, and doctors. Figure 1 shows a typical method of system construction and verification commonly applied currently. Control engineers construct models, run simulations, and predict the system behaviors. Biologists design and carry out the experiments and measure the output data. Control engineers revise and verify the models by comparing the predictions and experimental results.Figure 1
Typical analysis of biological systemsSystems biology is a cross-cutting research area connecting control engineering, biology, and medical science, as shown in Figure2. It provides a systematic view of the biological system and related medical interventions. It aims at understanding the bare function and integration function of the cell to reconstruct the biological systems with desired features. Control and automation play critical roles in this novel field not only by providing new technology and equipment for biologists to design and perform meticulous experiments, to take high-throughput measurements, and to analyze experimental data efficiently, but also by offering doctors new medical applications and improving the precision of medical manipulations. The wide range of aspects which control and automation have been applied to include, but are not limited to, gene regulation [3, 4], drug delivery [2, 5], and neuron networks [6, 7]. The equipment provided by control engineers includes, but is not limited to, nanodevices, biochips, cuvettes for electroporation, and gene guns. Biologists perform various biological experiments, such as protein synthesis and virus DNA modifications, to gather measurements for model revisions and verifications, to conclude theoretical and practical results from evidence, and to help medical practice. Doctors use both theoretical and practical results from biologists to perform tissue engineering, such as organ transplants and artificial tissue construction.Figure 2
Systems biology is a cross-cutting research area connecting control engineering, biology, and medical science. Sources: protein synthesishttp://www.anticancer.de/, liposome [8], corneal transplant http://www.avclinic.com/, microarray hybridization [9], cuvettes for electroportation http://www.en.wikipedia.org/, biochip http://www.clemson.edu/, nano robot http://www.molecularlab.it/.According to their scales, biological systems can be divided into three levels: the molecular level (nm), cellular level (μm), and tissue level (cm), analogous to the part, individual, and group, respectively. Molecular-level research focuses on how, when, where, and to what extent [10] a gene is expressed. The essential goal is to sketch a complete blueprint of genes by identifying the control sequences of coding DNA segments and their interactions. Cellular level research, in general, treats one cell as a plant in classical control theory and investigates the reactions of the cell to the changing environment, for instance, concentration changes of related chemicals. State-of-the-art medical therapies are primarily based on experimental results at the cellular level. Tissue-level research mainly concerns tissue reconstruction, artificial tissue substitutes, or tissue function recovery. The cell differentiation process is an important topic at the tissue level. Typical biological systems are collaboratively controlled at all three levels. Most current research work focuses on either cellular level or tissue level systems. Not much work has been done at the molecular-level. In contrast, understanding biological systems at the molecular level is crucial, since species have the same basic inheritance, DNA macromolecules, and follow a common rule in gene expression, the central dogma in molecular biology. Molecular level understanding of biological systems provides instrumental information about radical causes of many diseases and the genetic evidence of evolution. It also helps biologists to gain a better understanding of molecular level interactions, draw a complete blueprint of gene networks, improve existing means, create novel means to cure genetic diseases, and to elaborate on the theory of evolution. Recent technology in gene sequencing makes it possible to conquer the difficulty in measurement at the molecular level and to identify the nucleotide sequences of a particular DNA segment. Targeted sequencing is the most promising step toward maximizing the efficiency of the next-generation sequencing technology using polymerase chain reaction. The availability of DNA microarray makes it possible to accomplish tens of thousands of genetic tests for picomoles (10−12) of a specific DNA sequence.Researchers have applied various methods to model, simulate, and control the gene regulation processes. Early attempts to model and simulate gene regulatory systems are summarized in [10], including direct graphs, Bayesian networks, Boolean networks, ordinary and partial different equations, qualitative differential equations, quantative differential equations, stochastic equations, and role-base formalisms [10]. Other approaches include Petri nets [11], transformational grammars [12, 13], and process algebra [14]. Three important modeling methods in recent work are gene regulatory units viewed under compound control [4, 15–18], logic network models [19, 20], and base-to-base molecular-level formulation [21, 22]. The first modeling method quantitatively describes chemical concentration variations corresponding to external environmental changes at the cellular level. The second qualitatively illustrates the interactions among operons in gene regulatory units. The last modeling converts DNA segments to discrete vectors. In our paper, we adapt the base-to-base molecular level formulation to express state variables.Current obstacles in systems biology are obvious. The structure and dynamics of biological systems are sometimes unclear. Most existing models are constructed by data-driven or hypothesis-driven methods, with only partial information available. Due to the complexity of the systems and incomplete information, the mathematical models are usually formulated by modifying empirical equations or proposing heuristic equations. The parameters of proposed models are obtained by estimation methods. Although those models can disclose significant details of the system’s structure and dynamics, the inconsistency between theoretical and experimental results creates difficulties for control engineers to verify the models, develop optimal controls, and reconstruct systems with desired properties.In this paper, we use a novel approach to build up abstract mathematical models at the molecular level in Section2.2, based directly on biological theory. With reasonable assumptions, we can avoid the conventional obstacles mentioned above. Different from existing methods, focusing on a gene or changes in chemical concentration, we emphasize the base change in the nucleotide bases. The cost function, a summation of costs for applying mutagens and the off-trajectory penalty, together with the system equations, formulates the optimal control problem. The optimal control is then solved by dynamic programming algorithm in Section 2.3. Section 3 shows simulation results of the optimal control problem at different scales and is followed by several important propositions.


## 2. System Equations and Generalized Optimal Control Problem Formulation

The central dogma of molecular biology, first elaborated in [23] and restated in [24], illustrates the detailed residue-by-residue transfer of genetic sequential information. Nowadays, it is widely recognized as the backbone of molecular biology. It describes the genetic information flow among three kinds of biopolymers: DNA, RNA, and protein. In most living organisms, genetic information transfers from DNA to RNA, and then to protein. This process is usually irreversible, thus protein always acts as the sink of information flow. A codon consists of three consecutive nucleotide bases, corresponding to one amino acid according to the genetic codes. Since there are only 20 kinds of amino acids and 64 combinations of codons, there exists redundancy in genetic codes.We are particularly interested in mutations that happen during the process of DNA replication, as DNA serves as long-term genetic information storage and is the basis of genetic inheritance, the accuracy of which is particularly important to ensure the correct expression of genes. DNA molecules consist of four kinds of nucleotide acids,adenine (A), thymine (T), guanine (G), and cytosine (C), and a backbone made of sugars and phosphate. In 1953, James D. Watson and Francis Crick found the double helix structure of DNA and the rule of basepairing, known as Watson-Crick basepairing [25, 26]. A always pairs with T, G always pairs with C, and vice versa. In nature, replication errors occur at a very low rate, one error for every 107 nucleotides added [27]. The redundancy of information caused by the double-helix structure ensures the accuracy of DNA replication. Some DNA self-repair mechanisms, listed in [28], such as proofreading, also help to eliminate errors during the replication process.Gene mutations are changes in the nucleotide sequence of DNA or RNA. Usually, we focus only on mutations occurring in coding DNA sequences and RNA. Mutations are caused by various reasons. Induced mutations are caused by either chemical mutagens or radiation. In general, radiation induces higher randomness than chemical mutagens. Point mutation is the simplest form of mutation, involving only one base. Point mutations can be further divided into transitions (A↔G or C↔T) and transversions (A/G↔C/T). Transversions are theoretically expected to be twice as frequent as transitions, but transitions may be favored over transversions in coding DNA because they usually result in a more conserved polypeptide sequence [29].In this section, we first give the problem statement in Section2.1, and then we construct system equations for both deterministic and stochastic mutations in Section 2.2. At last, in Section 2.3, we formulate the optimal control problem and apply dynamic programming algorithm to solve it.


### 2.1. Problem Statement

Figure3 shows the system diagram of restoring an abnormal DNA segment back to a normal sequence by applying mutagens during the process of DNA replication. After we obtain a patient’s genome, we compare the coding DNA segments with normal DNA segments in our database to figure out the possible range of mutated segments. Due to the redundancy in genetic codes, as long as any two DNA segments can be transcribed and then translated to the same amino acid sequence, the distance reference between them is considered to be zero. Therefore, instead of having a predetermined final state or a neighborhood of a final state, our final state lies in a set where the distance reference between any sequence in this set and the desired sequence is zero. We name this set the final desired set. The prescription is then determined by comparing the current measurement and every sequence in the desired set. Both internal noises and external disturbance can be eliminated by the measurement. We treat DNA sequences as state variables, the ON/OFF controls of all available mutagens at every spot on the given DNA segment as inputs, the measurements as the outputs, and one cell cycle as the step increment in our system equations.Figure 3
System diagram of restoring an abnormal DNA segment back to a normal sequence by applying mutagens during the process of DNA replication.The objective function is defined as a summation of the costs (including risks) of applying mutagens and the off-trajectory penalty. The optimal control sequences are computed beforehand to let doctors make treatment plans according to the patient’s condition. In general, the optimal control sequence and the corresponding optimal trajectory are not unique because the bases mutate independently in most cases and the order of mutating different bases does not matter if the number of medical treatment sessions is not under a tight restriction. Additional measurements are taken before and after each treatment, if necessary. In deterministic cases, the purpose of taking additional measurements is to check the current sequence and to eliminate both internal and external disturbances. The treatment plan is adjusted if the measurement is not the same as expected. In stochastic settings, the measurements are taken to conquer the randomness caused by both mutagens and other noises. The treatment is then updated accordingly.To sum up, our system is a discrete-time dynamic system with finite state space and output space, and a set of ON/OFF switches as controls. Our goal is to optimally drive this system from a given initial state to a desired final set at the lowest cost.


### 2.2. System Equations Formulation

We mainly focus on applying chemical mutagens and radiation to restore the original amino acid sequence during the process of DNA replication. Other factors that may affect the gene mutation, including temperature and electroporation, are not within our consideration. In addition, we assume that chemical mutagens or radiation target one and only one nucleotide base at any preset site, despite the technical limitation, and the results of applying chemical mutagens and radiation are independent. For simplicity, we normalize the dose of mutagens to transfer one nucleotide base to another in one step to 1. In most cases, nucleotide bases mutate independently, therefore there is no chain effect caused by mutagens. To avoid reactions among different mutagens, we require that at most one chemical mutagen and one radiation be applied in each cell cycle. While constructing a generalized model, since the order of applying chemical mutagens and radiation does not affect the results, without loss of generality, we require they be applied in the order shown in Figure4. That is, chemical mutagens are always applied before the duplication process starts, radiation is always applied in the middle of the cell cycle, and the measurements are taken before every replication starts. Lastly, we assume that the measurements are always correct, and DNA replication error, background mutation rate, and other random noise can be eliminated from measurements by considering them as spontaneous mutation.Figure 4
The order of taking measurements and applying chemical mutagens and radiation in a cell cycle.


#### 2.2.1. Base-to-Base Deterministic Mutations

Denote the targeted DNA segment withn nucleotide bases at kth step by a column vector xk, as shown in Figure 5. xki is the ith element of xk. Let P be the transfer matrix from xk to xk+1, for all k, k∈ℤ+∪{0}, without mutation. Then, the perfect DNA replication process can be expressed as(1)xk+1=Pxk.Figure 5
Biological information flow in central dogma of molecular biology.Proposition 1.
P=-I.Proof.
As no mutation occurs,xk+1 is completely complementary to xk by Watson-Crick base pairing rule, and xk+2 is completely complementary to xk+1. Therefore, xk+2 is exactly the same as xk. thus,
(2)xk+2=Pxk+1=P2xk⟹P2=I.
Since every base mutates independently, every element of xk+1i only depends on the corresponding element of xki, thus P is diagonal. In addition, xk+1≠xk, we conclude P=-I.Based on Proposition1, we assign values to nucleotide bases set {A,G,C,T,O}, where O is an artificial nonsense base. Define an equivalence relationship between {A,G,C,T,O} and {1,2,-2,-1,0}, that is, {A,G,C,T,O}⇔{1,2,-2,-1,0}, with(3)xki={1,ifA,2,ifG,-2,ifC,-1,ifT,0,ifO.Proposition 2.
{1,2,-2,-1,0} is a field under proper definitions of addition and multiplication.Proof.
Defining the addition table and multiplication table as in Tables1 and 2, we check if the set {1,2,-2,-1,0} satisfies the definition of field.Table 1
Addition table for{1,2,-2,-1,0}.
+12−2−1012−2−1012−2−1012−2−1012−2−1012−2−1012−2−10Table 2
Multiplication table for{1,2,-2,-1,0}.
×12−2−10112−2−1022−11−20−2−21−120−1−1−2210000000Closed under Addition and Multiplication
Satisfied obviously from Tables1 and 2.Associativity of Addition and Multiplication
Implicitly satisfied by integer addition and multiplication.Commutativity of Addition and Multiplication
Satisfied as Tables1 and 2 are symmetric according to the diagonal.Additive and Multiplicative Identity
Additive identity is 0, and multiplicative identity is 1.Additive and Multiplicative Inverses
Additive inverses pair:1↔-1, 2↔-2, 0↔0.
Multiplicative inverses pair:1↔1, 2↔-2, -1↔-1.Distributivity of Multiplication over Addition
Implicitly satisfied by integer addition and multiplication.We conclude{0,1,2,-2,-1} is a field under addition and multiplication defined by Tables 1 and 2.From now on, we useℱ to denote the field {0,1,2,-2,-1}. And xk∈ℱn is the state vector representing a DNA segment with n nucleotide bases, where ℱn is the set of ℱ-valued vectors of dimension n.We start with the simplest form of mutations, point mutation. Suppose there is a point mutation, we write it mathematically as(4)xk+1=(-I+Δs)xk+Δw,
where xk+1,xk∈ℱ, and -I reduces to −1 as only one base is involved. The corresponding values of Δs and Δw, obtained by reverse engineering with all possible pairs of xk and xk+1, are listed in Table 3.Table 3
Possible values ofΔs and Δw. The corresponding values of Δs and Δw are obtained by substituting the value of xk and xk+1 into (4), with Δw≠0 only if xk=0.
kth(k+1)thAGCTOA2−2−101111}ΔsG−120−2C−202−1T0−1−22O12−2−10}ΔwHere,Δs represents the mutation from four normal nucleotide bases, and Δw corresponds to mutation from nonsense base, that is, Δw≠0 only if xk=0.Rewriting (4) by collecting all values of Δs and Δw in Table 3, we get(5a)xk+1=(-I+∑j=04ukjsj)xk+∑j=04ckjwj(5b)=(-I+uks)xk+ckw,where {s0,s1,s2,s3,s4}={w0,w1,w2,w3,w4}={0,1,2,-2,-1}, ukj,ckj∈{0,1}, representing the on/off controls, uk=[uk0uk1uk2uk3uk4], ck=[ck0ck1ck2ck3ck4], and s=w=[012-2-1]T. In (5a), sj and wj are constants for all k and j. ukj and ckj, the inputs of the system, are the on/off controls for chemical mutagens or radiation. Clearly, ∑j=04ckj=1 only if xk=0. Equation (5b) is a simplified version of (5a) as we put ukj,sj,ckj,wj into vector form uk,s,ck,w. s and w serve as vector basis for base-to-base deterministic model. uk and ck are now multi-input controls; each of them contains 5 on/off controls, corresponding to all possible transfer patterns. For a particular k, at most one of ukjs and ckjs can be 1, as stated in Proposition 3. This is consistent with the fact that every state can be transferred to only one of the five states in the state space ℱ with corresponding mutagens available.Proposition 3.
It is always1-1 transfer when mutation occurs, that is, one nucleotide base can only transfer to another one, therefore(i)
ifxk=0 and ck=0, or ck=[10000], then xk+1=0,(ii)
ifxk≠0, then ck=0 and uk is either 0 or a unit row vector,(iii)
ifxk=0, then uk=0 and ck is either 0 or a unit row vector,(iv)
uk+ck is either 0 or a unit row vector, for all k∈ℤ+∪{0}.Now, suppose for some reason we need to take an addition, measurement in the middle of the cell cycle, after the completion of thekth duplication and before the start of the (k+1)th. We name this kind of measurement an intermediate state, and denote by it xk′. Then, we have(6)xk+1=(I+Δs′)xk′+Δw′,
where the values of Δs′ and Δw′, listed in Table 4, are obtained in the same way as getting Δs and Δw in Table 3.Table 4
Possible values ofΔs′ and Δw′. The corresponding values of Δs′ and Δw′ are obtained by substituting the value of xk′ and xk+1 into (4), with Δw′≠0 only if xk′=0.
kth(k+1)thAGCTOA012−2-1-1-1-1}Δs′G20−21C1−202T−2210O12−2−10}Δw′Comparing Tables3 and 4, we find the collection of Δs and Δs′, Δw and Δw′, form the same set, respectively. Thus, we continue using s and w when rewriting (6) in the form of (5a) and (5b), that is,
(7)xk+1=(I+vks)xk′+ck′w,
where vk,ck′ are the counterparts of uk,ck, respectively, and s,w are the same as in (5b).Similar to Proposition3, we get Proposition 4.Proposition 4.
vk and ck′ in (7) need to satisfy the following conditions.(i)
Ifxk=0 and ck′=0, or ck′=[10000], then xk+1=0.(ii)
Ifxk≠0, then ck′=0 and vk is either 0 or a unit row vector.(iii)
Ifxk=0, then vk=0 and ck is either 0 or a unit row vector.(iv)
vk+ck′ is either 0 or a unit row vector, for all k∈ℤ+∪{0}.Now take, both chemical mutagens and radiative rays under our consideration and apply them in the order as shown in Figure4. Then, we can express our system equation as(8a)xk′=(-I+uks︸mutationscausedbychemicalmutagensfromnormalbases)xk+ckw︸mutationscausedbychemicalmutagensfromO,(8b)xk+1=(I+vks︸mutationscausedbyradiativeraysfromnormalbases)xk′+ck′w︸mutationscausedbyradiativeraysfromO,(8c)yk=xk,where uk and vk are the inputs of the system and yk is the measurement. Obviously, (8a) is modified from (5b) and (8b) from (7).The two-step mutation and the intermediate statexk′avoid the case xk is changed to different bases by radiation and chemical mutagens simultaneously, which causes confusion. Substituting (8a) and (8b), we get(9a)xk+1=(I+vks)(-I+uks)xk+(I+vks)ckw+ck′w,(9b)yk=xk.Obviously, Proposition 3 still holds for uk and ck, and Proposition 4 holds for vk and ck′ for (9a).For point mutations, we have 20 on/off controls in total for every stepk, 10 for chemical mutagens as described before, and the rest for radiation.


#### 2.2.2. Gene-to-Gene Deterministic Mutations

In general, mutations involve multiple bases. Therefore, large-scale deterministic model is necessary. Now, we show how to extend our model to large-scale systems.Suppose we have a coding DNA segment with lengthn, then xk∈ℱn. Since a coding DNA segment usually contains integer number of codons, which is made of three consecutive bases, n is a multiple of 3. Let xki denote the ith component of xk. This notation is consistent with the one in Section 2.2.1. Initiated by the base-to-base deterministic model from Section 2.2.1, we write our system equation for large-scale system as
(10)xk′=(-I+∑i=1nukiSki︸mutationscausedbychemicalmutagensfromnormalbases)xk+∑i∈OkckiWki︸mutationscausedbychemicalmutagensfromO,xk+1=(I+∑i=1nvkiQki︸mutationscausedbyradiativeraysfromnormalbases)xk′+∑i∈Ok′bkiRki︸mutationscausedbyradiativeraysfromO,yk=xk,
where uki,vki,cki,bki are on/off controls of the ith element,Ski,Qki are n×n square matrices corresponding to the mutations between normal bases or from normal bases by chemicals and radiation, respectively, Wki,Rki are n-dimensional column vectors representing mutations from nonsense bases by chemicals and radiation, respectively, and 𝒪k={i:xki=0,1≤i≤n}, 𝒪k′={i:xk′i=0,1≤i≤n}.Ski and Qki are diagonal matrices since each base mutates independently. The values in the first four rows of Tables 3 and 4 correspond to the diagonal elements of Ski and Qki, respectively. The last rows of Tables 3 and 4 are assigned to Wki and Rki, n-dimensional vectors, at nonsense base’s spots for xk.Define𝒮={sjeieiT,∀i,j,0≤j≤4,1≤i≤n}, a collection of n×n matrices, where sj is the same as in (5a) and (5b), ei is the unit column vector of length n with ith component equal to 1 and all other components equal to 0, and eieiT is the square matrix with only the ith element on the diagonal equals to 1, and 0 otherwise. Then, Ski,Qki can be written as linear combinations of all elements from 𝒮, with the coefficient of each element either 0 or 1 corresponding to the on/off control uk(i,j) and vk(i,j), respectively.Similarly, define𝒲={wjei,∀i,j,0≤j≤4,1≤i≤n}, where wj is the same as (5a) and (5b). Wki,Rki can be written as linear combinations of all components from 𝒲, with coefficient of every component either 0 or 1 corresponding to the on/off control ck(i,j) and ck′(i,j), respectively.Therefore, instead of using step-varyingSki,Ski,Wki,Rki, we find matrix basis for those four square matrices to make the controls to be the only variables depending on k, as we did for single-base cases. Then, we can, write (10) as(11a)xk′=(-I+∑i=1n∑j=04uk(i,j)sjeieiT)xk+∑i∈Ok∑j=04ck(i,j)wjei,(11b)xk+1=(I+∑i=1n∑j=04vk(i,j)sjeieiT)xk′+∑i∈Ok′∑j=04ck′(i,j)wjei,(11c)yk=xk,where uk(i,j),vk(i,j),ck(i,j),ck′(i,j)∈{0,1}.As shown in (11a), (11b), and (11c), multisites mutations contain 20n controls in total for every step k, where n is the number of nucleotide bases on the targeted gene. Similar to point mutations, every single site has 20 controls in each step, 10 for chemical mutagens and 10 for radiation.We can viewuk,vk,ck,ck′ as binary matrices of dimension n×5, and uk(i,j), vk(i,j)ck(i,j), ck′(i,j)are the corresponding element of ith row and jth column. Use uki,vki,cki,bki, binary row vectors of dimension 5, to denote the ith row of uk,vk,ck,ck′, respectively. Again, s=w=[012-2-1]T.Combining (11a) and (11b), and writing control variables in vector forms, we get(12)xk+1=(I+∑i=1nvkiseieiT)(-I+∑i=1nukiseieiT)xk+(I+∑i=1nvkiseieiT)∑i∈Okckiwei+∑i∈Ok′bkiwei,yk=xk.Proposition 5.
For large-scale deterministic system,uk,vk,ck,ck′ satisfy conditions below.(i)
IfeiTxk=0, then i∈𝒪k.(ii)
IfeiTxk=0, cki=0 or cki=[10000], then i∈𝒪k′.(iii)
For alli∉𝒪k, uki is either 0 or a row unit vector and cki=0.(iv)
For alli∈𝒪k, cki is either 0 a row unit vector and uki=0.(v)
For alli∉𝒪k′, vki is either 0 or a row unit vector and bki=0.(vi)
For alli∈𝒪k′, bki is either 0 or a row unit vector and vki=0.(vii)
For alli,k, 1≤i≤n,  k∈ℤ+∪{0}, uki+cki is either 0 or a unit row vector and vki+bki is either 0 or a unit row vector.The mathematical model (12) is quite flexible and can be easily extended to many cases, such as transcription process, multiple spot mutations within one-step or broken DNA strands.Take broken DNA strands as an example. DNA strand breaks due to various reasons. Our system equation can represent this phenomenon by dividing the whole system into small subsystems. Significant brokage of DNA strands is simply eliminated by cell mechanism to ensure the accuracy to DNA replication. Equation (13) shows the case of one single DNA strand breaking into two segments by chemical mutagens (13)(xk′(1)xk′(2))=(-Im+∑i=1mukiseieiT00-In-m+∑i=m+1nukiseieiT)×(xk(1)xk(2))+(∑i∈Ok,1≤i≤mckiwei∑i∈Ok,(m+1)≤i≤nckiwei),xk+1(1)=(Im+∑i=1mvkiseieiT)xk′(1)+∑i∈Ok′,1≤i≤mbkiwei,xk+1(2)=(In-m+∑i=m+1nvkiseieiT)xk′(2)+∑i∈Ok′,(m+1)≤i≤nbkiwei.


#### 2.2.3. Gene-to-Gene Stochastic Mutations

In reality, mutagens, no matter chemical or radiative, always cause randomness in mutation. Therefore, we need to derive the model for gene-to-gene stochastic mutations.Introduce new random variables,hk,l1(i,j),rk,l2(i,j),hk,l3′(i,j),rk,l4′(i,j)∈{0,1}, associated with probability pl1,j(h), pl2,j(r),pl3,j(h′),pl4,j(r), for all i,k, 1≤i≤n,  k∈ℤ+∪{0}, respectively, where k is the step index, l1,l2 are indices for chemical mutagens inducing mutation from normal bases and from O, respectively, l3,l4 are indices for radiation inducing mutation from normal bases and from O, respectively, i is the index of DNA segment, and the value of j corresponds to the transfer pattern, which can be found in Tables 3 and 4. Note different mutagens have different probability assignments, the probability assignments are only related to the type of mutagens, and the probability associated with every kind of mutagens sums up to 1, that is,
(14)∑j=04pl1,j(h)=1,∀l1,1≤l1≤l,∑j=04pl2,j(r)=1,∀l2,1≤l2≤m,∑j=04pl3,j(h′)=1,∀l3,1≤l3≤l′,∑j=04pl4,j(r′)=1,∀l4,1≤l4≤m′.
The controls are uk,l1i,ck,l2i,vk,l3i,bk,l4i∈{0,1}, with the fact that 1 representing mutagen with corresponding index is applied at ith spot of DNA segment at kth generation, and 0 representing mutagen with corresponding index is not applied at spot i at kth step, similar to Sections 2.2.1 and 2.2.2. The mutagen indices l1,l2,l3,l4 can be omitted in deterministic mutations since given the current state and control, the next state is unique. However, they are necessary for stochastic mutations, because there exist multiple possible states for the next stage given the control. In other words, the next state is determined by random variables hk,l1(i,j),rk,l2(i,j),hk,l3′(i,j),rk,l4′(i,j), given the values of uk,l1i,ck,l2i,vk,l3i,bk,l4i, and xk.Suppose we have(l+m) kinds of chemical mutagens available, with l kinds to induce mutations from normal bases and m kinds to induce mutations from O. And we have (l′+m′) kinds of radiation available, with l′ kinds to induce mutations from normal bases and m′ kinds to induce mutations from O. Therefore, we have total (l+m+l′+m′) controls for each spot i at each step k. We can write our system equation as(15a)xk′=(-I+∑l1=1l∑i=1nuk,l1i∑j=04hk,l1(i,j)sjeieiT︸mutationscausedbychemicalmutagensfromnormalbases)xk+∑l2=1m∑i∈Okck,l2i∑j=04rk,l2(i,j)wjei︸mutationscausedbychemicalmutagensfromO,(15b)xk+1=(I+∑l3=1l′∑i=1nvk,l3i∑j=04hk,l3′(i,j)sjeieiT︸mutationscausedbyradiativeraysfromnormalbases)xk′+∑l4=1m′∑i∈Ok′bk,l4i∑j=04rk,l4′(i,j)wjei︸mutationscausedbyradiativeraysfromO,(15c)yk=xk.Again, we define hk,l1(i,j),rk,l2(i,j),hk,l3′(i,j),rk,l4′(i,j) the elements at ith row and jth column of n×5 binary matrices hk,l1,rk,l2,hk,l3′,rk,l4′, respectively. hk,l1i,rk,l2i,hk,l3′i, rk,l4′i, and binary row vectors of dimension 5 denote the ith row of hk,l1,rk,l2,hk,l3′,rk,l4′, respectively. Then, we can simplify (15a), (15b), and (15c) and combine (15a) and (15b) as(16)xk+1=(I+∑l3=1l′∑i=1nvk,l3ihk,l3′iseieiT)×(-I+∑l1=1l∑i=1nuk,l1ihk,l1iseieiT)xk+(I+∑l3=1l′∑i=1nvk,l3ihk,l3′iseieiT)∑l2=1m∑i∈Okck,il2rk,l2iwei+∑l4=1m′∑i∈Ok′bk,l4irk,l4′iwei,yk=xk.Proposition 6.
For large-scale stochastic system,uk,l1i,hk,l1i,ck,l2i,rk,l2i,vk,l3i,hk,l3′i,bk,l4i,rk,l4′i follow the rules below.(i)
IfeiTxk=0, then i∈𝒪k.(ii)
IfeiTxk=0 and ∑l2=1mck,l2i=0, then i∈𝒪k′.(iii)
IfeiTxk=0, ∑l2=1mck,l2i=1 and rk,l2i=[10000], then i∈𝒪k′.(iv)
For alli,k,l1, 1≤i≤n,  k∈ℤ+∪{0},1≤l1≤l, if uk,l1i=1, then hk,l1i is a unit row vector.(v)
For alli,k,l2, 1≤i≤n,  k∈ℤ+∪{0}, 1≤l2≤m, if ck,l2i=1, then rk,l2i is a unit row vector.(vi)
For alli,k,l3, 1≤i≤n,  k∈ℤ+∪{0}, 1≤l3≤l′, if vk,l3i=1, then hk,l3′i is a unit row vector.(vii)
For alli,k,l4, 1≤i≤n,  k∈ℤ+∪{0}, 1≤l4≤m′, if bk,l4i=1, then rk,l4′i is a unit row vector.(viii)
For alli∉𝒪k, ∑l1=1luk,l1i=0or1 and ck,l2i=0, for all l2,1≤l2≤m.(ix)
For alli∈𝒪k, ∑l2=1mck,l2i=0 or 1 and uk,l1i=0, for all l1, 1≤l1≤l.(x)
For alli∉𝒪k′, ∑l3=1l′vk,l3i=0 or 1 and bk,l4i=0, for alll4, 1≤l4≤m′.(xi)
For alli∈𝒪k′, ∑l4=1m′bk,l4i=0 or 1 and vk,l3i=0, for all l3, 1≤l3≤l′′.(xii)
For alli,k, 1≤i≤n,  k∈ℤ+∪{0}, ∑l1=1luk,l1i+∑l2=1mck,l2i=0 or 1 and ∑l3=1l′vk,l3i+∑l4=1m′bk,l4i=0 or 1.We close this section with the definition of controllability to the system equations proposed above. DNA replication systems with system equations proposed as (9a), (9b), (12), and (16) are completely controllable if and only if  for all x0,x2k1,x2k2+1∈ℱ, k1,k2∈ℤ+∪{0},∃ at least one path from x0 to x2k1 and at least one path from x0 to x2k2+1 by applying proper mutagens in the correct order, with k1,k2 finite.


### 2.3. Generalized Optimal Control Problem Formulation

We first define our objective function that can be adapted to all kinds of systems proposed in Section2.2 with minor changes. Mathematically, in systems where the controllable parameters of interest are discrete, the objective function is usually a weighted sum representing the number of times that a piece of equipment is turned “on” or “off” or the number of resources needed to execute certain tasks in the frequent cases [30]. In our case, this summation is the total number of times that different mutagens are applied weighted by the corresponding cost (including the risk). Another key factor of objective function is the off-trajectory penalty. Designing a trajectory beforehand is necessary to avoid other hidden risks. If the measurement indicates that the current state is off the predefined trajectory, we include a distance reference between current state and desired state as penalty and change the treatment plan accordingly.Therefore, our objective function can be expressed as(17)J0(x0)=min⁡u,c,v,c′Eh,h′,r,r′[∑k=0N-1∑l1=1l∑i=1nαl1uk,l1i+∑k=0N-1∑l2=1m∑i=1nβl2ck,l2i︸costofapplyingchemicalmutagens+∑k=0N-1∑l3=1l′∑i=1nαl3′vk,l3i+∑k=0N-1∑l4=1m′∑i=1nβl4′bk,l4i︸costofapplyingradiativerays+∑k=0Nd(xk,{xkd})︸tracingcost],
with x0,xkd∈ℱn, 1≤k≤N, n≡0 (mod3) given. The physical meaning of uk,l1i, ck,l2i, vk,l3i, bk,l4i, l1, l2, l3, l4 is the same as in Section 2.2.3. αl1,βl2,αl3′,βl4′∈ℝ, for alll1,l2,l3,l4, 1≤l1≤l,  1≤l2≤m, 1≤l3≤l′,  1≤l4≤m′, are the corresponding cost of applying chemical mutagens and radiative rays indexed l1,l2,l3,l4, respectively. {xkd}:ℱn×ℱn→ℝ+∪{0} denotes the desired set at kth stage, generated by the DNA sequences representing the same amino acid sequence asxkd, the desired state at kth step. And d(xk,{xkd}) is the distance reference of the current state xk to the desired set {xkd} at kth step. The final penalty, the distance reference from the final state to the desired set at k=N, is included in the last term.In general,βl2,βl4′≪αl1,αl3′, for all l1,l2,l3,l4,  1≤l1≤l,  1≤l2≤m,  1≤l3≤l′,  1≤l4≤m′, because physically O is a set of nonsense bases and more details are necessary to convert an O back to normal bases, for instance, the cost to identify the exact element in the set O. Our goal is to drive our system optimally from initial state x0 to the desired final set {xNd} by applying a sequence of mutagens indexed with {l1,l2,l3,l4}, at problematic positions i, and in a correct order k.In (17), the first four terms inside the expectation do not depend on random variables hk,l1i, rk,l2i, hk,l3′i, and rk,l4′i, for all i,k,l1,l2,l3,l4 as the treatment plan is computed based on the initial state x0. Given yk, the updated treatment plan is computed accordingly but still not related to random variables. The last term inside expectation, ∑k=0Nd(xk,{xkd}), is the only term in summation that depends on the distribution of the random variables.The constraint of the optimal control problem, in general, is the system equation. We choose multidimensional stochastic system equation as the generalized constraints as it can be degenerated to one-dimensional and multidimensional deterministic cases by proper modifications.Therefore, we can rewrite our objective function and formulate our optimal control problem as(18)J0(x0)=min⁡{u,c,v,c′}0,1,…,N-1[∑k=0N-1∑l1=1l′∑i=1nαl1uk,l1i+∑k=0N-1∑l2=1m∑i=1nβl2ck,l2i+∑k=0N-1∑l3=1l′∑i=1nαl3′vk,l3i+∑k=0N-1∑l4=1m′∑i=1nβl4′bk,l4i+∑k=0NE{h,r,h′,r′}0,1,…,N-1[d(xk,{xkd})]],
subject to(19)xk+1=(I+∑l3=1l′∑i=1nvk,l3ihk,l3′iseieiT)×(-I+∑l1=1l∑i=1nuk,l1ihk,l1iseieiT)xk+(I+∑l3=1l′∑i=1nvk,l3ihk,l3′iseieiT)∑l2=1m∑i∈Okck,l2irk,l2iwei+∑l4=1m′∑i∈Ok′bk,l4irk,l4′iwei,yk=xk.
We need to choose a proper distance reference to quantitatively describe the relationship between DNA segments of same length. We first define the distance reference between codons, and the distance reference between DNA segments is a weighted sum of distance reference between every pair of codons.The distance reference between codons,d(φ1,φ2),φ1,φ2∈ℱ3, needs to fulfill the biological requirements as below.(1)
Nonnegativity: the distance reference between any two codons is either positive or zero. Mathematically,d:ℱ3×ℱ3→ℝ+∪{0}, d(φ1,φ2)≥0.(2)
The distance reference between two codons corresponding to the same amino acid is zero.(3)
Symmetry: the distance reference from codonφ1 to codon φ2 equals the distance reference from codon φ2 to codon φ1, that is, d(φ1,φ2)=d(φ2,φ1).(4)
The distance reference between two codons corresponding to different amino acids should reveal the chemical and physical differences between two amino acids.(5)
The distance reference from stop codons to all other codons is much larger than those between other codons as early termination of amino acid sequences is more harmful than other forms of mutations.All the existing metric defined on the finite field cannot achieve all the requirements above. The second requirement violates the identity of indiscernible, that is,d(φ1,φ2)=0 if and only if φ1=φ2. The redundancy in genetic codes implies d(φ1,φ2)=0 if those two amino acids, φ1 and φ2, are translated into the same amino acids. In addition, the triangular inequality is not necessarily true, according to the underlying physical meanings. We take the assumption that the stop codons are of the same distance reference from and to all other codons.Important physical and chemical properties are listed in Table5. We ignore codons containing O since their chemical and physical properties cannot be found in literature.Table 5
Properties of amino acids.
Amino AcidAbbrev.Codon(s)PolarityPHSizeAlanineAlaGCT,GCC,GCA,GCGNonpolar6.01TinyArginineArgCGA,CGG,CGC,CGT,
AGA,AGGPolar10.76NormalAsparagineAsnAAC,AATPolar5.41SmallAspartic acidAspGAT,GACPolar2.85SmallCysteineCysTGT,TGCNonpolar5.05SmallGlutamineGlnCAA,CAGPolar5.65NormalGlutamic acidGluGAA,GAGPolar3.15NormalGlycineGlyGGA,GGG,GGC,GGTNonpolar6.06TinyHistidineHisCAC,CATPolar7.60NormalIsoleucineIleATA,ATC,ATTNonpolar6.05NormalLeucineLeuTTA,TTG,CTA,CTG,
CTC,CTTNonpolar6.01NormalLysineLysAAA,AAGPolar9.60NormalMethionineMetATGNonpolar5.74NormalPhenylalaninePheTTC,TTTNonpolar5.49NormalProlineProCCA,CCG,CCC,CCTNonpolar6.30SmallSerineSerTCA,TCG,TCC,TCT,Polar5.68TinyAGT,AGCThreonineThrACT,ACC,ACA,ACGPolar5.60SmallTryptophanTrpTCCNonpolar5.89NormalTyrosineTyrTAC,TATPolar5.64NormalValineValGTA,GTG,GTC,GTTNonpolar6.00SmallStop codonTermTAA,TAG,TGA———From Table5, we can see all codons are divided into different sets with each set corresponding to one amino acid. The size and the elements in one codon set vary from one amino acid to another. This implies that the costs of driving one codon to the desired final set generated by the desired final state might be different from the costs of driving the complementary codon to the desired final set generated by the complementary of desired final state. More discussions about this issue are presented in Section 3.The distance reference between any two codons can be defined by a weighted sum of the differences between physical and chemical properties or other reasonable functions. And the distance reference between two DNA sequences is defined as the sum of distance reference between the corresponding pair of codons. The biological statics plays a crucial rule to define this distance function in practical.An example of the distance function can be expressed as(20)d(ξ1,ξ2)=ζpolaritypolarity(ξ1,ξ2)+ζPHPH(ξ1,ξ2)+ζsizesize(ξ1,ξ2),polarity(ξ1,ξ2)={0ifξ1,ξ2arebothpolarornon-polar,1ifoneofξ1,ξ2ispolar,andtheothernon-polar,PH(ξ1,ξ2)=|PHvalueofξ1-PHvalueofξ2|,size(ξ1,ξ2)={0,ifξ1,ξ2arebothtiny,small,ornormal,σ1,ifoneofξ1,ξ2istiny,andtheothersmall,σ2,ifoneofξ1,ξ2istiny,andtheothernormal,σ3,ifoneofξ1,ξ2issmall,andtheothernormal,
where ξ1,ξ2 are two amino acids.d(ξ1,ξ2) is then assigned to d(φ1,φ2) with φ1,φ2 corresponding to amino acids ξ1,ξ2, respectively.Since the generalized optimal control problem in (18) and (19) is a multistage problem that can be broken down into simpler steps at different time points. Therefore, we can solve it by dynamic programming.For dynamic programming, the optimal control policy is constructed backward. And Bellman’s principle of optimality states that the optimal policy forx0 to {xNd} is also the optimal policy for the tail problem, from xq to {xNd}.The tail problem is defined as(21)Jq(xq)=min⁡{u,c,v,c′}q,q+1,…,N-1{∑k=qN-1∑l1=1l∑i=1nαl1uk,l1i+∑k=qN-1∑l2=1m∑i=1nβl2ck,l2i+∑k=qN-1∑l3=1l′∑i=1nαl3′vk,l3i+∑k=qN-1∑l4=1m′∑i=1nβl4′bk,l4i+∑k=qNE{h,r,h′,r′}q,q+1,…,N-1[d(xk,{xkd})]}.
The iterative update equation to find optimal policy can be expressed by (22), according to the dynamic programming algorithm in [31]. 								(22)JN(xN)=d(xN,{xNd}),Jq(xq)=min⁡uq,cq,vq,cq′Ehq,rq,hq′,rq′[∑l1=1l∑i=1nαl1uq,l1i+∑l2=1m∑i=1nβl2cq,l2i+∑l3=1l′∑i=1nαl3′vq,l3i+∑l4=1m′∑i=1nβl4′cq,l4′i+d(xq,{xqd})+Jq+1(xq+1)]=min⁡uq,cq,vq,cq′{∑l1=1l∑i=1nαl1uq,l1i+∑l2=1m∑i=1nβl2cq,l2i+∑l3=1l′∑i=1nαl3′vq,l3i+∑l4=1m′∑i=1nβl4′cq,l4′i+Ehq,rq,hq′,rq′[d(xq,{xqd})+Jq+1(xq+1)]∑l1=1l∑i=1nαl1uq,l1i},q=0,1,…,N-1.


## 2.1. Problem Statement

Figure3 shows the system diagram of restoring an abnormal DNA segment back to a normal sequence by applying mutagens during the process of DNA replication. After we obtain a patient’s genome, we compare the coding DNA segments with normal DNA segments in our database to figure out the possible range of mutated segments. Due to the redundancy in genetic codes, as long as any two DNA segments can be transcribed and then translated to the same amino acid sequence, the distance reference between them is considered to be zero. Therefore, instead of having a predetermined final state or a neighborhood of a final state, our final state lies in a set where the distance reference between any sequence in this set and the desired sequence is zero. We name this set the final desired set. The prescription is then determined by comparing the current measurement and every sequence in the desired set. Both internal noises and external disturbance can be eliminated by the measurement. We treat DNA sequences as state variables, the ON/OFF controls of all available mutagens at every spot on the given DNA segment as inputs, the measurements as the outputs, and one cell cycle as the step increment in our system equations.Figure 3
System diagram of restoring an abnormal DNA segment back to a normal sequence by applying mutagens during the process of DNA replication.The objective function is defined as a summation of the costs (including risks) of applying mutagens and the off-trajectory penalty. The optimal control sequences are computed beforehand to let doctors make treatment plans according to the patient’s condition. In general, the optimal control sequence and the corresponding optimal trajectory are not unique because the bases mutate independently in most cases and the order of mutating different bases does not matter if the number of medical treatment sessions is not under a tight restriction. Additional measurements are taken before and after each treatment, if necessary. In deterministic cases, the purpose of taking additional measurements is to check the current sequence and to eliminate both internal and external disturbances. The treatment plan is adjusted if the measurement is not the same as expected. In stochastic settings, the measurements are taken to conquer the randomness caused by both mutagens and other noises. The treatment is then updated accordingly.To sum up, our system is a discrete-time dynamic system with finite state space and output space, and a set of ON/OFF switches as controls. Our goal is to optimally drive this system from a given initial state to a desired final set at the lowest cost.


## 2.2. System Equations Formulation

We mainly focus on applying chemical mutagens and radiation to restore the original amino acid sequence during the process of DNA replication. Other factors that may affect the gene mutation, including temperature and electroporation, are not within our consideration. In addition, we assume that chemical mutagens or radiation target one and only one nucleotide base at any preset site, despite the technical limitation, and the results of applying chemical mutagens and radiation are independent. For simplicity, we normalize the dose of mutagens to transfer one nucleotide base to another in one step to 1. In most cases, nucleotide bases mutate independently, therefore there is no chain effect caused by mutagens. To avoid reactions among different mutagens, we require that at most one chemical mutagen and one radiation be applied in each cell cycle. While constructing a generalized model, since the order of applying chemical mutagens and radiation does not affect the results, without loss of generality, we require they be applied in the order shown in Figure4. That is, chemical mutagens are always applied before the duplication process starts, radiation is always applied in the middle of the cell cycle, and the measurements are taken before every replication starts. Lastly, we assume that the measurements are always correct, and DNA replication error, background mutation rate, and other random noise can be eliminated from measurements by considering them as spontaneous mutation.Figure 4
The order of taking measurements and applying chemical mutagens and radiation in a cell cycle.


### 2.2.1. Base-to-Base Deterministic Mutations

Denote the targeted DNA segment withn nucleotide bases at kth step by a column vector xk, as shown in Figure 5. xki is the ith element of xk. Let P be the transfer matrix from xk to xk+1, for all k, k∈ℤ+∪{0}, without mutation. Then, the perfect DNA replication process can be expressed as(1)xk+1=Pxk.Figure 5
Biological information flow in central dogma of molecular biology.Proposition 1.
P=-I.Proof.
As no mutation occurs,xk+1 is completely complementary to xk by Watson-Crick base pairing rule, and xk+2 is completely complementary to xk+1. Therefore, xk+2 is exactly the same as xk. thus,
(2)xk+2=Pxk+1=P2xk⟹P2=I.
Since every base mutates independently, every element of xk+1i only depends on the corresponding element of xki, thus P is diagonal. In addition, xk+1≠xk, we conclude P=-I.Based on Proposition1, we assign values to nucleotide bases set {A,G,C,T,O}, where O is an artificial nonsense base. Define an equivalence relationship between {A,G,C,T,O} and {1,2,-2,-1,0}, that is, {A,G,C,T,O}⇔{1,2,-2,-1,0}, with(3)xki={1,ifA,2,ifG,-2,ifC,-1,ifT,0,ifO.Proposition 2.
{1,2,-2,-1,0} is a field under proper definitions of addition and multiplication.Proof.
Defining the addition table and multiplication table as in Tables1 and 2, we check if the set {1,2,-2,-1,0} satisfies the definition of field.Table 1
Addition table for{1,2,-2,-1,0}.
+12−2−1012−2−1012−2−1012−2−1012−2−1012−2−1012−2−10Table 2
Multiplication table for{1,2,-2,-1,0}.
×12−2−10112−2−1022−11−20−2−21−120−1−1−2210000000Closed under Addition and Multiplication
Satisfied obviously from Tables1 and 2.Associativity of Addition and Multiplication
Implicitly satisfied by integer addition and multiplication.Commutativity of Addition and Multiplication
Satisfied as Tables1 and 2 are symmetric according to the diagonal.Additive and Multiplicative Identity
Additive identity is 0, and multiplicative identity is 1.Additive and Multiplicative Inverses
Additive inverses pair:1↔-1, 2↔-2, 0↔0.
Multiplicative inverses pair:1↔1, 2↔-2, -1↔-1.Distributivity of Multiplication over Addition
Implicitly satisfied by integer addition and multiplication.We conclude{0,1,2,-2,-1} is a field under addition and multiplication defined by Tables 1 and 2.From now on, we useℱ to denote the field {0,1,2,-2,-1}. And xk∈ℱn is the state vector representing a DNA segment with n nucleotide bases, where ℱn is the set of ℱ-valued vectors of dimension n.We start with the simplest form of mutations, point mutation. Suppose there is a point mutation, we write it mathematically as(4)xk+1=(-I+Δs)xk+Δw,
where xk+1,xk∈ℱ, and -I reduces to −1 as only one base is involved. The corresponding values of Δs and Δw, obtained by reverse engineering with all possible pairs of xk and xk+1, are listed in Table 3.Table 3
Possible values ofΔs and Δw. The corresponding values of Δs and Δw are obtained by substituting the value of xk and xk+1 into (4), with Δw≠0 only if xk=0.
kth(k+1)thAGCTOA2−2−101111}ΔsG−120−2C−202−1T0−1−22O12−2−10}ΔwHere,Δs represents the mutation from four normal nucleotide bases, and Δw corresponds to mutation from nonsense base, that is, Δw≠0 only if xk=0.Rewriting (4) by collecting all values of Δs and Δw in Table 3, we get(5a)xk+1=(-I+∑j=04ukjsj)xk+∑j=04ckjwj(5b)=(-I+uks)xk+ckw,where {s0,s1,s2,s3,s4}={w0,w1,w2,w3,w4}={0,1,2,-2,-1}, ukj,ckj∈{0,1}, representing the on/off controls, uk=[uk0uk1uk2uk3uk4], ck=[ck0ck1ck2ck3ck4], and s=w=[012-2-1]T. In (5a), sj and wj are constants for all k and j. ukj and ckj, the inputs of the system, are the on/off controls for chemical mutagens or radiation. Clearly, ∑j=04ckj=1 only if xk=0. Equation (5b) is a simplified version of (5a) as we put ukj,sj,ckj,wj into vector form uk,s,ck,w. s and w serve as vector basis for base-to-base deterministic model. uk and ck are now multi-input controls; each of them contains 5 on/off controls, corresponding to all possible transfer patterns. For a particular k, at most one of ukjs and ckjs can be 1, as stated in Proposition 3. This is consistent with the fact that every state can be transferred to only one of the five states in the state space ℱ with corresponding mutagens available.Proposition 3.
It is always1-1 transfer when mutation occurs, that is, one nucleotide base can only transfer to another one, therefore(i)
ifxk=0 and ck=0, or ck=[10000], then xk+1=0,(ii)
ifxk≠0, then ck=0 and uk is either 0 or a unit row vector,(iii)
ifxk=0, then uk=0 and ck is either 0 or a unit row vector,(iv)
uk+ck is either 0 or a unit row vector, for all k∈ℤ+∪{0}.Now, suppose for some reason we need to take an addition, measurement in the middle of the cell cycle, after the completion of thekth duplication and before the start of the (k+1)th. We name this kind of measurement an intermediate state, and denote by it xk′. Then, we have(6)xk+1=(I+Δs′)xk′+Δw′,
where the values of Δs′ and Δw′, listed in Table 4, are obtained in the same way as getting Δs and Δw in Table 3.Table 4
Possible values ofΔs′ and Δw′. The corresponding values of Δs′ and Δw′ are obtained by substituting the value of xk′ and xk+1 into (4), with Δw′≠0 only if xk′=0.
kth(k+1)thAGCTOA012−2-1-1-1-1}Δs′G20−21C1−202T−2210O12−2−10}Δw′Comparing Tables3 and 4, we find the collection of Δs and Δs′, Δw and Δw′, form the same set, respectively. Thus, we continue using s and w when rewriting (6) in the form of (5a) and (5b), that is,
(7)xk+1=(I+vks)xk′+ck′w,
where vk,ck′ are the counterparts of uk,ck, respectively, and s,w are the same as in (5b).Similar to Proposition3, we get Proposition 4.Proposition 4.
vk and ck′ in (7) need to satisfy the following conditions.(i)
Ifxk=0 and ck′=0, or ck′=[10000], then xk+1=0.(ii)
Ifxk≠0, then ck′=0 and vk is either 0 or a unit row vector.(iii)
Ifxk=0, then vk=0 and ck is either 0 or a unit row vector.(iv)
vk+ck′ is either 0 or a unit row vector, for all k∈ℤ+∪{0}.Now take, both chemical mutagens and radiative rays under our consideration and apply them in the order as shown in Figure4. Then, we can express our system equation as(8a)xk′=(-I+uks︸mutationscausedbychemicalmutagensfromnormalbases)xk+ckw︸mutationscausedbychemicalmutagensfromO,(8b)xk+1=(I+vks︸mutationscausedbyradiativeraysfromnormalbases)xk′+ck′w︸mutationscausedbyradiativeraysfromO,(8c)yk=xk,where uk and vk are the inputs of the system and yk is the measurement. Obviously, (8a) is modified from (5b) and (8b) from (7).The two-step mutation and the intermediate statexk′avoid the case xk is changed to different bases by radiation and chemical mutagens simultaneously, which causes confusion. Substituting (8a) and (8b), we get(9a)xk+1=(I+vks)(-I+uks)xk+(I+vks)ckw+ck′w,(9b)yk=xk.Obviously, Proposition 3 still holds for uk and ck, and Proposition 4 holds for vk and ck′ for (9a).For point mutations, we have 20 on/off controls in total for every stepk, 10 for chemical mutagens as described before, and the rest for radiation.


### 2.2.2. Gene-to-Gene Deterministic Mutations

In general, mutations involve multiple bases. Therefore, large-scale deterministic model is necessary. Now, we show how to extend our model to large-scale systems.Suppose we have a coding DNA segment with lengthn, then xk∈ℱn. Since a coding DNA segment usually contains integer number of codons, which is made of three consecutive bases, n is a multiple of 3. Let xki denote the ith component of xk. This notation is consistent with the one in Section 2.2.1. Initiated by the base-to-base deterministic model from Section 2.2.1, we write our system equation for large-scale system as
(10)xk′=(-I+∑i=1nukiSki︸mutationscausedbychemicalmutagensfromnormalbases)xk+∑i∈OkckiWki︸mutationscausedbychemicalmutagensfromO,xk+1=(I+∑i=1nvkiQki︸mutationscausedbyradiativeraysfromnormalbases)xk′+∑i∈Ok′bkiRki︸mutationscausedbyradiativeraysfromO,yk=xk,
where uki,vki,cki,bki are on/off controls of the ith element,Ski,Qki are n×n square matrices corresponding to the mutations between normal bases or from normal bases by chemicals and radiation, respectively, Wki,Rki are n-dimensional column vectors representing mutations from nonsense bases by chemicals and radiation, respectively, and 𝒪k={i:xki=0,1≤i≤n}, 𝒪k′={i:xk′i=0,1≤i≤n}.Ski and Qki are diagonal matrices since each base mutates independently. The values in the first four rows of Tables 3 and 4 correspond to the diagonal elements of Ski and Qki, respectively. The last rows of Tables 3 and 4 are assigned to Wki and Rki, n-dimensional vectors, at nonsense base’s spots for xk.Define𝒮={sjeieiT,∀i,j,0≤j≤4,1≤i≤n}, a collection of n×n matrices, where sj is the same as in (5a) and (5b), ei is the unit column vector of length n with ith component equal to 1 and all other components equal to 0, and eieiT is the square matrix with only the ith element on the diagonal equals to 1, and 0 otherwise. Then, Ski,Qki can be written as linear combinations of all elements from 𝒮, with the coefficient of each element either 0 or 1 corresponding to the on/off control uk(i,j) and vk(i,j), respectively.Similarly, define𝒲={wjei,∀i,j,0≤j≤4,1≤i≤n}, where wj is the same as (5a) and (5b). Wki,Rki can be written as linear combinations of all components from 𝒲, with coefficient of every component either 0 or 1 corresponding to the on/off control ck(i,j) and ck′(i,j), respectively.Therefore, instead of using step-varyingSki,Ski,Wki,Rki, we find matrix basis for those four square matrices to make the controls to be the only variables depending on k, as we did for single-base cases. Then, we can, write (10) as(11a)xk′=(-I+∑i=1n∑j=04uk(i,j)sjeieiT)xk+∑i∈Ok∑j=04ck(i,j)wjei,(11b)xk+1=(I+∑i=1n∑j=04vk(i,j)sjeieiT)xk′+∑i∈Ok′∑j=04ck′(i,j)wjei,(11c)yk=xk,where uk(i,j),vk(i,j),ck(i,j),ck′(i,j)∈{0,1}.As shown in (11a), (11b), and (11c), multisites mutations contain 20n controls in total for every step k, where n is the number of nucleotide bases on the targeted gene. Similar to point mutations, every single site has 20 controls in each step, 10 for chemical mutagens and 10 for radiation.We can viewuk,vk,ck,ck′ as binary matrices of dimension n×5, and uk(i,j), vk(i,j)ck(i,j), ck′(i,j)are the corresponding element of ith row and jth column. Use uki,vki,cki,bki, binary row vectors of dimension 5, to denote the ith row of uk,vk,ck,ck′, respectively. Again, s=w=[012-2-1]T.Combining (11a) and (11b), and writing control variables in vector forms, we get(12)xk+1=(I+∑i=1nvkiseieiT)(-I+∑i=1nukiseieiT)xk+(I+∑i=1nvkiseieiT)∑i∈Okckiwei+∑i∈Ok′bkiwei,yk=xk.Proposition 5.
For large-scale deterministic system,uk,vk,ck,ck′ satisfy conditions below.(i)
IfeiTxk=0, then i∈𝒪k.(ii)
IfeiTxk=0, cki=0 or cki=[10000], then i∈𝒪k′.(iii)
For alli∉𝒪k, uki is either 0 or a row unit vector and cki=0.(iv)
For alli∈𝒪k, cki is either 0 a row unit vector and uki=0.(v)
For alli∉𝒪k′, vki is either 0 or a row unit vector and bki=0.(vi)
For alli∈𝒪k′, bki is either 0 or a row unit vector and vki=0.(vii)
For alli,k, 1≤i≤n,  k∈ℤ+∪{0}, uki+cki is either 0 or a unit row vector and vki+bki is either 0 or a unit row vector.The mathematical model (12) is quite flexible and can be easily extended to many cases, such as transcription process, multiple spot mutations within one-step or broken DNA strands.Take broken DNA strands as an example. DNA strand breaks due to various reasons. Our system equation can represent this phenomenon by dividing the whole system into small subsystems. Significant brokage of DNA strands is simply eliminated by cell mechanism to ensure the accuracy to DNA replication. Equation (13) shows the case of one single DNA strand breaking into two segments by chemical mutagens (13)(xk′(1)xk′(2))=(-Im+∑i=1mukiseieiT00-In-m+∑i=m+1nukiseieiT)×(xk(1)xk(2))+(∑i∈Ok,1≤i≤mckiwei∑i∈Ok,(m+1)≤i≤nckiwei),xk+1(1)=(Im+∑i=1mvkiseieiT)xk′(1)+∑i∈Ok′,1≤i≤mbkiwei,xk+1(2)=(In-m+∑i=m+1nvkiseieiT)xk′(2)+∑i∈Ok′,(m+1)≤i≤nbkiwei.


### 2.2.3. Gene-to-Gene Stochastic Mutations

In reality, mutagens, no matter chemical or radiative, always cause randomness in mutation. Therefore, we need to derive the model for gene-to-gene stochastic mutations.Introduce new random variables,hk,l1(i,j),rk,l2(i,j),hk,l3′(i,j),rk,l4′(i,j)∈{0,1}, associated with probability pl1,j(h), pl2,j(r),pl3,j(h′),pl4,j(r), for all i,k, 1≤i≤n,  k∈ℤ+∪{0}, respectively, where k is the step index, l1,l2 are indices for chemical mutagens inducing mutation from normal bases and from O, respectively, l3,l4 are indices for radiation inducing mutation from normal bases and from O, respectively, i is the index of DNA segment, and the value of j corresponds to the transfer pattern, which can be found in Tables 3 and 4. Note different mutagens have different probability assignments, the probability assignments are only related to the type of mutagens, and the probability associated with every kind of mutagens sums up to 1, that is,
(14)∑j=04pl1,j(h)=1,∀l1,1≤l1≤l,∑j=04pl2,j(r)=1,∀l2,1≤l2≤m,∑j=04pl3,j(h′)=1,∀l3,1≤l3≤l′,∑j=04pl4,j(r′)=1,∀l4,1≤l4≤m′.
The controls are uk,l1i,ck,l2i,vk,l3i,bk,l4i∈{0,1}, with the fact that 1 representing mutagen with corresponding index is applied at ith spot of DNA segment at kth generation, and 0 representing mutagen with corresponding index is not applied at spot i at kth step, similar to Sections 2.2.1 and 2.2.2. The mutagen indices l1,l2,l3,l4 can be omitted in deterministic mutations since given the current state and control, the next state is unique. However, they are necessary for stochastic mutations, because there exist multiple possible states for the next stage given the control. In other words, the next state is determined by random variables hk,l1(i,j),rk,l2(i,j),hk,l3′(i,j),rk,l4′(i,j), given the values of uk,l1i,ck,l2i,vk,l3i,bk,l4i, and xk.Suppose we have(l+m) kinds of chemical mutagens available, with l kinds to induce mutations from normal bases and m kinds to induce mutations from O. And we have (l′+m′) kinds of radiation available, with l′ kinds to induce mutations from normal bases and m′ kinds to induce mutations from O. Therefore, we have total (l+m+l′+m′) controls for each spot i at each step k. We can write our system equation as(15a)xk′=(-I+∑l1=1l∑i=1nuk,l1i∑j=04hk,l1(i,j)sjeieiT︸mutationscausedbychemicalmutagensfromnormalbases)xk+∑l2=1m∑i∈Okck,l2i∑j=04rk,l2(i,j)wjei︸mutationscausedbychemicalmutagensfromO,(15b)xk+1=(I+∑l3=1l′∑i=1nvk,l3i∑j=04hk,l3′(i,j)sjeieiT︸mutationscausedbyradiativeraysfromnormalbases)xk′+∑l4=1m′∑i∈Ok′bk,l4i∑j=04rk,l4′(i,j)wjei︸mutationscausedbyradiativeraysfromO,(15c)yk=xk.Again, we define hk,l1(i,j),rk,l2(i,j),hk,l3′(i,j),rk,l4′(i,j) the elements at ith row and jth column of n×5 binary matrices hk,l1,rk,l2,hk,l3′,rk,l4′, respectively. hk,l1i,rk,l2i,hk,l3′i, rk,l4′i, and binary row vectors of dimension 5 denote the ith row of hk,l1,rk,l2,hk,l3′,rk,l4′, respectively. Then, we can simplify (15a), (15b), and (15c) and combine (15a) and (15b) as(16)xk+1=(I+∑l3=1l′∑i=1nvk,l3ihk,l3′iseieiT)×(-I+∑l1=1l∑i=1nuk,l1ihk,l1iseieiT)xk+(I+∑l3=1l′∑i=1nvk,l3ihk,l3′iseieiT)∑l2=1m∑i∈Okck,il2rk,l2iwei+∑l4=1m′∑i∈Ok′bk,l4irk,l4′iwei,yk=xk.Proposition 6.
For large-scale stochastic system,uk,l1i,hk,l1i,ck,l2i,rk,l2i,vk,l3i,hk,l3′i,bk,l4i,rk,l4′i follow the rules below.(i)
IfeiTxk=0, then i∈𝒪k.(ii)
IfeiTxk=0 and ∑l2=1mck,l2i=0, then i∈𝒪k′.(iii)
IfeiTxk=0, ∑l2=1mck,l2i=1 and rk,l2i=[10000], then i∈𝒪k′.(iv)
For alli,k,l1, 1≤i≤n,  k∈ℤ+∪{0},1≤l1≤l, if uk,l1i=1, then hk,l1i is a unit row vector.(v)
For alli,k,l2, 1≤i≤n,  k∈ℤ+∪{0}, 1≤l2≤m, if ck,l2i=1, then rk,l2i is a unit row vector.(vi)
For alli,k,l3, 1≤i≤n,  k∈ℤ+∪{0}, 1≤l3≤l′, if vk,l3i=1, then hk,l3′i is a unit row vector.(vii)
For alli,k,l4, 1≤i≤n,  k∈ℤ+∪{0}, 1≤l4≤m′, if bk,l4i=1, then rk,l4′i is a unit row vector.(viii)
For alli∉𝒪k, ∑l1=1luk,l1i=0or1 and ck,l2i=0, for all l2,1≤l2≤m.(ix)
For alli∈𝒪k, ∑l2=1mck,l2i=0 or 1 and uk,l1i=0, for all l1, 1≤l1≤l.(x)
For alli∉𝒪k′, ∑l3=1l′vk,l3i=0 or 1 and bk,l4i=0, for alll4, 1≤l4≤m′.(xi)
For alli∈𝒪k′, ∑l4=1m′bk,l4i=0 or 1 and vk,l3i=0, for all l3, 1≤l3≤l′′.(xii)
For alli,k, 1≤i≤n,  k∈ℤ+∪{0}, ∑l1=1luk,l1i+∑l2=1mck,l2i=0 or 1 and ∑l3=1l′vk,l3i+∑l4=1m′bk,l4i=0 or 1.We close this section with the definition of controllability to the system equations proposed above. DNA replication systems with system equations proposed as (9a), (9b), (12), and (16) are completely controllable if and only if  for all x0,x2k1,x2k2+1∈ℱ, k1,k2∈ℤ+∪{0},∃ at least one path from x0 to x2k1 and at least one path from x0 to x2k2+1 by applying proper mutagens in the correct order, with k1,k2 finite.


## 2.2.1. Base-to-Base Deterministic Mutations

Denote the targeted DNA segment withn nucleotide bases at kth step by a column vector xk, as shown in Figure 5. xki is the ith element of xk. Let P be the transfer matrix from xk to xk+1, for all k, k∈ℤ+∪{0}, without mutation. Then, the perfect DNA replication process can be expressed as(1)xk+1=Pxk.Figure 5
Biological information flow in central dogma of molecular biology.Proposition 1.
P=-I.Proof.
As no mutation occurs,xk+1 is completely complementary to xk by Watson-Crick base pairing rule, and xk+2 is completely complementary to xk+1. Therefore, xk+2 is exactly the same as xk. thus,
(2)xk+2=Pxk+1=P2xk⟹P2=I.
Since every base mutates independently, every element of xk+1i only depends on the corresponding element of xki, thus P is diagonal. In addition, xk+1≠xk, we conclude P=-I.Based on Proposition1, we assign values to nucleotide bases set {A,G,C,T,O}, where O is an artificial nonsense base. Define an equivalence relationship between {A,G,C,T,O} and {1,2,-2,-1,0}, that is, {A,G,C,T,O}⇔{1,2,-2,-1,0}, with(3)xki={1,ifA,2,ifG,-2,ifC,-1,ifT,0,ifO.Proposition 2.
{1,2,-2,-1,0} is a field under proper definitions of addition and multiplication.Proof.
Defining the addition table and multiplication table as in Tables1 and 2, we check if the set {1,2,-2,-1,0} satisfies the definition of field.Table 1
Addition table for{1,2,-2,-1,0}.
+12−2−1012−2−1012−2−1012−2−1012−2−1012−2−1012−2−10Table 2
Multiplication table for{1,2,-2,-1,0}.
×12−2−10112−2−1022−11−20−2−21−120−1−1−2210000000Closed under Addition and Multiplication
Satisfied obviously from Tables1 and 2.Associativity of Addition and Multiplication
Implicitly satisfied by integer addition and multiplication.Commutativity of Addition and Multiplication
Satisfied as Tables1 and 2 are symmetric according to the diagonal.Additive and Multiplicative Identity
Additive identity is 0, and multiplicative identity is 1.Additive and Multiplicative Inverses
Additive inverses pair:1↔-1, 2↔-2, 0↔0.
Multiplicative inverses pair:1↔1, 2↔-2, -1↔-1.Distributivity of Multiplication over Addition
Implicitly satisfied by integer addition and multiplication.We conclude{0,1,2,-2,-1} is a field under addition and multiplication defined by Tables 1 and 2.From now on, we useℱ to denote the field {0,1,2,-2,-1}. And xk∈ℱn is the state vector representing a DNA segment with n nucleotide bases, where ℱn is the set of ℱ-valued vectors of dimension n.We start with the simplest form of mutations, point mutation. Suppose there is a point mutation, we write it mathematically as(4)xk+1=(-I+Δs)xk+Δw,
where xk+1,xk∈ℱ, and -I reduces to −1 as only one base is involved. The corresponding values of Δs and Δw, obtained by reverse engineering with all possible pairs of xk and xk+1, are listed in Table 3.Table 3
Possible values ofΔs and Δw. The corresponding values of Δs and Δw are obtained by substituting the value of xk and xk+1 into (4), with Δw≠0 only if xk=0.
kth(k+1)thAGCTOA2−2−101111}ΔsG−120−2C−202−1T0−1−22O12−2−10}ΔwHere,Δs represents the mutation from four normal nucleotide bases, and Δw corresponds to mutation from nonsense base, that is, Δw≠0 only if xk=0.Rewriting (4) by collecting all values of Δs and Δw in Table 3, we get(5a)xk+1=(-I+∑j=04ukjsj)xk+∑j=04ckjwj(5b)=(-I+uks)xk+ckw,where {s0,s1,s2,s3,s4}={w0,w1,w2,w3,w4}={0,1,2,-2,-1}, ukj,ckj∈{0,1}, representing the on/off controls, uk=[uk0uk1uk2uk3uk4], ck=[ck0ck1ck2ck3ck4], and s=w=[012-2-1]T. In (5a), sj and wj are constants for all k and j. ukj and ckj, the inputs of the system, are the on/off controls for chemical mutagens or radiation. Clearly, ∑j=04ckj=1 only if xk=0. Equation (5b) is a simplified version of (5a) as we put ukj,sj,ckj,wj into vector form uk,s,ck,w. s and w serve as vector basis for base-to-base deterministic model. uk and ck are now multi-input controls; each of them contains 5 on/off controls, corresponding to all possible transfer patterns. For a particular k, at most one of ukjs and ckjs can be 1, as stated in Proposition 3. This is consistent with the fact that every state can be transferred to only one of the five states in the state space ℱ with corresponding mutagens available.Proposition 3.
It is always1-1 transfer when mutation occurs, that is, one nucleotide base can only transfer to another one, therefore(i)
ifxk=0 and ck=0, or ck=[10000], then xk+1=0,(ii)
ifxk≠0, then ck=0 and uk is either 0 or a unit row vector,(iii)
ifxk=0, then uk=0 and ck is either 0 or a unit row vector,(iv)
uk+ck is either 0 or a unit row vector, for all k∈ℤ+∪{0}.Now, suppose for some reason we need to take an addition, measurement in the middle of the cell cycle, after the completion of thekth duplication and before the start of the (k+1)th. We name this kind of measurement an intermediate state, and denote by it xk′. Then, we have(6)xk+1=(I+Δs′)xk′+Δw′,
where the values of Δs′ and Δw′, listed in Table 4, are obtained in the same way as getting Δs and Δw in Table 3.Table 4
Possible values ofΔs′ and Δw′. The corresponding values of Δs′ and Δw′ are obtained by substituting the value of xk′ and xk+1 into (4), with Δw′≠0 only if xk′=0.
kth(k+1)thAGCTOA012−2-1-1-1-1}Δs′G20−21C1−202T−2210O12−2−10}Δw′Comparing Tables3 and 4, we find the collection of Δs and Δs′, Δw and Δw′, form the same set, respectively. Thus, we continue using s and w when rewriting (6) in the form of (5a) and (5b), that is,
(7)xk+1=(I+vks)xk′+ck′w,
where vk,ck′ are the counterparts of uk,ck, respectively, and s,w are the same as in (5b).Similar to Proposition3, we get Proposition 4.Proposition 4.
vk and ck′ in (7) need to satisfy the following conditions.(i)
Ifxk=0 and ck′=0, or ck′=[10000], then xk+1=0.(ii)
Ifxk≠0, then ck′=0 and vk is either 0 or a unit row vector.(iii)
Ifxk=0, then vk=0 and ck is either 0 or a unit row vector.(iv)
vk+ck′ is either 0 or a unit row vector, for all k∈ℤ+∪{0}.Now take, both chemical mutagens and radiative rays under our consideration and apply them in the order as shown in Figure4. Then, we can express our system equation as(8a)xk′=(-I+uks︸mutationscausedbychemicalmutagensfromnormalbases)xk+ckw︸mutationscausedbychemicalmutagensfromO,(8b)xk+1=(I+vks︸mutationscausedbyradiativeraysfromnormalbases)xk′+ck′w︸mutationscausedbyradiativeraysfromO,(8c)yk=xk,where uk and vk are the inputs of the system and yk is the measurement. Obviously, (8a) is modified from (5b) and (8b) from (7).The two-step mutation and the intermediate statexk′avoid the case xk is changed to different bases by radiation and chemical mutagens simultaneously, which causes confusion. Substituting (8a) and (8b), we get(9a)xk+1=(I+vks)(-I+uks)xk+(I+vks)ckw+ck′w,(9b)yk=xk.Obviously, Proposition 3 still holds for uk and ck, and Proposition 4 holds for vk and ck′ for (9a).For point mutations, we have 20 on/off controls in total for every stepk, 10 for chemical mutagens as described before, and the rest for radiation.


## 2.2.2. Gene-to-Gene Deterministic Mutations

In general, mutations involve multiple bases. Therefore, large-scale deterministic model is necessary. Now, we show how to extend our model to large-scale systems.Suppose we have a coding DNA segment with lengthn, then xk∈ℱn. Since a coding DNA segment usually contains integer number of codons, which is made of three consecutive bases, n is a multiple of 3. Let xki denote the ith component of xk. This notation is consistent with the one in Section 2.2.1. Initiated by the base-to-base deterministic model from Section 2.2.1, we write our system equation for large-scale system as
(10)xk′=(-I+∑i=1nukiSki︸mutationscausedbychemicalmutagensfromnormalbases)xk+∑i∈OkckiWki︸mutationscausedbychemicalmutagensfromO,xk+1=(I+∑i=1nvkiQki︸mutationscausedbyradiativeraysfromnormalbases)xk′+∑i∈Ok′bkiRki︸mutationscausedbyradiativeraysfromO,yk=xk,
where uki,vki,cki,bki are on/off controls of the ith element,Ski,Qki are n×n square matrices corresponding to the mutations between normal bases or from normal bases by chemicals and radiation, respectively, Wki,Rki are n-dimensional column vectors representing mutations from nonsense bases by chemicals and radiation, respectively, and 𝒪k={i:xki=0,1≤i≤n}, 𝒪k′={i:xk′i=0,1≤i≤n}.Ski and Qki are diagonal matrices since each base mutates independently. The values in the first four rows of Tables 3 and 4 correspond to the diagonal elements of Ski and Qki, respectively. The last rows of Tables 3 and 4 are assigned to Wki and Rki, n-dimensional vectors, at nonsense base’s spots for xk.Define𝒮={sjeieiT,∀i,j,0≤j≤4,1≤i≤n}, a collection of n×n matrices, where sj is the same as in (5a) and (5b), ei is the unit column vector of length n with ith component equal to 1 and all other components equal to 0, and eieiT is the square matrix with only the ith element on the diagonal equals to 1, and 0 otherwise. Then, Ski,Qki can be written as linear combinations of all elements from 𝒮, with the coefficient of each element either 0 or 1 corresponding to the on/off control uk(i,j) and vk(i,j), respectively.Similarly, define𝒲={wjei,∀i,j,0≤j≤4,1≤i≤n}, where wj is the same as (5a) and (5b). Wki,Rki can be written as linear combinations of all components from 𝒲, with coefficient of every component either 0 or 1 corresponding to the on/off control ck(i,j) and ck′(i,j), respectively.Therefore, instead of using step-varyingSki,Ski,Wki,Rki, we find matrix basis for those four square matrices to make the controls to be the only variables depending on k, as we did for single-base cases. Then, we can, write (10) as(11a)xk′=(-I+∑i=1n∑j=04uk(i,j)sjeieiT)xk+∑i∈Ok∑j=04ck(i,j)wjei,(11b)xk+1=(I+∑i=1n∑j=04vk(i,j)sjeieiT)xk′+∑i∈Ok′∑j=04ck′(i,j)wjei,(11c)yk=xk,where uk(i,j),vk(i,j),ck(i,j),ck′(i,j)∈{0,1}.As shown in (11a), (11b), and (11c), multisites mutations contain 20n controls in total for every step k, where n is the number of nucleotide bases on the targeted gene. Similar to point mutations, every single site has 20 controls in each step, 10 for chemical mutagens and 10 for radiation.We can viewuk,vk,ck,ck′ as binary matrices of dimension n×5, and uk(i,j), vk(i,j)ck(i,j), ck′(i,j)are the corresponding element of ith row and jth column. Use uki,vki,cki,bki, binary row vectors of dimension 5, to denote the ith row of uk,vk,ck,ck′, respectively. Again, s=w=[012-2-1]T.Combining (11a) and (11b), and writing control variables in vector forms, we get(12)xk+1=(I+∑i=1nvkiseieiT)(-I+∑i=1nukiseieiT)xk+(I+∑i=1nvkiseieiT)∑i∈Okckiwei+∑i∈Ok′bkiwei,yk=xk.Proposition 5.
For large-scale deterministic system,uk,vk,ck,ck′ satisfy conditions below.(i)
IfeiTxk=0, then i∈𝒪k.(ii)
IfeiTxk=0, cki=0 or cki=[10000], then i∈𝒪k′.(iii)
For alli∉𝒪k, uki is either 0 or a row unit vector and cki=0.(iv)
For alli∈𝒪k, cki is either 0 a row unit vector and uki=0.(v)
For alli∉𝒪k′, vki is either 0 or a row unit vector and bki=0.(vi)
For alli∈𝒪k′, bki is either 0 or a row unit vector and vki=0.(vii)
For alli,k, 1≤i≤n,  k∈ℤ+∪{0}, uki+cki is either 0 or a unit row vector and vki+bki is either 0 or a unit row vector.The mathematical model (12) is quite flexible and can be easily extended to many cases, such as transcription process, multiple spot mutations within one-step or broken DNA strands.Take broken DNA strands as an example. DNA strand breaks due to various reasons. Our system equation can represent this phenomenon by dividing the whole system into small subsystems. Significant brokage of DNA strands is simply eliminated by cell mechanism to ensure the accuracy to DNA replication. Equation (13) shows the case of one single DNA strand breaking into two segments by chemical mutagens (13)(xk′(1)xk′(2))=(-Im+∑i=1mukiseieiT00-In-m+∑i=m+1nukiseieiT)×(xk(1)xk(2))+(∑i∈Ok,1≤i≤mckiwei∑i∈Ok,(m+1)≤i≤nckiwei),xk+1(1)=(Im+∑i=1mvkiseieiT)xk′(1)+∑i∈Ok′,1≤i≤mbkiwei,xk+1(2)=(In-m+∑i=m+1nvkiseieiT)xk′(2)+∑i∈Ok′,(m+1)≤i≤nbkiwei.


## 2.2.3. Gene-to-Gene Stochastic Mutations

In reality, mutagens, no matter chemical or radiative, always cause randomness in mutation. Therefore, we need to derive the model for gene-to-gene stochastic mutations.Introduce new random variables,hk,l1(i,j),rk,l2(i,j),hk,l3′(i,j),rk,l4′(i,j)∈{0,1}, associated with probability pl1,j(h), pl2,j(r),pl3,j(h′),pl4,j(r), for all i,k, 1≤i≤n,  k∈ℤ+∪{0}, respectively, where k is the step index, l1,l2 are indices for chemical mutagens inducing mutation from normal bases and from O, respectively, l3,l4 are indices for radiation inducing mutation from normal bases and from O, respectively, i is the index of DNA segment, and the value of j corresponds to the transfer pattern, which can be found in Tables 3 and 4. Note different mutagens have different probability assignments, the probability assignments are only related to the type of mutagens, and the probability associated with every kind of mutagens sums up to 1, that is,
(14)∑j=04pl1,j(h)=1,∀l1,1≤l1≤l,∑j=04pl2,j(r)=1,∀l2,1≤l2≤m,∑j=04pl3,j(h′)=1,∀l3,1≤l3≤l′,∑j=04pl4,j(r′)=1,∀l4,1≤l4≤m′.
The controls are uk,l1i,ck,l2i,vk,l3i,bk,l4i∈{0,1}, with the fact that 1 representing mutagen with corresponding index is applied at ith spot of DNA segment at kth generation, and 0 representing mutagen with corresponding index is not applied at spot i at kth step, similar to Sections 2.2.1 and 2.2.2. The mutagen indices l1,l2,l3,l4 can be omitted in deterministic mutations since given the current state and control, the next state is unique. However, they are necessary for stochastic mutations, because there exist multiple possible states for the next stage given the control. In other words, the next state is determined by random variables hk,l1(i,j),rk,l2(i,j),hk,l3′(i,j),rk,l4′(i,j), given the values of uk,l1i,ck,l2i,vk,l3i,bk,l4i, and xk.Suppose we have(l+m) kinds of chemical mutagens available, with l kinds to induce mutations from normal bases and m kinds to induce mutations from O. And we have (l′+m′) kinds of radiation available, with l′ kinds to induce mutations from normal bases and m′ kinds to induce mutations from O. Therefore, we have total (l+m+l′+m′) controls for each spot i at each step k. We can write our system equation as(15a)xk′=(-I+∑l1=1l∑i=1nuk,l1i∑j=04hk,l1(i,j)sjeieiT︸mutationscausedbychemicalmutagensfromnormalbases)xk+∑l2=1m∑i∈Okck,l2i∑j=04rk,l2(i,j)wjei︸mutationscausedbychemicalmutagensfromO,(15b)xk+1=(I+∑l3=1l′∑i=1nvk,l3i∑j=04hk,l3′(i,j)sjeieiT︸mutationscausedbyradiativeraysfromnormalbases)xk′+∑l4=1m′∑i∈Ok′bk,l4i∑j=04rk,l4′(i,j)wjei︸mutationscausedbyradiativeraysfromO,(15c)yk=xk.Again, we define hk,l1(i,j),rk,l2(i,j),hk,l3′(i,j),rk,l4′(i,j) the elements at ith row and jth column of n×5 binary matrices hk,l1,rk,l2,hk,l3′,rk,l4′, respectively. hk,l1i,rk,l2i,hk,l3′i, rk,l4′i, and binary row vectors of dimension 5 denote the ith row of hk,l1,rk,l2,hk,l3′,rk,l4′, respectively. Then, we can simplify (15a), (15b), and (15c) and combine (15a) and (15b) as(16)xk+1=(I+∑l3=1l′∑i=1nvk,l3ihk,l3′iseieiT)×(-I+∑l1=1l∑i=1nuk,l1ihk,l1iseieiT)xk+(I+∑l3=1l′∑i=1nvk,l3ihk,l3′iseieiT)∑l2=1m∑i∈Okck,il2rk,l2iwei+∑l4=1m′∑i∈Ok′bk,l4irk,l4′iwei,yk=xk.Proposition 6.
For large-scale stochastic system,uk,l1i,hk,l1i,ck,l2i,rk,l2i,vk,l3i,hk,l3′i,bk,l4i,rk,l4′i follow the rules below.(i)
IfeiTxk=0, then i∈𝒪k.(ii)
IfeiTxk=0 and ∑l2=1mck,l2i=0, then i∈𝒪k′.(iii)
IfeiTxk=0, ∑l2=1mck,l2i=1 and rk,l2i=[10000], then i∈𝒪k′.(iv)
For alli,k,l1, 1≤i≤n,  k∈ℤ+∪{0},1≤l1≤l, if uk,l1i=1, then hk,l1i is a unit row vector.(v)
For alli,k,l2, 1≤i≤n,  k∈ℤ+∪{0}, 1≤l2≤m, if ck,l2i=1, then rk,l2i is a unit row vector.(vi)
For alli,k,l3, 1≤i≤n,  k∈ℤ+∪{0}, 1≤l3≤l′, if vk,l3i=1, then hk,l3′i is a unit row vector.(vii)
For alli,k,l4, 1≤i≤n,  k∈ℤ+∪{0}, 1≤l4≤m′, if bk,l4i=1, then rk,l4′i is a unit row vector.(viii)
For alli∉𝒪k, ∑l1=1luk,l1i=0or1 and ck,l2i=0, for all l2,1≤l2≤m.(ix)
For alli∈𝒪k, ∑l2=1mck,l2i=0 or 1 and uk,l1i=0, for all l1, 1≤l1≤l.(x)
For alli∉𝒪k′, ∑l3=1l′vk,l3i=0 or 1 and bk,l4i=0, for alll4, 1≤l4≤m′.(xi)
For alli∈𝒪k′, ∑l4=1m′bk,l4i=0 or 1 and vk,l3i=0, for all l3, 1≤l3≤l′′.(xii)
For alli,k, 1≤i≤n,  k∈ℤ+∪{0}, ∑l1=1luk,l1i+∑l2=1mck,l2i=0 or 1 and ∑l3=1l′vk,l3i+∑l4=1m′bk,l4i=0 or 1.We close this section with the definition of controllability to the system equations proposed above. DNA replication systems with system equations proposed as (9a), (9b), (12), and (16) are completely controllable if and only if  for all x0,x2k1,x2k2+1∈ℱ, k1,k2∈ℤ+∪{0},∃ at least one path from x0 to x2k1 and at least one path from x0 to x2k2+1 by applying proper mutagens in the correct order, with k1,k2 finite.


## 2.3. Generalized Optimal Control Problem Formulation

We first define our objective function that can be adapted to all kinds of systems proposed in Section2.2 with minor changes. Mathematically, in systems where the controllable parameters of interest are discrete, the objective function is usually a weighted sum representing the number of times that a piece of equipment is turned “on” or “off” or the number of resources needed to execute certain tasks in the frequent cases [30]. In our case, this summation is the total number of times that different mutagens are applied weighted by the corresponding cost (including the risk). Another key factor of objective function is the off-trajectory penalty. Designing a trajectory beforehand is necessary to avoid other hidden risks. If the measurement indicates that the current state is off the predefined trajectory, we include a distance reference between current state and desired state as penalty and change the treatment plan accordingly.Therefore, our objective function can be expressed as(17)J0(x0)=min⁡u,c,v,c′Eh,h′,r,r′[∑k=0N-1∑l1=1l∑i=1nαl1uk,l1i+∑k=0N-1∑l2=1m∑i=1nβl2ck,l2i︸costofapplyingchemicalmutagens+∑k=0N-1∑l3=1l′∑i=1nαl3′vk,l3i+∑k=0N-1∑l4=1m′∑i=1nβl4′bk,l4i︸costofapplyingradiativerays+∑k=0Nd(xk,{xkd})︸tracingcost],
with x0,xkd∈ℱn, 1≤k≤N, n≡0 (mod3) given. The physical meaning of uk,l1i, ck,l2i, vk,l3i, bk,l4i, l1, l2, l3, l4 is the same as in Section 2.2.3. αl1,βl2,αl3′,βl4′∈ℝ, for alll1,l2,l3,l4, 1≤l1≤l,  1≤l2≤m, 1≤l3≤l′,  1≤l4≤m′, are the corresponding cost of applying chemical mutagens and radiative rays indexed l1,l2,l3,l4, respectively. {xkd}:ℱn×ℱn→ℝ+∪{0} denotes the desired set at kth stage, generated by the DNA sequences representing the same amino acid sequence asxkd, the desired state at kth step. And d(xk,{xkd}) is the distance reference of the current state xk to the desired set {xkd} at kth step. The final penalty, the distance reference from the final state to the desired set at k=N, is included in the last term.In general,βl2,βl4′≪αl1,αl3′, for all l1,l2,l3,l4,  1≤l1≤l,  1≤l2≤m,  1≤l3≤l′,  1≤l4≤m′, because physically O is a set of nonsense bases and more details are necessary to convert an O back to normal bases, for instance, the cost to identify the exact element in the set O. Our goal is to drive our system optimally from initial state x0 to the desired final set {xNd} by applying a sequence of mutagens indexed with {l1,l2,l3,l4}, at problematic positions i, and in a correct order k.In (17), the first four terms inside the expectation do not depend on random variables hk,l1i, rk,l2i, hk,l3′i, and rk,l4′i, for all i,k,l1,l2,l3,l4 as the treatment plan is computed based on the initial state x0. Given yk, the updated treatment plan is computed accordingly but still not related to random variables. The last term inside expectation, ∑k=0Nd(xk,{xkd}), is the only term in summation that depends on the distribution of the random variables.The constraint of the optimal control problem, in general, is the system equation. We choose multidimensional stochastic system equation as the generalized constraints as it can be degenerated to one-dimensional and multidimensional deterministic cases by proper modifications.Therefore, we can rewrite our objective function and formulate our optimal control problem as(18)J0(x0)=min⁡{u,c,v,c′}0,1,…,N-1[∑k=0N-1∑l1=1l′∑i=1nαl1uk,l1i+∑k=0N-1∑l2=1m∑i=1nβl2ck,l2i+∑k=0N-1∑l3=1l′∑i=1nαl3′vk,l3i+∑k=0N-1∑l4=1m′∑i=1nβl4′bk,l4i+∑k=0NE{h,r,h′,r′}0,1,…,N-1[d(xk,{xkd})]],
subject to(19)xk+1=(I+∑l3=1l′∑i=1nvk,l3ihk,l3′iseieiT)×(-I+∑l1=1l∑i=1nuk,l1ihk,l1iseieiT)xk+(I+∑l3=1l′∑i=1nvk,l3ihk,l3′iseieiT)∑l2=1m∑i∈Okck,l2irk,l2iwei+∑l4=1m′∑i∈Ok′bk,l4irk,l4′iwei,yk=xk.
We need to choose a proper distance reference to quantitatively describe the relationship between DNA segments of same length. We first define the distance reference between codons, and the distance reference between DNA segments is a weighted sum of distance reference between every pair of codons.The distance reference between codons,d(φ1,φ2),φ1,φ2∈ℱ3, needs to fulfill the biological requirements as below.(1)
Nonnegativity: the distance reference between any two codons is either positive or zero. Mathematically,d:ℱ3×ℱ3→ℝ+∪{0}, d(φ1,φ2)≥0.(2)
The distance reference between two codons corresponding to the same amino acid is zero.(3)
Symmetry: the distance reference from codonφ1 to codon φ2 equals the distance reference from codon φ2 to codon φ1, that is, d(φ1,φ2)=d(φ2,φ1).(4)
The distance reference between two codons corresponding to different amino acids should reveal the chemical and physical differences between two amino acids.(5)
The distance reference from stop codons to all other codons is much larger than those between other codons as early termination of amino acid sequences is more harmful than other forms of mutations.All the existing metric defined on the finite field cannot achieve all the requirements above. The second requirement violates the identity of indiscernible, that is,d(φ1,φ2)=0 if and only if φ1=φ2. The redundancy in genetic codes implies d(φ1,φ2)=0 if those two amino acids, φ1 and φ2, are translated into the same amino acids. In addition, the triangular inequality is not necessarily true, according to the underlying physical meanings. We take the assumption that the stop codons are of the same distance reference from and to all other codons.Important physical and chemical properties are listed in Table5. We ignore codons containing O since their chemical and physical properties cannot be found in literature.Table 5
Properties of amino acids.
Amino AcidAbbrev.Codon(s)PolarityPHSizeAlanineAlaGCT,GCC,GCA,GCGNonpolar6.01TinyArginineArgCGA,CGG,CGC,CGT,
AGA,AGGPolar10.76NormalAsparagineAsnAAC,AATPolar5.41SmallAspartic acidAspGAT,GACPolar2.85SmallCysteineCysTGT,TGCNonpolar5.05SmallGlutamineGlnCAA,CAGPolar5.65NormalGlutamic acidGluGAA,GAGPolar3.15NormalGlycineGlyGGA,GGG,GGC,GGTNonpolar6.06TinyHistidineHisCAC,CATPolar7.60NormalIsoleucineIleATA,ATC,ATTNonpolar6.05NormalLeucineLeuTTA,TTG,CTA,CTG,
CTC,CTTNonpolar6.01NormalLysineLysAAA,AAGPolar9.60NormalMethionineMetATGNonpolar5.74NormalPhenylalaninePheTTC,TTTNonpolar5.49NormalProlineProCCA,CCG,CCC,CCTNonpolar6.30SmallSerineSerTCA,TCG,TCC,TCT,Polar5.68TinyAGT,AGCThreonineThrACT,ACC,ACA,ACGPolar5.60SmallTryptophanTrpTCCNonpolar5.89NormalTyrosineTyrTAC,TATPolar5.64NormalValineValGTA,GTG,GTC,GTTNonpolar6.00SmallStop codonTermTAA,TAG,TGA———From Table5, we can see all codons are divided into different sets with each set corresponding to one amino acid. The size and the elements in one codon set vary from one amino acid to another. This implies that the costs of driving one codon to the desired final set generated by the desired final state might be different from the costs of driving the complementary codon to the desired final set generated by the complementary of desired final state. More discussions about this issue are presented in Section 3.The distance reference between any two codons can be defined by a weighted sum of the differences between physical and chemical properties or other reasonable functions. And the distance reference between two DNA sequences is defined as the sum of distance reference between the corresponding pair of codons. The biological statics plays a crucial rule to define this distance function in practical.An example of the distance function can be expressed as(20)d(ξ1,ξ2)=ζpolaritypolarity(ξ1,ξ2)+ζPHPH(ξ1,ξ2)+ζsizesize(ξ1,ξ2),polarity(ξ1,ξ2)={0ifξ1,ξ2arebothpolarornon-polar,1ifoneofξ1,ξ2ispolar,andtheothernon-polar,PH(ξ1,ξ2)=|PHvalueofξ1-PHvalueofξ2|,size(ξ1,ξ2)={0,ifξ1,ξ2arebothtiny,small,ornormal,σ1,ifoneofξ1,ξ2istiny,andtheothersmall,σ2,ifoneofξ1,ξ2istiny,andtheothernormal,σ3,ifoneofξ1,ξ2issmall,andtheothernormal,
where ξ1,ξ2 are two amino acids.d(ξ1,ξ2) is then assigned to d(φ1,φ2) with φ1,φ2 corresponding to amino acids ξ1,ξ2, respectively.Since the generalized optimal control problem in (18) and (19) is a multistage problem that can be broken down into simpler steps at different time points. Therefore, we can solve it by dynamic programming.For dynamic programming, the optimal control policy is constructed backward. And Bellman’s principle of optimality states that the optimal policy forx0 to {xNd} is also the optimal policy for the tail problem, from xq to {xNd}.The tail problem is defined as(21)Jq(xq)=min⁡{u,c,v,c′}q,q+1,…,N-1{∑k=qN-1∑l1=1l∑i=1nαl1uk,l1i+∑k=qN-1∑l2=1m∑i=1nβl2ck,l2i+∑k=qN-1∑l3=1l′∑i=1nαl3′vk,l3i+∑k=qN-1∑l4=1m′∑i=1nβl4′bk,l4i+∑k=qNE{h,r,h′,r′}q,q+1,…,N-1[d(xk,{xkd})]}.
The iterative update equation to find optimal policy can be expressed by (22), according to the dynamic programming algorithm in [31]. 								(22)JN(xN)=d(xN,{xNd}),Jq(xq)=min⁡uq,cq,vq,cq′Ehq,rq,hq′,rq′[∑l1=1l∑i=1nαl1uq,l1i+∑l2=1m∑i=1nβl2cq,l2i+∑l3=1l′∑i=1nαl3′vq,l3i+∑l4=1m′∑i=1nβl4′cq,l4′i+d(xq,{xqd})+Jq+1(xq+1)]=min⁡uq,cq,vq,cq′{∑l1=1l∑i=1nαl1uq,l1i+∑l2=1m∑i=1nβl2cq,l2i+∑l3=1l′∑i=1nαl3′vq,l3i+∑l4=1m′∑i=1nβl4′cq,l4′i+Ehq,rq,hq′,rq′[d(xq,{xqd})+Jq+1(xq+1)]∑l1=1l∑i=1nαl1uq,l1i},q=0,1,…,N-1.


## 3. Results and Discussion

In the following examples, we consider applying chemical mutagens only because the randomness of applying radiation is much larger and more difficult to control. We also omit the mutations between a normal base andO because of high-cost βl2 and the unavailable chemical and physical properties for codons containing O.The distance reference between codons used in Sections3.2 and 3.3 is computed by (20) with ζpolarity=8, ζPH=3, ζsize=1, σ1=2, σ2=5 and σ3=3. We only keep the final penalty but omit the off-trajectory penalty along the trajectory.


### 3.1. Base-to-Base, Deterministic Optimal Control Problem

We define the distance reference between bases as(23)d(ψ1,ψ2)={0,ifxN=xNd,∞,ifxN≠xNd,
with ψ1,ψ2∈ℱ∖{0}, where ℱ∖{0} denotes the set ℱ excluding the element 0.Our optimal control problem for point mutations is(24)J0(x0)=min⁡uk,l1,0≤k≤N,1≤l1≤l{∑k=0N-1∑l1=1lαl1uk,l1},
subject to(25)xk+1=(-I+∑l1=1luk,l1s)xk,xN=xNd,
with x0 given, xk∈ℱ∖{0}.Suppose that there are 12 kinds of mutagens (l1=12), each corresponding to a specific transfer pattern as in Table 6, all available controls and the respective costs can be immediately listed as in Tables 7 and 8.Table 6
An example of chemical mutagens and their corresponding transfer patterns in deterministic mutations.
Index(l1)123456Transfer patternA→AA→GA→CG→AG→GG→TIndex(l1)789101112Transfer patternC→AC→CC→TT→GT→CT→TTable 7
Controls corresponding to transfer between bases within one step. The leftmost column denotes the statekth step, and the upmost row denotes the (k+1)th state.
kth(k+1)thAGCAuAAuAGuACuATGuGAuGGuGCuGTCuCAuCGuCCuCTTuTAuTGuTCuTTTable 8
Corresponding step cost of controls as shown in Table7.
kth(k+1)thAGCAαAAαAGαACαATGαGAαGGαGCαGTCαCAαCGαCCαCTTαTAαTGαTCαTTThe elements along the antidiagonal of Table7, uAT,uGC,uCG, and uTA, are artificially added, because the complementary transfers naturally happen and no mutagen is necessary. Thus, the costs along the antidiagonal of Table 8 are all zero, that is, αAT=αGC=αCG=αTC=0. The equivalence relationship between subscription in two nucleotide bases and subscription in integer l1(ψ1ψ2):{A,T,G,C}×{A,T,G,C}→{integersfrom1to12} is defined by Table 6.Under the above assumptions, we can write update equation for optimal policy explicitly as(26)Jq(xq)=min⁡uq,l1{αxqψ+Jq+1(ψ),∀ψ∈{A,T,G,C}⟺F∖{0}}.Proposition 7.
For the samexNd, Jq(ψ)≤Jq+1(ψ¯),  for all q, 0≤q≤N-1,  for all ψ∈{A,T,G,C}, where ψ¯ denotes the complementary base of ψ. If, in addition, the system is completely controllable, ∃M, s.t. JM(ψ) is the global minimum and for all q≤M, Jq(ψ)=JM(ψ¯) if M-q≡ 1 (mod 2), and Jq(ψ)=JM(ψ) if M-q≡0 (mod 2). In our example, M≥N-6.Proof.
This first part is due to the zero cost for the transfers between complementary bases in the consecutive steps.
For any0≤q≤N-1, the relationship between minimal costs in consecutive steps is shown in (26). Since ψ∈{A,T,G,C}, αψψ¯+Jq+1(ψ¯) is one of the four elements in the set from which the Jq(ψ) is picked. Moreover, αψψ¯=0. Therefore, Jq+1(ψ¯) is one of the four elements in the set. Since Jq(ψ) is the minimum picking for a set containing Jq+1(ψ¯), we conclude that Jq(ψ)≤Jq+1(ψ¯).
TheM value in our example is proved by brute force method, that is, JN-6(ψ) is a guaranteed global minimum. For completely controllable systems, this M always exists. The existence of M implies that for without limitation in the number of steps, we can reach the global optimal in N-M steps, 6 steps in our example.
Suppose thatq=M, JM(ψ) is the global minimum, thus JM-1(ψ¯)≥JM(ψ). However,JM-1(ψ¯)≤JM(ψ) according to the first part of the proposition. Therefore, JM-1(ψ¯)=JM(ψ) for the same xNd. Therefore, JM-1(ψ¯) is also a global minimum.
By backward induction, suppose forq=q1, the statement is true, that is, Jq1-1(ψ¯)=Jq1(ψ) is the global optimal either from xq1-1=ψ¯ or xq1=ψ to xNd. Obviously, for q=q1-1, the statement is still true. Therefore, Jq(ψ)=Jq-2(ψ)=Jq-1(ψ¯),ψ∈{A,T,G,C}, for all q,  2≤q≤M.In the proof of global minimum that can be reached in the finite step in Proposition7, we also discover Proposition 8. Here, Jq(xq,xNd) denotes the optimal cost from xq to xNd.Proposition 8.
Given two single base mutation optimal control problems, with the same fixedN, with and desired final states complementary to each other. If JM(ψ,xNd) is the global minimum, then JM(ψ¯,xNd¯) is also the global minimum, that is, the global minimum of both systems is reach at the same stage M. Moreover, for all q,  0≤q≤M,
(27)Jq(ψ,xNd)=Jq(ψ¯,xNd¯),ψ,xNd∈{A,T,G,C}.Physically, Proposition8 states that the optimal can be achieve at the same step from a pair of complementary bases to another pair of complementary bases at the same cost. However, this fact is true only for base-to-base deterministic mutations, because the distance reference is well defined by (23).Now, we show an example with simulation results.The costs of applying different mutagens are listed in Table9. It is a numerical assignment to Table 8. Since we apply mutagens before the replication starts, uAA actually transfer A to T and then to A by replication. For simplicity, we just use the kth and (k+1)th step states as subscripts to represent the corresponding control and cost. The costs of transitions are lower than the costs of transversions. Therefore, αAC,αCA,αGT,αTG is smaller than other mutagens, except artificial ones.Table 9
Sample step costs.
xkxk+1AGCTA5.216.602.330G6.158.9503.82C4.6109.177.24T00.645.0910.28If we useχ to denote the costs of mutagens as listed in Table 9, then (28)χ=[5.216.602.3306.158.9503.824.6109.177.2400.645.0910.28]=[αAAαAGαACαATαGAαGGαGCαGTαCAαCGαCCαCTαTAαTGαTCαTT]⟺[α1α2α30α4α50α6α70α8α90α10α11α12].Running the dynamic programming for every pair of(xq,xNd)∈{A,T,G,C}×{A,T,G,C}, N=9. Here, we sightly modify our notation. We use Jq(xq,xNd) to denote the optimal cost from xq to xNd. Then, (29)Jq=[Jq(A,A)Jq(A,G)Jq(A,C)Jq(A,T)Jq(G,A)Jq(G,G)Jq(G,C)Jq(G,T)Jq(C,A)Jq(C,G)Jq(C,C)Jq(C,T)Jq(T,A)Jq(T,G)Jq(C,C)Jq(T,T)].The path to reach the optimal cost is denoted by(30)Pq=[Pq(A,A)Pq(A,G)Pq(A,C)Pq(A,T)Pq(G,A)Pq(G,G)Pq(G,C)Pq(G,T)Pq(C,A)Pq(C,G)Pq(C,C)Pq(C,T)Pq(T,A)Pq(T,G)Pq(T,C)Pq(T,T)],wherePq(xq,xNd) is the (q+1)th state from xq to xNd, that is, xq+1=Pq(xq,xNd).The simulation results are shown as below, including optimal costs for all possible transfer pairs(xq,xNd)∈{A,T,G,C}×{A,T,G,C}, Jq, 0≤q≤8, graphical representation in Figure 6, and optimal path Pq,  0≤q≤7. 								(31)J0=[5.215.090.6406.156.7903.823.8206.796.1500.645.095.21],J1=[00.645.095.213.8206.796.156.156.7903.825.215.090.640],J2=[5.215.090.6406.156.7903.823.8206.796.1500.645.095.21],J3=[00.645.095.213.8206.796.156.156.7903.825.215.090.640],J4=[5.215.090.6406.156.7903.823.8206.796.1500.645.095.21],J5=[00.645.095.213.8206.796.156.156.7903.825.215.090.640],J6=[5.215.090.6406.156.7903.823.8207.886.1500.645.095.21],J7=[00.645.095.213.8208.486.156.157.8803.825.215.090.640],J8=[5.216.602.3306.158.9503.824.6109.177.2400.645.0910.28],(32)P0=[A,TTTTA,CA,CCC,TGGGGAA,GA,CA],P1=[TTTA,TC,TCA,CA,CGGGGAA,CA,GA],P2=[A,TTTTA,CA,CCC,TGGGGAA,GA,CA],P3=[TTTA,TC,TCA,CA,CGGGGAA,CA,GA],P4=[A,TTTTA,CA,CCC,TGGGGAA,GA,CA],P5=[TTTA,TC,TCAA,CGGGGAA,CA,GA],P6=[A,TTTTA,CACC,TGGTGAA,GA,CA],P7=[TTTATCAAGTGGACGA].Figure 6
Graphical representation ofJq,  0≤q≤8,  N=9, in single-base deterministic mutation example. The x-axis and y-axis represent xq and xNd, respectively. Jq(xq,xNd) are represented by 16 isolated points. Those discrete points are connected together to show the surface.For simplicity, we use 1 to representA, 2 to G, 3 to C, and 4 to T in graphical interpretation. From Figure 6, we can see clearly that the optimal cost decreases as q decreases in the first few steps, and then optimal cost remains at the global minimum. This phenomenon obeys Proposition 7. In this example, global optimal is reached at M=5 for all pairs of initial and final states as J7≠J5=J3 and J6≠J4=J2. So, the global minimum is achieved before we reach N-5=4 in this particular case. This also implies that with N free we can reach desired final state in 4 steps from given initial state.Observing closely toJq,  0≤q≤5, Jq-1 equals to Jq by exchanging the first and the last columns, and the second and the third columns, which is consistent with Proposition 8. Or we can exchange the first and the last rows, and the second and the third rows of Jq to obtain Jq-1. Jq1 and Jq2 are the same for q1,q2≤M=5 for q1-q2=0 (mod 2). This obeys Proposition 7.The optimal trajectories are generated fromPq(xq,xNd). For example, given x2=T, and the final state x9d=G, we want to generate the optimal trajectories.x3=P2(T,G)=A,G.Ifx3=A, x4=P3(A,G)=T;  if x3=G, x4=P3(A,G)=C.Ifx4=T, x5=P4(T,G)=A,G; if x4=C, x5=P4(C,G)=G.Ifx5=A, x6=P5(A,G)=T; if x5=G, x6=P5(G,G)=C.Ifx6=T, x7=P6(T,G)=A,G; if x6=C, x7=P6(C,G)=G.Ifx7=A, x8=P7(A,G)=T; if x7=G, x8=P7(G,G)=C.So, the optimal routes areT→A→T→A→T→A→T→αTGuTGG;T→A→T→A→T→αTGuTGG→C→G;T→A→T→αTGuTGG→C→G→C→G;T→αTGuTGG→C→G→C→G→C→G.Consequently, the optimal cost isJ2(T,G)=0.64=αTG.Optimal trajectories for other pairs of initial and final states can be obtained in the same manner.It takes less than 1 second to generate optimal trajectories for all pairs of initial and final states withN=9 on a regular desktop. Since we have already proven by the brute force method that the global optimal can be achieved with M≤6, the computation time can be further reduced by taking N=6 with all the results necessary for this example.


### 3.2. Codon-to-Codon, Deterministic Optimal Control Problem

For codon-to-codon deterministic mutations, we formulate our optimal control problem as(33)J0(x0)=min⁡uk,l1i,0≤k≤N,1≤l1≤l,1≤i≤3{∑k=0N-1∑l1=1l∑i=13αl1uk,l1i+d(xN,{xNd})},
subject to(34)xk+1=(-I+∑l1=1l∑i=13uk,l1iseieiT)xk,
with x0,xNd∈ℱ∖{0}3 given, xk∈ℱ∖{0}3, for all k,  0≤k≤N, and d(φ1,φ2),φ1,φ2∈ℱ∖{0}3 as defined in Section 2.3.If we take the same assumption on available mutagens as in Section3.1, then we can write update equation for optimal control policy explicitly as
(35)Jq(xq)=min⁡uq,l1i,1≤l1≤l,1≤i≤3×{αxq1ψ1+Jq+1([ψ1xq2¯xq3¯]),αxq2ψ2+Jq+1([xq1¯ψ2xq3¯]),αxq3ψ3+Jq+1([xq1¯xq2¯ψ3]),([xq1¯xq2¯ψ3]),ψ1,ψ2,ψ3∈{A,T,G,C}⟺F∖{0}},
with xqi∈{A,T,G,C}⇔ℱ∖{0},  1≤i≤3 denotes the ith element of xq∈ℱ∖{0}3, and xqi¯ denotes the complementary base of xqi.The optimal control sequences depend on the numerical values ofαl1s and d(φ1,φ2), φ1,φ2∈ℱ∖{0}3. Though we do not have real values of αl1s and d(φ1,φ2), we can always obtain simulation results to compare the result differences by assigning different numerical values to those parameters.Therefore, we use three different assignments ofαl1s and the same d(φ1,φ2) to generate our simulation results. Those three assignments of αl1s are χ,  5χ, and 0.5χ, respectively, with χ the same as assigned in Section 3.1. In every particular example, it takes approximately 2 seconds on a regular desktop to generate the optimal path table for all pairs of initial and final states with N=19, and the dynamic programming algorithm ensures that the optimal control for tail problems is generated at the same time.The graphical interpretation of three assignments are shown in Figures7, 8, and 9, respectively. The x-axis and y-axis denote xq and xNd, respectively. For a codon [ψ1ψ2ψ3]T,ψ1,ψ2, ψ3∈{A,T,G,C}⇔ℱ∖{0}, its index is calculated by
(36)42(ψ1-1)+4(ψ2-1)+ψ3,where ψi=1 if A, ψi=2 if G, ψi=3 if C, and ψi=4 if T, 1≤i≤3, for the simplicity of graphical interpretation. A codon has 43=64 combinations. Thus, there are 642 pairs of initial and final desired states, and there are 64×21 pairs of initial state and final desired set. The surface is generated by connecting 64×64 discrete points together. Jq is calculated following the same procedure as in base-to-base deterministic cases. The value of optimal cost can be read directly from graphical interpretation, and the optimal path can be generated from path matrix Pq, similar to base-to-base deterministic case. Both Jq and Pq,forallq,0≤q≤N, are of 64×64 dimension.Figure 7
Graphical representation ofJq,  q=0,1,2,3,15,16,17,18 for codon-to-codon deterministic mutations, with αl1=χ, N=19.Figure 8
Graphical representation ofJq,  q=0,1,2,3,15,16,17,18 for codon-to-codon deterministic mutations, with αl1=5χ, N=19.Figure 9
Graphical representation ofJq,  q=0,1,2,3,15,16,17,18 for codon-to-codon deterministic mutations, with αl1=0.5χ, N=19.From the graphical interpretation and Table10, we find that the value of q where the global minimum is reached at the first time decreases as αl1 decreases. And J0 is more similar to J18 with α=5χ than with α=χ or α=0.5χ. This implies that if d(φ1,φ2) are the deterministic term in our objective function, then the treatment plan is made to drive the final state as close to the desired set as possible; if the costs of applying mutagens is the deterministic term in the objective function, then the treatment plan tends to stay in the original state and applies less mutagens; if they are of equal weight, then the treatment plan deals with this tradeoff.Table 10
Simulation results with differentαl1 assignments and the first q where the global optimal is reached.
Nαl1M (first q when the  global minimum is reached)Global minimum190.5χq=12J12(x12)19χq=13J13(x13)195χq=15J15(x15)Moreover, no matter how the numerical values of final penalty and the costs of applying mutagens changes in our objective function as shown in (33), there is always a M, M≤N-18, JM(xM) is global minimum. Proposition 7 can be extended to codon-to-codon deterministic mutations as stated in Proposition 9.Proposition 9.
Given an optimal control problem with objective function in the form of (33), constraints in the form of (34), and all available chemical mutagens and their corresponding transfer pairs and costs as listed in Tables 6, 7, and 8, Jq(φ¯)≤Jq+1(φ),φ∈ℱ∖{0}3. If, in addition, the system is completely controllable, ∃M, s.t. JM(φ) is the global minimum and for all q≤M, Jq(φ)=JM(φ¯) if M-q≡1 (mod 2), and Jq(φ)=JM(φ) if M-q≡0 (mod2). In our example, M≥N-18.Proof.
We only prove thatM≥N-18 here since the rest is similar to the proof of Proposition 7.
The objective function in (33) can be written as the summation of three separate single-base mutation systems and the distance reference between final states and the final desired set, that is, (37)Jq(xq)=min⁡n1,n2,n3≥02N≤n1+n2+n3≤3N-1{Jn1(xn-1)(xq1,ψ1)︸optimalcostsofbase-to-basedeterministicoptimalcontrolproblemformedbythe1stbase+Jn2(xn-2)(xq2,ψ2)︸optimalcostsofbase-to-basedeterministicoptimalcontrolproblemformedby-the2ndbase+Jn3(xn-3)(xq3,ψ3)︸optimalcostsofbase-to-basedeterministicoptimalcontrolproblemformedbythe3rdbase+d([ψ1ψ2ψ3],{xNd})Jn1(xn-1)(xq1,ψ1)︸optimalcostsofbase-to-basedeterministicoptimalcontrolproblemformedbythe1stbase},
where N-q=(N-n1)+(N-n2)+(N-n3).
According to Proposition7, JN-6(xN-6) is guaranteed to be the global optimal for single-base mutations. Therefore, optimal costs corresponding to three single-base mutation systems, Jn1(xn-1)(xq1,ψ1), Jn2(xn-2)(xq2,ψ2), Jn3(xn-3)(xq3,ψ3) are guaranteed to reach their own global optimal at n1=n2=n3=N-6 with all possible combinations of ψ1,ψ2,ψ3∈{A,T,G,C}. Therefore, q=N-18 is a guaranteed global optimal.Indeed, Proposition9 is a quite loose condition. In real examples above, the M value where the first global optimal is reached at earlier stage as listed in Table 10.Graphically, the indices of complementary codonsφ1=[ψ1ψ2ψ3]T and φ2=[ψ1¯ψ2¯ψ3¯]T sum up to 65, that is,  (38)(16(ψ1-1)+4(ψ2-1)+ψ3)+(16(ψ1¯-1)+4(ψ2¯-1)+ψ3¯)=(16(ψ1-1)+4(ψ2-1)+ψ3)+(16((5-ψ1)-1)+4((5-ψ2)-1)+(5-ψ3))=65.Therefore,Jq and Jq-1, 1≤q≤M are symmetric about the plane x=32.5, Jq-2 and Jq, 2≤q≤M are the same, as shown in Figures 7, 8, and 9.However, Proposition8 cannot be extended to codon-to-codon deterministic case due to the redundancy of genetic codes, that is, the set of codons translated to the same amino acid varies from one amino acid to another as shown in Table 5. The simulation results show that the costs, a pair of complementary codons, to two final desired set generated by a pair of complementary final desired codons are different, that is, Jq(xq,{xNd})≠Jq(xq¯,{xNd¯}), in general, for any q. Graphically, the optimal cost profile Jq is not symmetric about the plane y=32.5 for Jq-1, 1≤q≤M. Therefore, the doctors need to pick the strand with lower cost to make the treatment plan. This also implies that in large-scale cases, for instance, a gene containing hundreds of nucleotide bases, the doctors should make the treatment plan based on the strand the total cost of which is lower than the other.


### 3.3. Codon-to-Codon, Stochastic Optimal Control Problem

The optimal control problem of codon-to-codon stochastic mutations can be written as(39)J0(x0)=min⁡uk,l1i,0≤k≤N-11≤l1≤l,1≤i≤3{+Ehk,l1i,0≤k≤N-11≤l1≤l,1≤i≤3[d(xN,{xNd})][d(xN,{xNd})]∑k=0N-1∑l1=1l∑i=13αl1uk,l1i+Ehk,l1i,0≤k≤N-11≤l1≤l,1≤i≤3[d(xN,{xNd})]},
subject to(40)xk+1=-Ixk+∑l1=1l∑i=13uk,l1ihk,l1iseieiTxk,
with x0,xNd∈ℱ∖{0}3 given, xk∈ℱ∖{0}3.The major difference between deterministic and stochastic systems is that we impose the random binary vector,hk,l1i, in our system equation (40). We denote the probability associated with hk,l1i,j to be  pl1,ψ1ψ2(h) with ψ1,ψ2∈{A,T,G,C}. The equivalence relationship between j and ψ1ψ2 can be found in Table 3.Again, we assume that we havel1=12 kinds of mutagens, each corresponding to one major transfer pattern, associated with probability assignments, as listed in Table 11.Table 11
12 kinds of mutagens, each corresponding to major transfer patterns, and the probabilities of different mutagens on different transfer patterns.
Index (l1)FromToMajor transfer patternAGCT1Ap1,AA(h)p1,AG(h)p1,AC(h)p1,AT(h)A→A2Ap2,AA(h)p2,AG(h)p2,AC(h)p2,AT(h)A→G3Ap3,AA(h)p3,AG(h)p3,AC(h)p3,AT(h)A→C4Gp4,GA(h)p4,GG(h)p4,GC(h)p4,GT(h)G→A5Gp5,GA(h)p5,GG(h)p5,GC(h)p5,GT(h)G→G6Gp6,GA(h)p6,GG(h)p6,GC(h)p6,GT(h)G→T7Cp7,CA(h)p7,CG(h)p7,CC(h)p7,CT(h)C→A8Cp8,CA(h)p8,CG(h)p8,CC(h)p8,CT(h)C→C9Cp9,CA(h)p9,CG(h)p9,CC(h)p9,CT(h)C→T10Tp10,TA(h)p10,TG(h)p10,TC(h)p10,TT(h)T→G11Tp11,TA(h)p11,TG(h)p11,TC(h)p11,TT(h)T→C12Tp12,TA(h)p12,TG(h)p12,TC(h)p12,TT(h)T→TThen, we can write updated formula for optimal control policy explicitly as(41)Jq(xq)=min⁡uq,l1i,1≤l1≤l,1≤i≤3{αxq1ψ1+Ehq,l1(xq1ψ1)1[Jq+1([⋅xq2¯xq3¯])],αxq2ψ2+Ehq,l1(xq2ψ2)2[Jq+1([xq1¯⋅xq3¯])],αxq3ψ3+Ehq,l1(xq3ψ3)3[Jq+1([xq1¯xq2¯⋅])],ψ1,ψ2,ψ3∈{A,T,G,C}⟺F∖{0}xq1¯xq2¯⋅}
where(42)Ehq,l1(xq1ψ1)1[Jq+1([⋅xq2¯xq3¯])]=pl1(xq1ψ1),xq1A(h)[Jq+1([Axq2¯xq3¯])]+pl1(xq1ψ1),xq1G(h)[Jq+1([Gxq2¯xq3¯])]+pl1(xq1ψ1),xq1C(h)[Jq+1([Cxq2¯xq3¯])]+pl1(xq1ψ1),xq1T(h)[Jq+1([Txq2¯xq3¯])],wherexqi∈{A,T,G,C}⇔ℱ∖{0},  0≤q≤N-1, 1≤i≤3 denotes the ith element of xq∈ℱ∖{0}3, xqi¯ denotes the complementary base of xqi, and l1:ψ1ψ2∈{A,T,G,C}×{A,T,G,C}→{integersfrom1to12}, the mapping from major transfer pattern ψ1→ψ2 to mutagen index, as shown in Table 11. The mathematical expression of 𝔼hq,l1(xq2ψ2)2[Jq+1([xq1¯·xq3¯]T)] and 𝔼hq,l1(xq3ψ3)3[Jq+1([xq1¯xq2¯·]T)] is similar to 𝔼hq,l1(xq1ψ1)1[Jq+1([·xq2¯xq3¯]T)] as shown above.In order to run the simulation, we assign numerical values to probabilities in Table11, as illustrated in Table 12.Table 12
Sample probabilities with respect to different mutagens and different transfer patterns.
Index (l1)FromToAGCT1A0.900.050.030.022A0.110.580.210.103A0.140.160.420.284G0.850.070.030.055G0.020.020.920.046G0.100.090.220.597C0.790.130.040.048C0.010.020.870.109C0.040.120.090.7510T0.130.760.050.0611T0.070.030.620.2812T0.080.040.250.63As in Section3.2, we use three different assignments for αl1s, χ,5χ, and 0.5χ, respectively. The optimal cost profile Jq with selected q values, for every pair of (xq,xNd), is graphically interpreted in Figures 10, 11, and 12, respectively, with N=29, with computation time of approximately 7 seconds on a regular desktop.Figure 10
Graphical representation ofJq(xq),  q=0,1,2,3,15,19,24,28 for codon-to-codon stochastic mutations, with αl1=χ, probability assignment as in Table 12, N=29.Figure 11
Graphical representation ofJq(xq),  q=0,1,2,3,15,19,24,28 for codon-to-codon stochastic mutations, with αl1=5χ, N=29.Figure 12
Graphical representation ofJq(xq),  q=0,1,2,3,15,19,24,28 for codon-to-codon stochastic mutations, with αl1=0.5χ, N=29.The simulation results are similar to those in Section3.2. The profile of J0 is more similar to J29 when αl1s are assigned 5χ than χ or 0.5χ. This implies in codon-to-codon stochastic mutations; the optimal control sequence behaves as codon-to-codon deterministic cases, that is, the system tends to getting as close as possible to the final desired set if αl1s are much smaller than d(φ1,φ2), and the system tends to remain in the same state with minor mutations when αl1s are relatively larger than d(φ1,φ2),φ1,φ2∈ℱ∖{0}3. And Jq(φ¯)≤Jq+1(φ),  φ∈ℱ∖{0}3 is still valid in codon-to-codon stochastic case.However, for stochastic cases, we cannot reach a global minimum because of the randomness caused by mutagens. Since, in usual cases, there exists no stationary global minimum, we need to define error toleranceϵ, that is, if |Jq(ψ)-Jq-1(ψ¯)|≤ϵ with the same {xNd}, then we can stop at Jq(xq). Otherwise, we need to proceed to calculate Jq-2(ψ). The value of ϵ is decided based on the doctors, experience. Obviously, the smaller ϵ is, the better treatment plan is. However, we can still observe Figures 10, 11, and 12 to conclude that J0 and J2 are almost of the same shape in all three different parameter assignments. Higher dimensional optimal control problems, gene-to-gene stochastic mutations, can be solved as a series of cascade codon-to-codon stochastic problems.


## 3.1. Base-to-Base, Deterministic Optimal Control Problem

We define the distance reference between bases as(23)d(ψ1,ψ2)={0,ifxN=xNd,∞,ifxN≠xNd,
with ψ1,ψ2∈ℱ∖{0}, where ℱ∖{0} denotes the set ℱ excluding the element 0.Our optimal control problem for point mutations is(24)J0(x0)=min⁡uk,l1,0≤k≤N,1≤l1≤l{∑k=0N-1∑l1=1lαl1uk,l1},
subject to(25)xk+1=(-I+∑l1=1luk,l1s)xk,xN=xNd,
with x0 given, xk∈ℱ∖{0}.Suppose that there are 12 kinds of mutagens (l1=12), each corresponding to a specific transfer pattern as in Table 6, all available controls and the respective costs can be immediately listed as in Tables 7 and 8.Table 6
An example of chemical mutagens and their corresponding transfer patterns in deterministic mutations.
Index(l1)123456Transfer patternA→AA→GA→CG→AG→GG→TIndex(l1)789101112Transfer patternC→AC→CC→TT→GT→CT→TTable 7
Controls corresponding to transfer between bases within one step. The leftmost column denotes the statekth step, and the upmost row denotes the (k+1)th state.
kth(k+1)thAGCAuAAuAGuACuATGuGAuGGuGCuGTCuCAuCGuCCuCTTuTAuTGuTCuTTTable 8
Corresponding step cost of controls as shown in Table7.
kth(k+1)thAGCAαAAαAGαACαATGαGAαGGαGCαGTCαCAαCGαCCαCTTαTAαTGαTCαTTThe elements along the antidiagonal of Table7, uAT,uGC,uCG, and uTA, are artificially added, because the complementary transfers naturally happen and no mutagen is necessary. Thus, the costs along the antidiagonal of Table 8 are all zero, that is, αAT=αGC=αCG=αTC=0. The equivalence relationship between subscription in two nucleotide bases and subscription in integer l1(ψ1ψ2):{A,T,G,C}×{A,T,G,C}→{integersfrom1to12} is defined by Table 6.Under the above assumptions, we can write update equation for optimal policy explicitly as(26)Jq(xq)=min⁡uq,l1{αxqψ+Jq+1(ψ),∀ψ∈{A,T,G,C}⟺F∖{0}}.Proposition 7.
For the samexNd, Jq(ψ)≤Jq+1(ψ¯),  for all q, 0≤q≤N-1,  for all ψ∈{A,T,G,C}, where ψ¯ denotes the complementary base of ψ. If, in addition, the system is completely controllable, ∃M, s.t. JM(ψ) is the global minimum and for all q≤M, Jq(ψ)=JM(ψ¯) if M-q≡ 1 (mod 2), and Jq(ψ)=JM(ψ) if M-q≡0 (mod 2). In our example, M≥N-6.Proof.
This first part is due to the zero cost for the transfers between complementary bases in the consecutive steps.
For any0≤q≤N-1, the relationship between minimal costs in consecutive steps is shown in (26). Since ψ∈{A,T,G,C}, αψψ¯+Jq+1(ψ¯) is one of the four elements in the set from which the Jq(ψ) is picked. Moreover, αψψ¯=0. Therefore, Jq+1(ψ¯) is one of the four elements in the set. Since Jq(ψ) is the minimum picking for a set containing Jq+1(ψ¯), we conclude that Jq(ψ)≤Jq+1(ψ¯).
TheM value in our example is proved by brute force method, that is, JN-6(ψ) is a guaranteed global minimum. For completely controllable systems, this M always exists. The existence of M implies that for without limitation in the number of steps, we can reach the global optimal in N-M steps, 6 steps in our example.
Suppose thatq=M, JM(ψ) is the global minimum, thus JM-1(ψ¯)≥JM(ψ). However,JM-1(ψ¯)≤JM(ψ) according to the first part of the proposition. Therefore, JM-1(ψ¯)=JM(ψ) for the same xNd. Therefore, JM-1(ψ¯) is also a global minimum.
By backward induction, suppose forq=q1, the statement is true, that is, Jq1-1(ψ¯)=Jq1(ψ) is the global optimal either from xq1-1=ψ¯ or xq1=ψ to xNd. Obviously, for q=q1-1, the statement is still true. Therefore, Jq(ψ)=Jq-2(ψ)=Jq-1(ψ¯),ψ∈{A,T,G,C}, for all q,  2≤q≤M.In the proof of global minimum that can be reached in the finite step in Proposition7, we also discover Proposition 8. Here, Jq(xq,xNd) denotes the optimal cost from xq to xNd.Proposition 8.
Given two single base mutation optimal control problems, with the same fixedN, with and desired final states complementary to each other. If JM(ψ,xNd) is the global minimum, then JM(ψ¯,xNd¯) is also the global minimum, that is, the global minimum of both systems is reach at the same stage M. Moreover, for all q,  0≤q≤M,
(27)Jq(ψ,xNd)=Jq(ψ¯,xNd¯),ψ,xNd∈{A,T,G,C}.Physically, Proposition8 states that the optimal can be achieve at the same step from a pair of complementary bases to another pair of complementary bases at the same cost. However, this fact is true only for base-to-base deterministic mutations, because the distance reference is well defined by (23).Now, we show an example with simulation results.The costs of applying different mutagens are listed in Table9. It is a numerical assignment to Table 8. Since we apply mutagens before the replication starts, uAA actually transfer A to T and then to A by replication. For simplicity, we just use the kth and (k+1)th step states as subscripts to represent the corresponding control and cost. The costs of transitions are lower than the costs of transversions. Therefore, αAC,αCA,αGT,αTG is smaller than other mutagens, except artificial ones.Table 9
Sample step costs.
xkxk+1AGCTA5.216.602.330G6.158.9503.82C4.6109.177.24T00.645.0910.28If we useχ to denote the costs of mutagens as listed in Table 9, then (28)χ=[5.216.602.3306.158.9503.824.6109.177.2400.645.0910.28]=[αAAαAGαACαATαGAαGGαGCαGTαCAαCGαCCαCTαTAαTGαTCαTT]⟺[α1α2α30α4α50α6α70α8α90α10α11α12].Running the dynamic programming for every pair of(xq,xNd)∈{A,T,G,C}×{A,T,G,C}, N=9. Here, we sightly modify our notation. We use Jq(xq,xNd) to denote the optimal cost from xq to xNd. Then, (29)Jq=[Jq(A,A)Jq(A,G)Jq(A,C)Jq(A,T)Jq(G,A)Jq(G,G)Jq(G,C)Jq(G,T)Jq(C,A)Jq(C,G)Jq(C,C)Jq(C,T)Jq(T,A)Jq(T,G)Jq(C,C)Jq(T,T)].The path to reach the optimal cost is denoted by(30)Pq=[Pq(A,A)Pq(A,G)Pq(A,C)Pq(A,T)Pq(G,A)Pq(G,G)Pq(G,C)Pq(G,T)Pq(C,A)Pq(C,G)Pq(C,C)Pq(C,T)Pq(T,A)Pq(T,G)Pq(T,C)Pq(T,T)],wherePq(xq,xNd) is the (q+1)th state from xq to xNd, that is, xq+1=Pq(xq,xNd).The simulation results are shown as below, including optimal costs for all possible transfer pairs(xq,xNd)∈{A,T,G,C}×{A,T,G,C}, Jq, 0≤q≤8, graphical representation in Figure 6, and optimal path Pq,  0≤q≤7. 								(31)J0=[5.215.090.6406.156.7903.823.8206.796.1500.645.095.21],J1=[00.645.095.213.8206.796.156.156.7903.825.215.090.640],J2=[5.215.090.6406.156.7903.823.8206.796.1500.645.095.21],J3=[00.645.095.213.8206.796.156.156.7903.825.215.090.640],J4=[5.215.090.6406.156.7903.823.8206.796.1500.645.095.21],J5=[00.645.095.213.8206.796.156.156.7903.825.215.090.640],J6=[5.215.090.6406.156.7903.823.8207.886.1500.645.095.21],J7=[00.645.095.213.8208.486.156.157.8803.825.215.090.640],J8=[5.216.602.3306.158.9503.824.6109.177.2400.645.0910.28],(32)P0=[A,TTTTA,CA,CCC,TGGGGAA,GA,CA],P1=[TTTA,TC,TCA,CA,CGGGGAA,CA,GA],P2=[A,TTTTA,CA,CCC,TGGGGAA,GA,CA],P3=[TTTA,TC,TCA,CA,CGGGGAA,CA,GA],P4=[A,TTTTA,CA,CCC,TGGGGAA,GA,CA],P5=[TTTA,TC,TCAA,CGGGGAA,CA,GA],P6=[A,TTTTA,CACC,TGGTGAA,GA,CA],P7=[TTTATCAAGTGGACGA].Figure 6
Graphical representation ofJq,  0≤q≤8,  N=9, in single-base deterministic mutation example. The x-axis and y-axis represent xq and xNd, respectively. Jq(xq,xNd) are represented by 16 isolated points. Those discrete points are connected together to show the surface.For simplicity, we use 1 to representA, 2 to G, 3 to C, and 4 to T in graphical interpretation. From Figure 6, we can see clearly that the optimal cost decreases as q decreases in the first few steps, and then optimal cost remains at the global minimum. This phenomenon obeys Proposition 7. In this example, global optimal is reached at M=5 for all pairs of initial and final states as J7≠J5=J3 and J6≠J4=J2. So, the global minimum is achieved before we reach N-5=4 in this particular case. This also implies that with N free we can reach desired final state in 4 steps from given initial state.Observing closely toJq,  0≤q≤5, Jq-1 equals to Jq by exchanging the first and the last columns, and the second and the third columns, which is consistent with Proposition 8. Or we can exchange the first and the last rows, and the second and the third rows of Jq to obtain Jq-1. Jq1 and Jq2 are the same for q1,q2≤M=5 for q1-q2=0 (mod 2). This obeys Proposition 7.The optimal trajectories are generated fromPq(xq,xNd). For example, given x2=T, and the final state x9d=G, we want to generate the optimal trajectories.x3=P2(T,G)=A,G.Ifx3=A, x4=P3(A,G)=T;  if x3=G, x4=P3(A,G)=C.Ifx4=T, x5=P4(T,G)=A,G; if x4=C, x5=P4(C,G)=G.Ifx5=A, x6=P5(A,G)=T; if x5=G, x6=P5(G,G)=C.Ifx6=T, x7=P6(T,G)=A,G; if x6=C, x7=P6(C,G)=G.Ifx7=A, x8=P7(A,G)=T; if x7=G, x8=P7(G,G)=C.So, the optimal routes areT→A→T→A→T→A→T→αTGuTGG;T→A→T→A→T→αTGuTGG→C→G;T→A→T→αTGuTGG→C→G→C→G;T→αTGuTGG→C→G→C→G→C→G.Consequently, the optimal cost isJ2(T,G)=0.64=αTG.Optimal trajectories for other pairs of initial and final states can be obtained in the same manner.It takes less than 1 second to generate optimal trajectories for all pairs of initial and final states withN=9 on a regular desktop. Since we have already proven by the brute force method that the global optimal can be achieved with M≤6, the computation time can be further reduced by taking N=6 with all the results necessary for this example.


## 3.2. Codon-to-Codon, Deterministic Optimal Control Problem

For codon-to-codon deterministic mutations, we formulate our optimal control problem as(33)J0(x0)=min⁡uk,l1i,0≤k≤N,1≤l1≤l,1≤i≤3{∑k=0N-1∑l1=1l∑i=13αl1uk,l1i+d(xN,{xNd})},
subject to(34)xk+1=(-I+∑l1=1l∑i=13uk,l1iseieiT)xk,
with x0,xNd∈ℱ∖{0}3 given, xk∈ℱ∖{0}3, for all k,  0≤k≤N, and d(φ1,φ2),φ1,φ2∈ℱ∖{0}3 as defined in Section 2.3.If we take the same assumption on available mutagens as in Section3.1, then we can write update equation for optimal control policy explicitly as
(35)Jq(xq)=min⁡uq,l1i,1≤l1≤l,1≤i≤3×{αxq1ψ1+Jq+1([ψ1xq2¯xq3¯]),αxq2ψ2+Jq+1([xq1¯ψ2xq3¯]),αxq3ψ3+Jq+1([xq1¯xq2¯ψ3]),([xq1¯xq2¯ψ3]),ψ1,ψ2,ψ3∈{A,T,G,C}⟺F∖{0}},
with xqi∈{A,T,G,C}⇔ℱ∖{0},  1≤i≤3 denotes the ith element of xq∈ℱ∖{0}3, and xqi¯ denotes the complementary base of xqi.The optimal control sequences depend on the numerical values ofαl1s and d(φ1,φ2), φ1,φ2∈ℱ∖{0}3. Though we do not have real values of αl1s and d(φ1,φ2), we can always obtain simulation results to compare the result differences by assigning different numerical values to those parameters.Therefore, we use three different assignments ofαl1s and the same d(φ1,φ2) to generate our simulation results. Those three assignments of αl1s are χ,  5χ, and 0.5χ, respectively, with χ the same as assigned in Section 3.1. In every particular example, it takes approximately 2 seconds on a regular desktop to generate the optimal path table for all pairs of initial and final states with N=19, and the dynamic programming algorithm ensures that the optimal control for tail problems is generated at the same time.The graphical interpretation of three assignments are shown in Figures7, 8, and 9, respectively. The x-axis and y-axis denote xq and xNd, respectively. For a codon [ψ1ψ2ψ3]T,ψ1,ψ2, ψ3∈{A,T,G,C}⇔ℱ∖{0}, its index is calculated by
(36)42(ψ1-1)+4(ψ2-1)+ψ3,where ψi=1 if A, ψi=2 if G, ψi=3 if C, and ψi=4 if T, 1≤i≤3, for the simplicity of graphical interpretation. A codon has 43=64 combinations. Thus, there are 642 pairs of initial and final desired states, and there are 64×21 pairs of initial state and final desired set. The surface is generated by connecting 64×64 discrete points together. Jq is calculated following the same procedure as in base-to-base deterministic cases. The value of optimal cost can be read directly from graphical interpretation, and the optimal path can be generated from path matrix Pq, similar to base-to-base deterministic case. Both Jq and Pq,forallq,0≤q≤N, are of 64×64 dimension.Figure 7
Graphical representation ofJq,  q=0,1,2,3,15,16,17,18 for codon-to-codon deterministic mutations, with αl1=χ, N=19.Figure 8
Graphical representation ofJq,  q=0,1,2,3,15,16,17,18 for codon-to-codon deterministic mutations, with αl1=5χ, N=19.Figure 9
Graphical representation ofJq,  q=0,1,2,3,15,16,17,18 for codon-to-codon deterministic mutations, with αl1=0.5χ, N=19.From the graphical interpretation and Table10, we find that the value of q where the global minimum is reached at the first time decreases as αl1 decreases. And J0 is more similar to J18 with α=5χ than with α=χ or α=0.5χ. This implies that if d(φ1,φ2) are the deterministic term in our objective function, then the treatment plan is made to drive the final state as close to the desired set as possible; if the costs of applying mutagens is the deterministic term in the objective function, then the treatment plan tends to stay in the original state and applies less mutagens; if they are of equal weight, then the treatment plan deals with this tradeoff.Table 10
Simulation results with differentαl1 assignments and the first q where the global optimal is reached.
Nαl1M (first q when the  global minimum is reached)Global minimum190.5χq=12J12(x12)19χq=13J13(x13)195χq=15J15(x15)Moreover, no matter how the numerical values of final penalty and the costs of applying mutagens changes in our objective function as shown in (33), there is always a M, M≤N-18, JM(xM) is global minimum. Proposition 7 can be extended to codon-to-codon deterministic mutations as stated in Proposition 9.Proposition 9.
Given an optimal control problem with objective function in the form of (33), constraints in the form of (34), and all available chemical mutagens and their corresponding transfer pairs and costs as listed in Tables 6, 7, and 8, Jq(φ¯)≤Jq+1(φ),φ∈ℱ∖{0}3. If, in addition, the system is completely controllable, ∃M, s.t. JM(φ) is the global minimum and for all q≤M, Jq(φ)=JM(φ¯) if M-q≡1 (mod 2), and Jq(φ)=JM(φ) if M-q≡0 (mod2). In our example, M≥N-18.Proof.
We only prove thatM≥N-18 here since the rest is similar to the proof of Proposition 7.
The objective function in (33) can be written as the summation of three separate single-base mutation systems and the distance reference between final states and the final desired set, that is, (37)Jq(xq)=min⁡n1,n2,n3≥02N≤n1+n2+n3≤3N-1{Jn1(xn-1)(xq1,ψ1)︸optimalcostsofbase-to-basedeterministicoptimalcontrolproblemformedbythe1stbase+Jn2(xn-2)(xq2,ψ2)︸optimalcostsofbase-to-basedeterministicoptimalcontrolproblemformedby-the2ndbase+Jn3(xn-3)(xq3,ψ3)︸optimalcostsofbase-to-basedeterministicoptimalcontrolproblemformedbythe3rdbase+d([ψ1ψ2ψ3],{xNd})Jn1(xn-1)(xq1,ψ1)︸optimalcostsofbase-to-basedeterministicoptimalcontrolproblemformedbythe1stbase},
where N-q=(N-n1)+(N-n2)+(N-n3).
According to Proposition7, JN-6(xN-6) is guaranteed to be the global optimal for single-base mutations. Therefore, optimal costs corresponding to three single-base mutation systems, Jn1(xn-1)(xq1,ψ1), Jn2(xn-2)(xq2,ψ2), Jn3(xn-3)(xq3,ψ3) are guaranteed to reach their own global optimal at n1=n2=n3=N-6 with all possible combinations of ψ1,ψ2,ψ3∈{A,T,G,C}. Therefore, q=N-18 is a guaranteed global optimal.Indeed, Proposition9 is a quite loose condition. In real examples above, the M value where the first global optimal is reached at earlier stage as listed in Table 10.Graphically, the indices of complementary codonsφ1=[ψ1ψ2ψ3]T and φ2=[ψ1¯ψ2¯ψ3¯]T sum up to 65, that is,  (38)(16(ψ1-1)+4(ψ2-1)+ψ3)+(16(ψ1¯-1)+4(ψ2¯-1)+ψ3¯)=(16(ψ1-1)+4(ψ2-1)+ψ3)+(16((5-ψ1)-1)+4((5-ψ2)-1)+(5-ψ3))=65.Therefore,Jq and Jq-1, 1≤q≤M are symmetric about the plane x=32.5, Jq-2 and Jq, 2≤q≤M are the same, as shown in Figures 7, 8, and 9.However, Proposition8 cannot be extended to codon-to-codon deterministic case due to the redundancy of genetic codes, that is, the set of codons translated to the same amino acid varies from one amino acid to another as shown in Table 5. The simulation results show that the costs, a pair of complementary codons, to two final desired set generated by a pair of complementary final desired codons are different, that is, Jq(xq,{xNd})≠Jq(xq¯,{xNd¯}), in general, for any q. Graphically, the optimal cost profile Jq is not symmetric about the plane y=32.5 for Jq-1, 1≤q≤M. Therefore, the doctors need to pick the strand with lower cost to make the treatment plan. This also implies that in large-scale cases, for instance, a gene containing hundreds of nucleotide bases, the doctors should make the treatment plan based on the strand the total cost of which is lower than the other.


## 3.3. Codon-to-Codon, Stochastic Optimal Control Problem

The optimal control problem of codon-to-codon stochastic mutations can be written as(39)J0(x0)=min⁡uk,l1i,0≤k≤N-11≤l1≤l,1≤i≤3{+Ehk,l1i,0≤k≤N-11≤l1≤l,1≤i≤3[d(xN,{xNd})][d(xN,{xNd})]∑k=0N-1∑l1=1l∑i=13αl1uk,l1i+Ehk,l1i,0≤k≤N-11≤l1≤l,1≤i≤3[d(xN,{xNd})]},
subject to(40)xk+1=-Ixk+∑l1=1l∑i=13uk,l1ihk,l1iseieiTxk,
with x0,xNd∈ℱ∖{0}3 given, xk∈ℱ∖{0}3.The major difference between deterministic and stochastic systems is that we impose the random binary vector,hk,l1i, in our system equation (40). We denote the probability associated with hk,l1i,j to be  pl1,ψ1ψ2(h) with ψ1,ψ2∈{A,T,G,C}. The equivalence relationship between j and ψ1ψ2 can be found in Table 3.Again, we assume that we havel1=12 kinds of mutagens, each corresponding to one major transfer pattern, associated with probability assignments, as listed in Table 11.Table 11
12 kinds of mutagens, each corresponding to major transfer patterns, and the probabilities of different mutagens on different transfer patterns.
Index (l1)FromToMajor transfer patternAGCT1Ap1,AA(h)p1,AG(h)p1,AC(h)p1,AT(h)A→A2Ap2,AA(h)p2,AG(h)p2,AC(h)p2,AT(h)A→G3Ap3,AA(h)p3,AG(h)p3,AC(h)p3,AT(h)A→C4Gp4,GA(h)p4,GG(h)p4,GC(h)p4,GT(h)G→A5Gp5,GA(h)p5,GG(h)p5,GC(h)p5,GT(h)G→G6Gp6,GA(h)p6,GG(h)p6,GC(h)p6,GT(h)G→T7Cp7,CA(h)p7,CG(h)p7,CC(h)p7,CT(h)C→A8Cp8,CA(h)p8,CG(h)p8,CC(h)p8,CT(h)C→C9Cp9,CA(h)p9,CG(h)p9,CC(h)p9,CT(h)C→T10Tp10,TA(h)p10,TG(h)p10,TC(h)p10,TT(h)T→G11Tp11,TA(h)p11,TG(h)p11,TC(h)p11,TT(h)T→C12Tp12,TA(h)p12,TG(h)p12,TC(h)p12,TT(h)T→TThen, we can write updated formula for optimal control policy explicitly as(41)Jq(xq)=min⁡uq,l1i,1≤l1≤l,1≤i≤3{αxq1ψ1+Ehq,l1(xq1ψ1)1[Jq+1([⋅xq2¯xq3¯])],αxq2ψ2+Ehq,l1(xq2ψ2)2[Jq+1([xq1¯⋅xq3¯])],αxq3ψ3+Ehq,l1(xq3ψ3)3[Jq+1([xq1¯xq2¯⋅])],ψ1,ψ2,ψ3∈{A,T,G,C}⟺F∖{0}xq1¯xq2¯⋅}
where(42)Ehq,l1(xq1ψ1)1[Jq+1([⋅xq2¯xq3¯])]=pl1(xq1ψ1),xq1A(h)[Jq+1([Axq2¯xq3¯])]+pl1(xq1ψ1),xq1G(h)[Jq+1([Gxq2¯xq3¯])]+pl1(xq1ψ1),xq1C(h)[Jq+1([Cxq2¯xq3¯])]+pl1(xq1ψ1),xq1T(h)[Jq+1([Txq2¯xq3¯])],wherexqi∈{A,T,G,C}⇔ℱ∖{0},  0≤q≤N-1, 1≤i≤3 denotes the ith element of xq∈ℱ∖{0}3, xqi¯ denotes the complementary base of xqi, and l1:ψ1ψ2∈{A,T,G,C}×{A,T,G,C}→{integersfrom1to12}, the mapping from major transfer pattern ψ1→ψ2 to mutagen index, as shown in Table 11. The mathematical expression of 𝔼hq,l1(xq2ψ2)2[Jq+1([xq1¯·xq3¯]T)] and 𝔼hq,l1(xq3ψ3)3[Jq+1([xq1¯xq2¯·]T)] is similar to 𝔼hq,l1(xq1ψ1)1[Jq+1([·xq2¯xq3¯]T)] as shown above.In order to run the simulation, we assign numerical values to probabilities in Table11, as illustrated in Table 12.Table 12
Sample probabilities with respect to different mutagens and different transfer patterns.
Index (l1)FromToAGCT1A0.900.050.030.022A0.110.580.210.103A0.140.160.420.284G0.850.070.030.055G0.020.020.920.046G0.100.090.220.597C0.790.130.040.048C0.010.020.870.109C0.040.120.090.7510T0.130.760.050.0611T0.070.030.620.2812T0.080.040.250.63As in Section3.2, we use three different assignments for αl1s, χ,5χ, and 0.5χ, respectively. The optimal cost profile Jq with selected q values, for every pair of (xq,xNd), is graphically interpreted in Figures 10, 11, and 12, respectively, with N=29, with computation time of approximately 7 seconds on a regular desktop.Figure 10
Graphical representation ofJq(xq),  q=0,1,2,3,15,19,24,28 for codon-to-codon stochastic mutations, with αl1=χ, probability assignment as in Table 12, N=29.Figure 11
Graphical representation ofJq(xq),  q=0,1,2,3,15,19,24,28 for codon-to-codon stochastic mutations, with αl1=5χ, N=29.Figure 12
Graphical representation ofJq(xq),  q=0,1,2,3,15,19,24,28 for codon-to-codon stochastic mutations, with αl1=0.5χ, N=29.The simulation results are similar to those in Section3.2. The profile of J0 is more similar to J29 when αl1s are assigned 5χ than χ or 0.5χ. This implies in codon-to-codon stochastic mutations; the optimal control sequence behaves as codon-to-codon deterministic cases, that is, the system tends to getting as close as possible to the final desired set if αl1s are much smaller than d(φ1,φ2), and the system tends to remain in the same state with minor mutations when αl1s are relatively larger than d(φ1,φ2),φ1,φ2∈ℱ∖{0}3. And Jq(φ¯)≤Jq+1(φ),  φ∈ℱ∖{0}3 is still valid in codon-to-codon stochastic case.However, for stochastic cases, we cannot reach a global minimum because of the randomness caused by mutagens. Since, in usual cases, there exists no stationary global minimum, we need to define error toleranceϵ, that is, if |Jq(ψ)-Jq-1(ψ¯)|≤ϵ with the same {xNd}, then we can stop at Jq(xq). Otherwise, we need to proceed to calculate Jq-2(ψ). The value of ϵ is decided based on the doctors, experience. Obviously, the smaller ϵ is, the better treatment plan is. However, we can still observe Figures 10, 11, and 12 to conclude that J0 and J2 are almost of the same shape in all three different parameter assignments. Higher dimensional optimal control problems, gene-to-gene stochastic mutations, can be solved as a series of cascade codon-to-codon stochastic problems.


## 4. Conclusion

In this paper, we present a mathematical model to deal with mutations in the process of DNA replication in the view of control systems. Different from the existing models, our model is constructed directly from the basic biological theories, the central dogma in molecular biology, and the complementary base pair for DNA molecules with double helix structure. It precisely describes how the induced mutations affect the targeted DNA segments at molecular level. It provides instrumental information of molecule interactions in gene mutation for biologists and doctors to gain a better understanding of cellular and tissue level systems’ behavior. Our model is adaptive to point and multisite, deterministic and stochastic mutations. Though we emphasize that we target at induced mutations during the process of DNA replication in our work, this model can be extended to other biological processes at molecular level, such as transcription process and DNA brokage.In our optimal control problem, the objective function includes two factors: the risk/cost of applying mutagens and the off-trajectory penalty. Under optimal control policy, the summation of those two factors are minimized, by dynamic programming, to propose a low-risk treatment plan. We define the distance reference following the chemical and physical properties of amino acids, representing the penalty. Our objective is to drive the system from given initial state to the final desired set generated by the final desired state at the lowest cost. We define the final desired set since redundancy in genetic codes gives us additional options of final desired state to further lower the cost. Dynamic programming algorithm ensures the optimality of the solution. We also discuss three different small-scale system, and show the simulation results of examples. The optimal control problems of base-to-base deterministic mutations and codon-to-codon deterministic mutations are of theoretical importance. As shown in the propositions, the global optimal can be reached within finite steps. If the step limit is larger than the number of steps that global optimal can achieve, then we have some flexibility in our treatment plan. In addition, there exist multiple optimal paths with the same total cost, given the initial state and the final desired set. The optimal control problem of codon-to-codon stochastic mutations is of practical importance, since codon is the basic component forming long DNA sequences. The step limitN is decided by doctors according to patients’ conditions, and the treatment plan is made according to the initial state, the final desired set, and the step limit. Since the doctors constantly take measurement to see the result of treatments at current stage, the treatment plan is updated accordingly. Solving codon-to-codon stochastic optimal control problem is a key step to realize the optimal control to gene-to-gene stochastic mutations in the real world.Our work contributes to several aspects in systems biology. The optimal control sequences generated by dynamic programming make it possible for biologists and doctors to mutate certain sections of genes on purpose in laboratory, at a relative low cost and low risk, which is an essential step to identify the functional units, to examine the interactions among different segments, and to find healthy, harmful, and lethal nucleotide sequences. All those results are beneficial in gene network construction. In addition, the fundamental details of gene mutations at molecular level help biologists to elaborate on the biological theories at the cellular and tissue levels, such as the theory of evolution. Moreover, by our method, biologists can distinguish the harmful and beneficial mutations and induce beneficial mutations during the evolution process in a proper way, which greatly helps to save rare species in danger. Furthermore, our solution to the optimal control problem proposed provides a new medical intervention to genetic diseases. Compared to the existing gene therapy, treatments by mutagens are safer by avoiding the side effect of virus infection. Lastly, our work also contributes to the construction of DNA computers. Calculation errors, the mispairings in the process of two single-stranded DNA segments, can be corrected at lowest cost by applying a correct mutagen sequence.Further work can be done by extending codon-to-codon stochastic optimal control problem to gene-to-gene stochastic mutations. The distance reference between DNA segment with equal length can be defined as a weighted sum of the distance references between codons. Since certain combinations of amino acids are lethal, those high-risk states should be avoided. This goal can be achieved by either defining a preset trajectory or eliminating high-risk sequences in the state space. Another possibility is to examine system’s behavior under noisy measurements. Under this situation, the spontaneous mutations can be modeled as an additional random factor in our state updating equations, and another random noises should be added to the output equation.



---

*Source: 743172-2012-01-22.xml*