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.

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.

Systems biology is a cross-cutting research area connecting control engineering, biology, and medical science, as shown in Figure 2. 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.

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, 1518], 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 Section 2.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 (𝐴), thymine (𝑇), guanine (𝐺), and cytosine (𝐶), 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]. 𝐴 always pairs with 𝑇, 𝐺 always pairs with 𝐶, 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 (𝐴𝐺 or 𝐶𝑇) and transversions (𝐴/𝐺𝐶/𝑇). 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 Section 2.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

Figure 3 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.

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 Figure 4. 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.

2.2.1. Base-to-Base Deterministic Mutations

Denote the targeted DNA segment with 𝑛 nucleotide bases at 𝑘th step by a column vector 𝑥𝑘, as shown in Figure 5. 𝑥𝑖𝑘 is the 𝑖th element of 𝑥𝑘. Let 𝑃 be the transfer matrix from 𝑥𝑘 to 𝑥𝑘+1, for all 𝑘, 𝑘+{0}, without mutation. Then, the perfect DNA replication process can be expressed as𝑥𝑘+1=𝑃𝑥𝑘.(1)

Proposition 1. 𝑃=𝐼.

Proof. As no mutation occurs, 𝑥𝑘+1 is completely complementary to 𝑥𝑘 by Watson-Crick base pairing rule, and 𝑥𝑘+2 is completely complementary to 𝑥𝑘+1. Therefore, 𝑥𝑘+2 is exactly the same as 𝑥𝑘. thus, 𝑥𝑘+2=𝑃𝑥𝑘+1=𝑃2𝑥𝑘𝑃2=𝐼.(2) Since every base mutates independently, every element of 𝑥𝑖𝑘+1 only depends on the corresponding element of 𝑥𝑖𝑘, thus 𝑃 is diagonal. In addition, 𝑥𝑘+1𝑥𝑘, we conclude 𝑃=𝐼.

Based on Proposition 1, we assign values to nucleotide bases set {𝐴,𝐺,𝐶,𝑇,𝑂}, where 𝑂 is an artificial nonsense base. Define an equivalence relationship between {𝐴,𝐺,𝐶,𝑇,𝑂} and {1,2,2,1,0}, that is, {𝐴,𝐺,𝐶,𝑇,𝑂}{1,2,2,1,0}, with𝑥𝑖𝑘=1,if𝐴,2,if𝐺,2,if𝐶,1,if𝑇,0,if𝑂.(3)

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 Tables 1 and 2, we check if the set {1,2,2,1,0} satisfies the definition of field.

Closed under Addition and Multiplication
Satisfied obviously from Tables 1 and 2.

Associativity of Addition and Multiplication
Implicitly satisfied by integer addition and multiplication.

Commutativity of Addition and Multiplication
Satisfied as Tables 1 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: 11, 22, 00.
Multiplicative inverses pair: 11, 22, 11.

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 𝑥𝑘𝑛 is the state vector representing a DNA segment with 𝑛 nucleotide bases, where 𝑛 is the set of -valued vectors of dimension 𝑛.

We start with the simplest form of mutations, point mutation. Suppose there is a point mutation, we write it mathematically as𝑥𝑘+1=(𝐼+Δ𝑠)𝑥𝑘+Δ𝑤,(4) where 𝑥𝑘+1,𝑥𝑘, and 𝐼 reduces to −1 as only one base is involved. The corresponding values of Δ𝑠 and Δ𝑤, obtained by reverse engineering with all possible pairs of 𝑥𝑘 and 𝑥𝑘+1, are listed in Table 3.

Here, Δ𝑠 represents the mutation from four normal nucleotide bases, and Δ𝑤 corresponds to mutation from nonsense base, that is, Δ𝑤0 only if 𝑥𝑘=0.

Rewriting (4) by collecting all values of Δ𝑠 and Δ𝑤 in Table 3, we get𝑥𝑘+1=𝐼+4𝑗=0𝑢𝑗𝑘𝑠𝑗𝑥𝑘+4𝑗=0𝑐𝑗𝑘𝑤𝑗(5a)=𝐼+𝑢𝑘𝑠𝑥𝑘+𝑐𝑘𝑤,(5b)where {𝑠0,𝑠1,𝑠2,𝑠3,𝑠4}={𝑤0,𝑤1,𝑤2,𝑤3,𝑤4}={0,1,2,2,1}, 𝑢𝑗𝑘,𝑐𝑗𝑘{0,1}, representing the on/off controls, 𝑢𝑘=[𝑢0𝑘𝑢1𝑘𝑢2𝑘𝑢3𝑘𝑢4𝑘], 𝑐𝑘=[𝑐0𝑘𝑐1𝑘𝑐2𝑘𝑐3𝑘𝑐4𝑘], and 𝑠=𝑤=[01221]𝑇. In (5a), 𝑠𝑗 and 𝑤𝑗 are constants for all 𝑘 and 𝑗. 𝑢𝑗𝑘 and 𝑐𝑗𝑘, the inputs of the system, are the on/off controls for chemical mutagens or radiation. Clearly, 4𝑗=0𝑐𝑗𝑘=1 only if 𝑥𝑘=0. Equation (5b) is a simplified version of (5a) as we put 𝑢𝑗𝑘,𝑠𝑗,𝑐𝑗𝑘,𝑤𝑗 into vector form 𝑢𝑘,𝑠,𝑐𝑘,𝑤. 𝑠 and 𝑤 serve as vector basis for base-to-base deterministic model. 𝑢𝑘 and 𝑐𝑘 are now multi-input controls; each of them contains 5 on/off controls, corresponding to all possible transfer patterns. For a particular 𝑘, at most one of 𝑢𝑗𝑘s and 𝑐𝑗𝑘s 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 always 11 transfer when mutation occurs, that is, one nucleotide base can only transfer to another one, therefore(i)if 𝑥𝑘=0 and 𝑐𝑘=0, or 𝑐𝑘=[10000], then 𝑥𝑘+1=0,(ii)if 𝑥𝑘0, then 𝑐𝑘=0 and 𝑢𝑘 is either 0 or a unit row vector,(iii)if 𝑥𝑘=0, then 𝑢𝑘=0 and 𝑐𝑘 is either 0 or a unit row vector,(iv)𝑢𝑘+𝑐𝑘 is either 0 or a unit row vector, for all 𝑘+{0}.

Now, suppose for some reason we need to take an addition, measurement in the middle of the cell cycle, after the completion of the 𝑘th duplication and before the start of the (𝑘+1)th. We name this kind of measurement an intermediate state, and denote by it 𝑥𝑘. Then, we have𝑥𝑘+1=𝐼+Δ𝑠𝑥𝑘+Δ𝑤,(6) where the values of Δ𝑠 and Δ𝑤, listed in Table 4, are obtained in the same way as getting Δ𝑠 and Δ𝑤 in Table 3.

Comparing Tables 3 and 4, we find the collection of Δ𝑠 and Δ𝑠, Δ𝑤 and Δ𝑤, form the same set, respectively. Thus, we continue using 𝑠 and 𝑤 when rewriting (6) in the form of (5a) and (5b), that is, 𝑥𝑘+1=𝐼+𝑣𝑘𝑠𝑥𝑘+𝑐𝑘𝑤,(7) where 𝑣𝑘,𝑐𝑘 are the counterparts of 𝑢𝑘,𝑐𝑘, respectively, and 𝑠,𝑤 are the same as in (5b).

Similar to Proposition 3, we get Proposition 4.

Proposition 4. 𝑣𝑘 and 𝑐𝑘 in (7) need to satisfy the following conditions.(i)If 𝑥𝑘=0 and 𝑐𝑘=0, or 𝑐𝑘=[10000], then 𝑥𝑘+1=0.(ii)If 𝑥𝑘0, then 𝑐𝑘=0 and 𝑣𝑘 is either 0 or a unit row vector.(iii)If 𝑥𝑘=0, then 𝑣𝑘=0 and 𝑐𝑘 is either 0 or a unit row vector.(iv)𝑣𝑘+𝑐𝑘 is either 0 or a unit row vector, for all 𝑘+{0}.

Now take, both chemical mutagens and radiative rays under our consideration and apply them in the order as shown in Figure 4. Then, we can express our system equation as𝑥𝑘=𝐼+𝑢𝑘𝑠mutationscausedbychemicalmutagensfromnormalbases𝑥𝑘+𝑐𝑘𝑤mutationscausedbychemicalmutagensfrom𝑂,(8a)𝑥𝑘+1=𝐼+𝑣𝑘𝑠mutationscausedbyradiativeraysfromnormalbases𝑥𝑘+𝑐𝑘𝑤mutationscausedbyradiativeraysfrom𝑂,(8b)𝑦𝑘=𝑥𝑘,(8c)where 𝑢𝑘 and 𝑣𝑘 are the inputs of the system and 𝑦𝑘 is the measurement. Obviously, (8a) is modified from (5b) and (8b) from (7).

The two-step mutation and the intermediate state 𝑥𝑘avoid the case 𝑥𝑘 is changed to different bases by radiation and chemical mutagens simultaneously, which causes confusion. Substituting (8a) and (8b), we get𝑥𝑘+1=𝐼+𝑣𝑘𝑠𝐼+𝑢𝑘𝑠𝑥𝑘+𝐼+𝑣𝑘𝑠𝑐𝑘𝑤+𝑐𝑘𝑤,(9a)𝑦𝑘=𝑥𝑘.(9b)Obviously, Proposition 3 still holds for 𝑢𝑘 and 𝑐𝑘, and Proposition 4 holds for 𝑣𝑘 and 𝑐𝑘 for (9a).

For point mutations, we have 20 on/off controls in total for every step 𝑘, 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 length 𝑛, then 𝑥𝑘𝑛. Since a coding DNA segment usually contains integer number of codons, which is made of three consecutive bases, 𝑛 is a multiple of 3. Let 𝑥𝑖𝑘 denote the 𝑖th component of 𝑥𝑘. 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 𝑥𝑘=𝐼+𝑛𝑖=1𝑢𝑖𝑘𝑆𝑖𝑘mutationscausedbychemicalmutagensfromnormalbases𝑥𝑘+𝑖𝒪𝑘𝑐𝑖𝑘𝑊𝑖𝑘mutationscausedbychemicalmutagensfrom𝑂,𝑥𝑘+1=𝐼+𝑛𝑖=1𝑣𝑖𝑘𝑄𝑖𝑘mutationscausedbyradiativeraysfromnormalbases𝑥𝑘+𝑖𝒪𝑘𝑏𝑖𝑘𝑅𝑖𝑘mutationscausedbyradiativeraysfrom𝑂,𝑦𝑘=𝑥𝑘,(10) where 𝑢𝑖𝑘,𝑣𝑖𝑘,𝑐𝑖𝑘,𝑏𝑖𝑘 are on/off controls of the 𝑖th element,𝑆𝑖𝑘,𝑄𝑖𝑘 are 𝑛×𝑛 square matrices corresponding to the mutations between normal bases or from normal bases by chemicals and radiation, respectively, 𝑊𝑖𝑘,𝑅𝑖𝑘 are 𝑛-dimensional column vectors representing mutations from nonsense bases by chemicals and radiation, respectively, and 𝒪𝑘={𝑖𝑥𝑖𝑘=0,1𝑖𝑛}, 𝒪𝑘={𝑖𝑥𝑖𝑘=0,1𝑖𝑛}.

𝑆𝑖𝑘 and 𝑄𝑖𝑘 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 𝑆𝑖𝑘 and 𝑄𝑖𝑘, respectively. The last rows of Tables 3 and 4 are assigned to 𝑊𝑖𝑘 and 𝑅𝑖𝑘, 𝑛-dimensional vectors, at nonsense base’s spots for 𝑥𝑘.

Define 𝒮={𝑠𝑗𝑒𝑖𝑒𝑇𝑖,𝑖,𝑗,0𝑗4,1𝑖𝑛}, a collection of 𝑛×𝑛 matrices, where 𝑠𝑗 is the same as in (5a) and (5b), 𝑒𝑖 is the unit column vector of length 𝑛 with 𝑖th component equal to 1 and all other components equal to 0, and 𝑒𝑖𝑒𝑇𝑖 is the square matrix with only the 𝑖th element on the diagonal equals to 1, and 0 otherwise. Then, 𝑆𝑖𝑘,𝑄𝑖𝑘 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 𝑢(𝑖,𝑗)𝑘 and 𝑣(𝑖,𝑗)𝑘, respectively.

Similarly, define 𝒲={𝑤𝑗𝑒𝑖,𝑖,𝑗,0𝑗4,1𝑖𝑛}, where 𝑤𝑗 is the same as (5a) and (5b). 𝑊𝑖𝑘,𝑅𝑖𝑘 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 𝑐(𝑖,𝑗)𝑘 and 𝑐(𝑖,𝑗)𝑘, respectively.

Therefore, instead of using step-varying 𝑆𝑖𝑘,𝑆𝑖𝑘,𝑊𝑖𝑘,𝑅𝑖𝑘, we find matrix basis for those four square matrices to make the controls to be the only variables depending on 𝑘, as we did for single-base cases. Then, we can, write (10) as𝑥𝑘=𝐼+𝑛𝑖=14𝑗=0𝑢(𝑖,𝑗)𝑘𝑠𝑗𝑒𝑖𝑒𝑇𝑖𝑥𝑘+𝑖𝒪𝑘4𝑗=0𝑐(𝑖,𝑗)𝑘𝑤𝑗𝑒𝑖,(11a)𝑥𝑘+1=𝐼+𝑛𝑖=14𝑗=0𝑣(𝑖,𝑗)𝑘𝑠𝑗𝑒𝑖𝑒𝑇𝑖𝑥𝑘+𝑖𝒪𝑘4𝑗=0𝑐(𝑖,𝑗)𝑘𝑤𝑗𝑒𝑖,(11b)𝑦𝑘=𝑥𝑘,(11c)where 𝑢(𝑖,𝑗)𝑘,𝑣(𝑖,𝑗)𝑘,𝑐(𝑖,𝑗)𝑘,𝑐(𝑖,𝑗)𝑘{0,1}.

As shown in (11a), (11b), and (11c), multisites mutations contain 20𝑛 controls in total for every step 𝑘, where 𝑛 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 view 𝑢𝑘,𝑣𝑘,𝑐𝑘,𝑐𝑘 as binary matrices of dimension 𝑛×5, and 𝑢(𝑖,𝑗)𝑘, 𝑣(𝑖,𝑗)𝑘𝑐(𝑖,𝑗)𝑘, 𝑐(𝑖,𝑗)𝑘are the corresponding element of 𝑖th row and 𝑗th column. Use 𝑢𝑖𝑘,𝑣𝑖𝑘,𝑐𝑖𝑘,𝑏𝑖𝑘, binary row vectors of dimension 5, to denote the 𝑖th row of 𝑢𝑘,𝑣𝑘,𝑐𝑘,𝑐𝑘, respectively. Again, 𝑠=𝑤=[01221]𝑇.

Combining (11a) and (11b), and writing control variables in vector forms, we get𝑥𝑘+1=𝐼+𝑛𝑖=1𝑣𝑖𝑘𝑠𝑒𝑖𝑒𝑇𝑖𝐼+𝑛𝑖=1𝑢𝑖𝑘𝑠𝑒𝑖𝑒𝑇𝑖𝑥𝑘+𝐼+𝑛𝑖=1𝑣𝑖𝑘𝑠𝑒𝑖𝑒𝑇𝑖𝑖𝒪𝑘𝑐𝑖𝑘𝑤𝑒𝑖+𝑖𝒪𝑘𝑏𝑖𝑘𝑤𝑒𝑖,𝑦𝑘=𝑥𝑘.(12)

Proposition 5. For large-scale deterministic system, 𝑢𝑘,𝑣𝑘,𝑐𝑘,𝑐𝑘 satisfy conditions below.(i)If 𝑒𝑇𝑖𝑥𝑘=0, then 𝑖𝒪𝑘.(ii)If 𝑒𝑇𝑖𝑥𝑘=0, 𝑐i𝑘=0 or 𝑐𝑖𝑘=[10000], then 𝑖𝒪𝑘.(iii)For all 𝑖𝒪𝑘, 𝑢𝑖𝑘 is either 0 or a row unit vector and 𝑐𝑖𝑘=0.(iv)For all 𝑖𝒪𝑘, 𝑐𝑖𝑘 is either 0 a row unit vector and 𝑢𝑖𝑘=0.(v)For all 𝑖𝒪𝑘, 𝑣𝑖𝑘 is either 0 or a row unit vector and 𝑏𝑖𝑘=0.(vi)For all 𝑖𝒪𝑘, 𝑏𝑖𝑘 is either 0 or a row unit vector and 𝑣𝑖𝑘=0.(vii)For all 𝑖,𝑘, 1𝑖𝑛, 𝑘+{0}, 𝑢𝑖𝑘+𝑐𝑖𝑘 is either 0 or a unit row vector and 𝑣𝑖𝑘+𝑏𝑖𝑘 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 𝑥𝑘(1)𝑥𝑘(2)=𝐼𝑚+𝑚𝑖=1𝑢𝑖𝑘𝑠𝑒𝑖𝑒𝑇𝑖00𝐼𝑛𝑚+𝑛𝑖=𝑚+1𝑢𝑖𝑘𝑠𝑒𝑖𝑒𝑇𝑖×𝑥𝑘(1)𝑥𝑘(2)+𝑖𝒪𝑘,1𝑖𝑚𝑐𝑖𝑘𝑤𝑒𝑖𝑖𝒪k,(𝑚+1)𝑖𝑛𝑐𝑖𝑘𝑤𝑒𝑖,𝑥𝑘+1(1)=𝐼𝑚+𝑚𝑖=1𝑣𝑖𝑘𝑠𝑒𝑖𝑒𝑇𝑖𝑥𝑘(1)+𝑖𝒪𝑘,1𝑖𝑚𝑏𝑖𝑘𝑤𝑒𝑖,𝑥𝑘+1(2)=𝐼𝑛𝑚+𝑛𝑖=𝑚+1𝑣𝑖𝑘𝑠𝑒𝑖𝑒𝑇𝑖𝑥𝑘(2)+𝑖𝒪𝑘,(𝑚+1)𝑖𝑛𝑏𝑖𝑘𝑤𝑒𝑖.(13)

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, (𝑖,𝑗)𝑘,𝑙1,𝑟(𝑖,𝑗)𝑘,𝑙2,(𝑖,𝑗)𝑘,𝑙3,𝑟(𝑖,𝑗)𝑘,𝑙4{0,1}, associated with probability 𝑝()𝑙1,𝑗, 𝑝(𝑟)𝑙2,𝑗,𝑝()𝑙3,𝑗,𝑝(𝑟)𝑙4,𝑗, for all 𝑖,𝑘, 1𝑖𝑛, 𝑘+{0}, respectively, where 𝑘 is the step index, 𝑙1,𝑙2 are indices for chemical mutagens inducing mutation from normal bases and from 𝑂, respectively, 𝑙3,𝑙4 are indices for radiation inducing mutation from normal bases and from 𝑂, respectively, 𝑖 is the index of DNA segment, and the value of 𝑗 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, 4𝑗=0𝑝()𝑙1,𝑗=1,𝑙1,1𝑙1𝑙,4𝑗=0𝑝(𝑟)𝑙2,𝑗=1,𝑙2,1𝑙2𝑚,4𝑗=0𝑝()𝑙3,𝑗=1,𝑙3,1𝑙3𝑙,4𝑗=0𝑝(𝑟)𝑙4,𝑗=1,𝑙4,1𝑙4𝑚.(14) The controls are 𝑢𝑖𝑘,𝑙1,𝑐𝑖𝑘,𝑙2,𝑣𝑖𝑘,𝑙3,𝑏𝑖𝑘,l4{0,1}, with the fact that 1 representing mutagen with corresponding index is applied at 𝑖th spot of DNA segment at 𝑘th generation, and 0 representing mutagen with corresponding index is not applied at spot 𝑖 at 𝑘th step, similar to Sections 2.2.1 and 2.2.2. The mutagen indices 𝑙1,𝑙2,𝑙3,𝑙4 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 (𝑖,𝑗)𝑘,𝑙1,𝑟(𝑖,𝑗)𝑘,𝑙2,(𝑖,𝑗)𝑘,𝑙3,𝑟(𝑖,𝑗)𝑘,𝑙4, given the values of 𝑢𝑖𝑘,𝑙1,𝑐𝑖𝑘,𝑙2,𝑣𝑖𝑘,𝑙3,𝑏𝑖𝑘,𝑙4, and 𝑥𝑘.

Suppose we have (𝑙+𝑚) kinds of chemical mutagens available, with 𝑙 kinds to induce mutations from normal bases and 𝑚 kinds to induce mutations from 𝑂. And we have (𝑙+𝑚) kinds of radiation available, with 𝑙 kinds to induce mutations from normal bases and 𝑚 kinds to induce mutations from 𝑂. Therefore, we have total (𝑙+𝑚+𝑙+𝑚) controls for each spot 𝑖 at each step 𝑘. We can write our system equation as𝑥𝑘=𝐼+𝑙𝑙1=1𝑛𝑖=1𝑢𝑖𝑘,𝑙14𝑗=0(𝑖,𝑗)𝑘,𝑙1𝑠𝑗𝑒𝑖𝑒𝑇𝑖mutationscausedbychemicalmutagensfromnormalbases𝑥𝑘+𝑚𝑙2=1𝑖𝒪𝑘𝑐𝑖𝑘,𝑙24𝑗=0𝑟(𝑖,𝑗)𝑘,𝑙2𝑤𝑗𝑒𝑖mutationscausedbychemicalmutagensfrom𝑂,(15a)𝑥𝑘+1=𝐼+𝑙𝑙3=1𝑛𝑖=1𝑣𝑖𝑘,𝑙34𝑗=0(𝑖,𝑗)𝑘,𝑙3𝑠𝑗𝑒𝑖𝑒𝑇𝑖mutationscausedbyradiativeraysfromnormalbases𝑥𝑘+𝑚𝑙4=1𝑖𝒪𝑘𝑏𝑖𝑘,𝑙44𝑗=0𝑟(𝑖,𝑗)𝑘,𝑙4𝑤𝑗𝑒𝑖mutationscausedbyradiativeraysfrom𝑂,(15b)𝑦𝑘=𝑥𝑘.(15c)Again, we define (𝑖,𝑗)𝑘,𝑙1,𝑟(𝑖,𝑗)𝑘,𝑙2,(𝑖,𝑗)𝑘,𝑙3,𝑟(𝑖,𝑗)𝑘,𝑙4 the elements at 𝑖th row and 𝑗th column of 𝑛×5 binary matrices 𝑘,𝑙1,𝑟𝑘,𝑙2,𝑘,𝑙3,𝑟𝑘,𝑙4, respectively. 𝑖𝑘,𝑙1,𝑟𝑖𝑘,𝑙2,𝑖𝑘,𝑙3, 𝑟𝑖𝑘,𝑙4, and binary row vectors of dimension 5 denote the 𝑖th row of 𝑘,𝑙1,𝑟𝑘,𝑙2,𝑘,𝑙3,𝑟𝑘,𝑙4, respectively. Then, we can simplify (15a), (15b), and (15c) and combine (15a) and (15b) as𝑥𝑘+1=𝐼+𝑙𝑙3=1𝑛𝑖=1𝑣𝑖𝑘,𝑙3𝑖𝑘,𝑙3𝑠𝑒𝑖𝑒𝑇𝑖×𝐼+𝑙𝑙1=1𝑛𝑖=1𝑢𝑖𝑘,𝑙1𝑖𝑘,𝑙1𝑠𝑒𝑖𝑒𝑇𝑖𝑥𝑘+𝐼+𝑙𝑙3=1𝑛𝑖=1𝑣𝑖𝑘,𝑙3𝑖𝑘,𝑙3𝑠𝑒𝑖𝑒𝑇𝑖𝑚𝑙2=1𝑖𝒪𝑘𝑐𝑖𝑘,𝑙2𝑟𝑖𝑘,𝑙2𝑤𝑒𝑖+𝑚𝑙4=1𝑖𝒪𝑘𝑏𝑖𝑘,𝑙4𝑟𝑖𝑘,𝑙4𝑤𝑒𝑖,𝑦𝑘=𝑥𝑘.(16)

Proposition 6. For large-scale stochastic system, 𝑢𝑖𝑘,𝑙1,𝑖𝑘,𝑙1,𝑐𝑖𝑘,𝑙2,𝑟𝑖𝑘,𝑙2,𝑣𝑖𝑘,𝑙3,𝑖𝑘,𝑙3,𝑏𝑖𝑘,𝑙4,𝑟𝑖𝑘,𝑙4 follow the rules below.(i)If 𝑒𝑇𝑖𝑥𝑘=0, then 𝑖𝒪𝑘.(ii)If 𝑒𝑇𝑖𝑥𝑘=0 and 𝑚𝑙2=1𝑐𝑖𝑘,𝑙2=0, then 𝑖𝒪𝑘.(iii)If 𝑒𝑇𝑖𝑥𝑘=0, 𝑚𝑙2=1𝑐𝑖𝑘,𝑙2=1 and 𝑟𝑖𝑘,𝑙2=[10000], then 𝑖𝒪𝑘.(iv)For all 𝑖,𝑘,𝑙1, 1𝑖𝑛, 𝑘+{0},1𝑙1𝑙, if 𝑢𝑖𝑘,𝑙1=1, then 𝑖𝑘,𝑙1 is a unit row vector.(v)For all 𝑖,𝑘,𝑙2, 1𝑖𝑛, 𝑘+{0}, 1𝑙2𝑚, if 𝑐𝑖𝑘,𝑙2=1, then 𝑟𝑖𝑘,𝑙2 is a unit row vector.(vi)For all 𝑖,𝑘,𝑙3, 1𝑖𝑛, 𝑘+{0}, 1𝑙3𝑙, if 𝑣𝑖𝑘,𝑙3=1, then 𝑖𝑘,𝑙3 is a unit row vector.(vii)For all 𝑖,𝑘,𝑙4, 1𝑖𝑛, 𝑘+{0}, 1𝑙4𝑚, if 𝑏𝑖𝑘,𝑙4=1, then r𝑖𝑘,𝑙4 is a unit row vector.(viii)For all 𝑖𝒪𝑘, 𝑙𝑙1=1𝑢𝑖𝑘,𝑙1=0or1 and 𝑐𝑖𝑘,𝑙2=0, for all 𝑙2,1𝑙2𝑚.(ix)For all 𝑖𝒪𝑘, 𝑚𝑙2=1𝑐𝑖𝑘,𝑙2=0 or 1 and 𝑢𝑖𝑘,𝑙1=0, for all 𝑙1, 1𝑙1𝑙.(x)For all 𝑖𝒪𝑘, 𝑙𝑙3=1𝑣𝑖𝑘,𝑙3=0 or 1 and 𝑏𝑖𝑘,𝑙4=0, for all𝑙4, 1𝑙4𝑚.(xi)For all 𝑖𝒪𝑘, 𝑚𝑙4=1𝑏𝑖𝑘,𝑙4=0 or 1 and 𝑣𝑖𝑘,𝑙3=0, for all 𝑙3, 1𝑙3𝑙.(xii)For all 𝑖,𝑘, 1𝑖𝑛, 𝑘+{0}, 𝑙𝑙1=1𝑢𝑖𝑘,𝑙1+𝑚𝑙2=1𝑐𝑖𝑘,𝑙2=0 or 1 and 𝑙𝑙3=1𝑣𝑖𝑘,𝑙3+𝑚𝑙4=1𝑏𝑖𝑘,𝑙4=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 𝑥0,𝑥2𝑘1,𝑥2𝑘2+1, 𝑘1,𝑘2+{0}, at least one path from 𝑥0 to 𝑥2𝑘1 and at least one path from 𝑥0 to 𝑥2𝑘2+1 by applying proper mutagens in the correct order, with 𝑘1,𝑘2 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 Section 2.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 𝐽0𝑥0=min𝑢,𝑐,𝑣,𝑐𝔼,,𝑟,𝑟𝑁1𝑘=0𝑙𝑙1=1𝑛𝑖=1𝛼𝑙1𝑢𝑖𝑘,𝑙1+𝑁1𝑘=0𝑚𝑙2=1𝑛𝑖=1𝛽𝑙2𝑐𝑖𝑘,𝑙2costofapplyingchemicalmutagens+𝑁1𝑘=0𝑙𝑙3=1𝑛𝑖=1𝛼𝑙3𝑣𝑖𝑘,𝑙3+𝑁1𝑘=0𝑚𝑙4=1𝑛𝑖=1𝛽𝑙4𝑏𝑖𝑘,𝑙4costofapplyingradiativerays+𝑁𝑘=0𝑑𝑥𝑘,𝑥𝑑𝑘tracingcost,(17) with 𝑥0,𝑥𝑑𝑘𝑛, 1𝑘𝑁, 𝑛0 (mod3) given. The physical meaning of 𝑢𝑖𝑘,𝑙1, 𝑐𝑖𝑘,𝑙2, 𝑣𝑖𝑘,𝑙3, 𝑏𝑖𝑘,𝑙4, 𝑙1, 𝑙2, 𝑙3, 𝑙4 is the same as in Section 2.2.3. 𝛼𝑙1,𝛽𝑙2,𝛼𝑙3,𝛽𝑙4, for all𝑙1,𝑙2,𝑙3,𝑙4, 1𝑙1𝑙, 1𝑙2𝑚, 1𝑙3𝑙, 1𝑙4𝑚, are the corresponding cost of applying chemical mutagens and radiative rays indexed 𝑙1,𝑙2,𝑙3,𝑙4, respectively. {𝑥𝑑𝑘}𝑛×𝑛+{0} denotes the desired set at 𝑘th stage, generated by the DNA sequences representing the same amino acid sequence as𝑥𝑑𝑘, the desired state at 𝑘th step. And 𝑑(𝑥𝑘,{𝑥𝑑𝑘}) is the distance reference of the current state 𝑥𝑘 to the desired set {𝑥𝑑𝑘} at 𝑘th step. The final penalty, the distance reference from the final state to the desired set at 𝑘=𝑁, is included in the last term.

In general, 𝛽𝑙2,𝛽𝑙4𝛼𝑙1,𝛼𝑙3, for all 𝑙1,𝑙2,𝑙3,𝑙4, 1𝑙1𝑙, 1𝑙2𝑚, 1𝑙3𝑙, 1𝑙4𝑚, because physically 𝑂 is a set of nonsense bases and more details are necessary to convert an 𝑂 back to normal bases, for instance, the cost to identify the exact element in the set 𝑂. Our goal is to drive our system optimally from initial state 𝑥0 to the desired final set {𝑥𝑑𝑁} by applying a sequence of mutagens indexed with {𝑙1,𝑙2,𝑙3,𝑙4}, at problematic positions 𝑖, and in a correct order 𝑘.

In (17), the first four terms inside the expectation do not depend on random variables 𝑖k,𝑙1, 𝑟𝑖𝑘,𝑙2, 𝑖𝑘,𝑙3, and 𝑟𝑖𝑘,𝑙4, for all 𝑖,𝑘,𝑙1,𝑙2,𝑙3,𝑙4 as the treatment plan is computed based on the initial state 𝑥0. Given 𝑦𝑘, the updated treatment plan is computed accordingly but still not related to random variables. The last term inside expectation, 𝑁𝑘=0𝑑(𝑥𝑘,{𝑥𝑑𝑘}), 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𝐽0𝑥0=min{𝑢,𝑐,𝑣,𝑐}0,1,,𝑁1𝑁1𝑘=0𝑙𝑙1=1𝑛𝑖=1𝛼𝑙1𝑢𝑖𝑘,𝑙1+𝑁1𝑘=0𝑚𝑙2=1𝑛𝑖=1𝛽𝑙2𝑐𝑖𝑘,𝑙2+𝑁1𝑘=0𝑙𝑙3=1𝑛𝑖=1𝛼𝑙3𝑣𝑖𝑘,𝑙3+𝑁1𝑘=0𝑚𝑙4=1𝑛𝑖=1𝛽𝑙4𝑏𝑖𝑘,𝑙4+𝑁𝑘=0𝔼{,𝑟,,𝑟}0,1,,𝑁1𝑑𝑥𝑘,𝑥𝑑𝑘,(18) subject to𝑥𝑘+1=𝐼+𝑙𝑙3=1𝑛𝑖=1𝑣𝑖𝑘,𝑙3𝑖𝑘,𝑙3𝑠𝑒𝑖𝑒𝑇𝑖×𝐼+𝑙𝑙1=1𝑛𝑖=1𝑢𝑖𝑘,𝑙1𝑖𝑘,𝑙1𝑠𝑒𝑖𝑒𝑇𝑖𝑥𝑘+𝐼+𝑙𝑙3=1𝑛𝑖=1𝑣𝑖𝑘,𝑙3𝑖𝑘,𝑙3𝑠𝑒𝑖𝑒𝑇𝑖𝑚𝑙2=1𝑖𝒪𝑘𝑐𝑖𝑘,𝑙2𝑟𝑖𝑘,𝑙2𝑤𝑒𝑖+𝑚𝑙4=1𝑖𝒪𝑘𝑏𝑖𝑘,𝑙4𝑟𝑖𝑘,𝑙4𝑤𝑒𝑖,𝑦𝑘=𝑥𝑘.(19) 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, 𝑑(𝜑1,𝜑2),𝜑1,𝜑23, needs to fulfill the biological requirements as below.(1)Nonnegativity: the distance reference between any two codons is either positive or zero. Mathematically, 𝑑3×3+{0}, 𝑑(𝜑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, 𝑑(𝜑1,𝜑2)=𝑑(𝜑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, 𝑑(𝜑1,𝜑2)=0 if and only if 𝜑1=𝜑2. The redundancy in genetic codes implies 𝑑(𝜑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 Table 5. We ignore codons containing 𝑂 since their chemical and physical properties cannot be found in literature.

From Table 5, 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𝑑𝜉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𝜉1PHvalueof𝜉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,(20) where 𝜉1,𝜉2 are two amino acids.

𝑑(𝜉1,𝜉2) is then assigned to 𝑑(𝜑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 for 𝑥0 to {𝑥𝑑𝑁} is also the optimal policy for the tail problem, from 𝑥𝑞 to {𝑥𝑑𝑁}.

The tail problem is defined as𝐽𝑞𝑥𝑞=min{𝑢,𝑐,𝑣,𝑐}𝑞,𝑞+1,,𝑁1𝑁1𝑘=𝑞𝑙𝑙1=1𝑛𝑖=1𝛼𝑙1𝑢𝑖𝑘,𝑙1+𝑁1𝑘=𝑞𝑚𝑙2=1𝑛𝑖=1𝛽𝑙2𝑐𝑖𝑘,𝑙2+𝑁1𝑘=𝑞𝑙𝑙3=1𝑛𝑖=1𝛼𝑙3𝑣𝑖𝑘,𝑙3+𝑁1𝑘=𝑞𝑚𝑙4=1𝑛𝑖=1𝛽𝑙4𝑏𝑖𝑘,𝑙4+𝑁𝑘=𝑞𝔼{,𝑟,,𝑟}𝑞,𝑞+1,,𝑁1𝑑𝑥𝑘,𝑥𝑑𝑘.(21) The iterative update equation to find optimal policy can be expressed by (22), according to the dynamic programming algorithm in [31]. 𝐽𝑁𝑥𝑁=𝑑𝑥𝑁,𝑥𝑑𝑁,𝐽𝑞𝑥𝑞=min𝑢𝑞,𝑐𝑞,𝑣𝑞,𝑐𝑞𝔼𝑞,𝑟𝑞,𝑞,𝑟𝑞𝑙𝑙1=1𝑛𝑖=1𝛼𝑙1𝑢𝑖𝑞,𝑙1+𝑚𝑙2=1𝑛𝑖=1𝛽𝑙2𝑐𝑖𝑞,𝑙2+𝑙𝑙3=1𝑛𝑖=1𝛼𝑙3𝑣𝑖𝑞,𝑙3+𝑚𝑙4=1𝑛𝑖=1𝛽𝑙4𝑐𝑖𝑞,𝑙4+𝑑𝑥𝑞,𝑥𝑑𝑞+𝐽𝑞+1𝑥𝑞+1=min𝑢𝑞,𝑐𝑞,𝑣𝑞,𝑐𝑞𝑙𝑙1=1𝑛𝑖=1𝛼𝑙1𝑢𝑖𝑞,𝑙1+𝑚𝑙2=1𝑛𝑖=1𝛽𝑙2𝑐𝑖𝑞,𝑙2+𝑙𝑙3=1𝑛𝑖=1𝛼𝑙3𝑣𝑖𝑞,𝑙3+𝑚𝑙4=1𝑛𝑖=1𝛽𝑙4𝑐𝑖𝑞,𝑙4+𝔼𝑞,𝑟𝑞,𝑞,𝑟𝑞𝑑𝑥𝑞,𝑥𝑑𝑞+𝐽𝑞+1𝑥𝑞+1,𝑞=0,1,,𝑁1.(22)

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 and 𝑂 because of high-cost 𝛽𝑙2 and the unavailable chemical and physical properties for codons containing 𝑂.

The distance reference between codons used in Sections 3.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𝑑𝜓1,𝜓2=0,if𝑥𝑁=𝑥𝑑𝑁,,if𝑥𝑁𝑥𝑑𝑁,(23) with 𝜓1,𝜓2{0}, where {0} denotes the set excluding the element 0.

Our optimal control problem for point mutations is𝐽0𝑥0=min𝑢𝑘,𝑙1,0𝑘𝑁,1𝑙1𝑙𝑁1𝑘=0𝑙𝑙1=1𝛼𝑙1𝑢𝑘,𝑙1,(24) subject to𝑥𝑘+1=𝐼+𝑙𝑙1=1𝑢𝑘,𝑙1𝑠𝑥𝑘,𝑥𝑁=𝑥𝑑𝑁,(25) with 𝑥0 given, 𝑥𝑘{0}.

Suppose that there are 12 kinds of mutagens (𝑙1=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.

The elements along the antidiagonal of Table 7, 𝑢𝐴𝑇,𝑢𝐺𝐶,𝑢𝐶𝐺, and 𝑢𝑇𝐴, 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, 𝛼𝐴𝑇=𝛼𝐺𝐶=𝛼𝐶𝐺=𝛼𝑇𝐶=0. The equivalence relationship between subscription in two nucleotide bases and subscription in integer 𝑙1(𝜓1𝜓2){𝐴,𝑇,𝐺,𝐶}×{𝐴,𝑇,𝐺,𝐶}{integersfrom1to12} is defined by Table 6.

Under the above assumptions, we can write update equation for optimal policy explicitly as𝐽𝑞𝑥𝑞=min𝑢𝑞,𝑙1𝛼𝑥𝑞𝜓+𝐽𝑞+1(𝜓),𝜓{𝐴,𝑇,𝐺,𝐶}{0}.(26)

Proposition 7. For the same 𝑥𝑑𝑁, 𝐽𝑞(𝜓)𝐽𝑞+1(𝜓), for all 𝑞, 0𝑞𝑁1, for all 𝜓{𝐴,𝑇,𝐺,𝐶}, where 𝜓 denotes the complementary base of 𝜓. If, in addition, the system is completely controllable, 𝑀, s.t. 𝐽𝑀(𝜓) is the global minimum and for all 𝑞𝑀, 𝐽𝑞(𝜓)=𝐽𝑀(𝜓) if 𝑀𝑞 1 (mod 2), and 𝐽𝑞(𝜓)=𝐽𝑀(𝜓) if 𝑀𝑞0 (mod 2). In our example, 𝑀𝑁6.

Proof. This first part is due to the zero cost for the transfers between complementary bases in the consecutive steps.
For any 0𝑞𝑁1, the relationship between minimal costs in consecutive steps is shown in (26). Since 𝜓{𝐴,𝑇,𝐺,𝐶}, 𝛼𝜓𝜓+𝐽𝑞+1(𝜓) is one of the four elements in the set from which the 𝐽𝑞(𝜓) is picked. Moreover, 𝛼𝜓𝜓=0. Therefore, 𝐽𝑞+1(𝜓) is one of the four elements in the set. Since 𝐽𝑞(𝜓) is the minimum picking for a set containing 𝐽𝑞+1(𝜓), we conclude that 𝐽𝑞(𝜓)𝐽𝑞+1(𝜓).
The 𝑀 value in our example is proved by brute force method, that is, 𝐽𝑁6(𝜓) is a guaranteed global minimum. For completely controllable systems, this 𝑀 always exists. The existence of 𝑀 implies that for without limitation in the number of steps, we can reach the global optimal in 𝑁𝑀 steps, 6 steps in our example.
Suppose that 𝑞=𝑀, 𝐽𝑀(𝜓) is the global minimum, thus 𝐽𝑀1(𝜓)𝐽𝑀(𝜓). However,𝐽𝑀1(𝜓)𝐽𝑀(𝜓) according to the first part of the proposition. Therefore, 𝐽𝑀1(𝜓)=𝐽𝑀(𝜓) for the same 𝑥𝑑𝑁. Therefore, 𝐽𝑀1(𝜓) is also a global minimum.
By backward induction, suppose for 𝑞=𝑞1, the statement is true, that is, 𝐽𝑞11(𝜓)=𝐽𝑞1(𝜓) is the global optimal either from 𝑥𝑞11=𝜓 or 𝑥𝑞1=𝜓 to 𝑥𝑑𝑁. Obviously, for 𝑞=𝑞11, the statement is still true. Therefore, 𝐽𝑞(𝜓)=𝐽𝑞2(𝜓)=𝐽𝑞1(𝜓),𝜓{𝐴,𝑇,𝐺,𝐶}, for all 𝑞, 2𝑞𝑀.

In the proof of global minimum that can be reached in the finite step in Proposition 7, we also discover Proposition 8. Here, 𝐽𝑞(𝑥𝑞,𝑥𝑑𝑁) denotes the optimal cost from 𝑥𝑞 to 𝑥𝑑𝑁.

Proposition 8. Given two single base mutation optimal control problems, with the same fixed 𝑁, with and desired final states complementary to each other. If 𝐽𝑀(𝜓,𝑥𝑑𝑁) is the global minimum, then 𝐽𝑀(𝜓,𝑥𝑑𝑁) is also the global minimum, that is, the global minimum of both systems is reach at the same stage 𝑀. Moreover, for all 𝑞, 0𝑞𝑀, 𝐽𝑞𝜓,𝑥𝑑𝑁=𝐽𝑞𝜓,𝑥𝑑𝑁,𝜓,𝑥𝑑𝑁{𝐴,𝑇,𝐺,𝐶}.(27)

Physically, Proposition 8 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 Table 9. It is a numerical assignment to Table 8. Since we apply mutagens before the replication starts, 𝑢𝐴𝐴 actually transfer 𝐴 to 𝑇 and then to 𝐴 by replication. For simplicity, we just use the 𝑘th and (𝑘+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, 𝛼𝐴𝐶,𝛼𝐶𝐴,𝛼𝐺𝑇,𝛼𝑇𝐺 is smaller than other mutagens, except artificial ones.

If we use 𝜒 to denote the costs of mutagens as listed in Table 9, then 𝜒=5.216.602.3306.158.9503.824.6109.177.2400.645.0910.28=𝛼𝐴𝐴𝛼𝐴𝐺𝛼𝐴𝐶𝛼𝐴𝑇𝛼𝐺𝐴𝛼𝐺𝐺𝛼𝐺𝐶𝛼𝐺𝑇𝛼𝐶𝐴𝛼𝐶𝐺𝛼𝐶𝐶𝛼𝐶𝑇𝛼𝑇𝐴𝛼𝑇𝐺𝛼𝑇𝐶𝛼𝑇𝑇𝛼1𝛼2𝛼30𝛼4𝛼50𝛼6𝛼70𝛼8𝛼90𝛼10𝛼11𝛼12.(28)

Running the dynamic programming for every pair of (𝑥𝑞,𝑥𝑑𝑁){𝐴,𝑇,𝐺,𝐶}×{𝐴,𝑇,𝐺,𝐶}, 𝑁=9. Here, we sightly modify our notation. We use 𝐽𝑞(𝑥𝑞,𝑥𝑑𝑁) to denote the optimal cost from 𝑥𝑞 to 𝑥𝑑𝑁. Then, 𝐽𝑞=𝐽𝑞(𝐴,𝐴)𝐽𝑞(𝐴,𝐺)𝐽𝑞(𝐴,𝐶)𝐽𝑞(𝐴,𝑇)𝐽𝑞(𝐺,𝐴)𝐽𝑞(𝐺,𝐺)𝐽𝑞(𝐺,𝐶)𝐽𝑞(𝐺,𝑇)𝐽𝑞(𝐶,𝐴)𝐽𝑞(𝐶,𝐺)𝐽𝑞(𝐶,𝐶)𝐽𝑞(𝐶,𝑇)𝐽𝑞(𝑇,𝐴)𝐽𝑞(𝑇,𝐺)𝐽𝑞(𝐶,𝐶)𝐽𝑞(𝑇,𝑇).(29)

The path to reach the optimal cost is denoted by 𝑃𝑞=𝑃𝑞(𝐴,𝐴)𝑃𝑞(𝐴,𝐺)P𝑞(𝐴,𝐶)𝑃𝑞(𝐴,𝑇)𝑃𝑞(𝐺,𝐴)𝑃𝑞(𝐺,𝐺)𝑃𝑞(𝐺,𝐶)𝑃𝑞(𝐺,𝑇)𝑃𝑞(𝐶,𝐴)𝑃𝑞(𝐶,𝐺)𝑃𝑞(𝐶,𝐶)𝑃𝑞(𝐶,𝑇)𝑃𝑞(𝑇,𝐴)𝑃𝑞(𝑇,𝐺)𝑃𝑞(𝑇,𝐶)𝑃𝑞(𝑇,𝑇),(30)

where 𝑃𝑞(𝑥𝑞,𝑥𝑑𝑁) is the (𝑞+1)th state from 𝑥𝑞 to 𝑥𝑑𝑁, that is, 𝑥𝑞+1=𝑃𝑞(𝑥𝑞,𝑥𝑑𝑁).

The simulation results are shown as below, including optimal costs for all possible transfer pairs (𝑥𝑞,𝑥𝑑𝑁){𝐴,𝑇,𝐺,𝐶}×{𝐴,𝑇,𝐺,𝐶}, 𝐽𝑞, 0𝑞8, graphical representation in Figure 6, and optimal path 𝑃𝑞, 0𝑞7. 𝐽0=,𝐽1=00.645.095.213.8206.796.156.156.7903.825.215.090.640,𝐽2=,𝐽3=00.645.095.213.8206.796.156.156.7903.825.215.090.640,𝐽4=,𝐽5=00.645.095.213.8206.796.156.156.7903.825.215.090.640,𝐽6=,𝐽7=00.645.095.213.8208.486.156.157.8803.825.215.090.640,𝐽8=5.216.602.3306.158.9503.824.6109.177.2400.645.0910.28,(31)𝑃0=𝐴,𝑇𝑇𝑇𝑇𝐴,𝐶𝐴,𝐶𝐶𝐶,𝑇𝐺𝐺𝐺𝐺𝐴𝐴,𝐺𝐴,𝐶𝐴,𝑃1=𝑇𝑇𝑇𝐴,𝑇𝐶,𝑇𝐶𝐴,𝐶𝐴,𝐶𝐺𝐺𝐺𝐺𝐴𝐴,𝐶𝐴,𝐺𝐴,𝑃2=𝐴,𝑇𝑇𝑇𝑇𝐴,𝐶𝐴,𝐶𝐶𝐶,𝑇𝐺𝐺𝐺𝐺𝐴𝐴,𝐺𝐴,𝐶𝐴,𝑃3=𝑇𝑇𝑇𝐴,𝑇𝐶,𝑇𝐶𝐴,𝐶𝐴,𝐶𝐺𝐺𝐺𝐺𝐴𝐴,𝐶𝐴,𝐺𝐴,𝑃4=𝐴,𝑇𝑇𝑇𝑇𝐴,𝐶𝐴,𝐶𝐶𝐶,𝑇𝐺𝐺𝐺𝐺𝐴𝐴,𝐺𝐴,𝐶𝐴,𝑃5=𝑇𝑇𝑇𝐴,𝑇𝐶,𝑇𝐶𝐴𝐴,𝐶𝐺𝐺𝐺𝐺𝐴𝐴,𝐶𝐴,𝐺𝐴,𝑃6=