Time-Adaptive Determination of Drug Efficacy in Mathematical Model of HIV Infection

The paper considers a time-adaptive finite element method for determination of drug efficacy in a parameter identification problem (PIP) for a system of ordinary differential equations (ODE) that describes dynamics of the primary human immunodeficiency virus (HIV) infection with drug therapy. Tikhonov’s regularization method, optimization approach and finite element method to solve this problem are presented. A posteriori error estimates in the Tikhonov’s functional and reconstructed parameter are derived. Based on these estimates a time adaptive algorithm is formulated and numerically tested for different scenarios of noisy observations of virus population function. Numerical results show a significant improvement of reconstruction of drug efficacy parameter using the local time-adaptive mesh refinement method compared to the gradient method applied on a uniform time mesh. Introduction Mathematical modeling of the immune processes is an essential part of the research in immunology [15, 16, 22]. Despite the emergence of a great amount of new high-performance methods for experimental analysis of the immunity, the results of mathematical modeling are relevant, in particular, in clinical practice in order to work out optimal, individually customized strategies for treatment of the pathological process (bacterial/viral infections or tumor growth). It is well known that physiological parameters vary between individual patients and thus, personalized treatment approaches require the development of robust and efficient * L. Beilina larisa@chalmers.se M. Eriksson martin.eriksson.2@gu.se I. Gainova gajnova@math.nsc.ru 1 Department of Mathematical Sciences, Chalmers University of Technology and University of Gothenburg, 42196 Göteborg, Sweden 2 Department of Marine Sciences, University of Gothenburg, 40530 Göteborg, Sweden 3 Sobolev Institute of Mathematics, Novosibirsk 630090, Russia Differential Equations and Dynamical Systems 1 3 parameter estimation methods to assimilate individual data with mathematical models [6]. The problems of parameter identification in mathematical models are often non-linear and ill-posed, and thus challenging to solve numerically [45, 46]. The computational algorithms for parameter identification that we work on, are based on an adaptive time-mesh refinement [13] for coefficient inverse problems (CIPs). The main idea of our work is adoptation of results for space mesh refinement developed in [13] for solution of hyperbolic CIPs, to the parameter identification problem (PIP) for reconstruction of parameters on the time mesh. More precisely, first, we determine candidate parameter at known initial (coarse) time partition. Then we refine time-mesh locally only at a such time intervals where a posteriori error indicator is large and compute new time-dependent control functions on a new time mesh until the error in the computed residual is reduced to the desired accuracy. The adaptive finite element method has shown that it significantly improves reconstruction of parameters when solving the coefficient inverse problems for hyperbolic PDE [5, 7, 8, 12, 13]. We note that the main goal of our work is to present mathematical framework of a posteriori error estimation for solution of PIP’s and to show usefulness of the time-adaptive error control for determination of parameters in PIP which we demonstrate on the example of the model of HIV infection. More than 35 years have passed since the discovery of the etiological agent of AIDS—human immunodeficiency virus (HIV), nevertheless, the problem of the spread of HIV infection, treatment and quality of life for people living with HIV still remains relevant: the number of newly infected does not decrease. Appearing of highly active antiretroviral therapy (HAART) in 1996-1997 has led to a significant improvement in the quality of life of patients, has caused a clear decrease in AIDS-related diseases and mortality. HAART provides treatment protocols with using combinations (two or more) of antiviral drugs which affect both the different stages of viral replication and prevent HIV from entering the host cell. We took the simplified model of HIV infection proposed in [42] as an illustrative example to show effectiveness and robustness of identification time-dependent parameter drug efficacy using a local time-mesh refinement algorithm. This work is a continuation of the works by authors [9, 10] where was introduced the time-adaptive finite element method for parameter identification problems. Compared to other optimal control algorithms for solution of PIP, see, for example, works [1, 2, 26] and references therein, our time-adaptive algorithms are based on rigorous finite element a posteriori error analysis for the error in the Tikhonov’s functional or for the error in the reconstructed parameter. The same approach can be applied to the solution of any other PIP and particularly, for more complicated models of HIV infection which involve more unknown functions and parameters [3, 27, 34–36, 40, 41, 48]. However, these models are much complicated compared with studied here model of [42], and can be considered as topic for a future research. As was mentioned above, in [10] was studied the optimal control problem of reconstruction of drug efficacy in the model of HIV infection when measurements of all functions in time of the model ODE system were used what, actually, is not realistic problem. Moreover, numerical simulations were not presented in [10]. In the current work we fill this gap and study a more realistic case when instead of measuring all four functions in ODE system, we take measurements only of the virus population function. New a posteriori error estimate between regularized and computed drug efficacy is derived. Based on this estimate, a time adaptive algorithm is formulated and numerically tested on the optimal determination of drug efficacy in time domain from noisy measurements of virus population function. Extended numerical studies for different noise levels in the virus population function and for different values of time for initial observation of this function are presented in Differential Equations and Dynamical Systems 1 3 the preprint version of this work [11]. We note that compared to the current work proofs of Theorems 1, 2, and 3 are not presented in [11], as well as study on stability of solution of the forward problem and study of ill-posedness of the PIP problem. The time-adaptive method proposed in this paper can eventually be used by clinicians to determine the drug-response for each treated individual. The exact knowledge of the personal drug efficacy can aid in the determination of the most suitable drug as well as the most optimal dose to an individual, in the long run resulting in a personalized treatment with maximum efficacy and minimum adverse drug reactions. The outline of the paper is as follows. The biological description of the mathematical model is given in “The Mathematical Model and Its Biological Description”. “Inverse Problems and Ill-Posedness” is based on material of the Master’s thesis [25] and discusses ill-posedness of parameter identification problems and Tikhonov’s regularization method for their solutions. In “The Parameter Identification Problem” the parameter identification problem is formulated. The optimization method to solve the parameter identification problem is presented in “Optimization Method”. The finite element method is formulated in “Finite Element Discretization” and a posteriori error estimates are derived in “A Posteriori Error Estimates”. An adaptive algorithm for solution of PIP is presented in “Algorithms for Solution of PIP”. Finally, in “Numerical Results” numerical examples illustrate effectiveness of the proposed time-adaptive algorithm. The Mathematical Model and Its Biological Description The main cellular targets of HIV are the immune system cells that have CD4 receptors at their surface, called CD4+ T-cells. HIV differs from other viruses by a high mutation level, the mutation rate is 10–10 per nucleotide during one replication cycle [30]. High genetic variability allows HIV to skillfully avoid humoral and cellular defense factors and the effects of drugs. The constant presence of large reservoirs of latently infected cells is one of the important feature of the pathogenesis of HIV infection. Another paradoxical feature of HIV is that activation of the immune system does not lead to a suppression of virus multiplication, but rather to activation of latently infected cells, which start to produce new viral particles. These factors are the major obstacles for antiviral therapy and the development of efficient vaccines [18, 19]. According to the latest data on HIV (UNAIDS, 2020) [47], there are currently 38 million people living with HIV globally and 1.7 million people became newly infected with HIV in 2019. The HIV life cycle starts with attachment of the viral envelope protein gp120 to the cell surface via interaction with CD4 receptor. In general case for living organisms the genetic information goes from the storage in DNA through messenger RNA (mRNA) to protein synthesis in the ribosomes. The process of converting the genetic information from DNA to mRNA is called transcription [25]. In the case of retroviruses, such as HIV, HIV’s genetic information is encoded in form of RNA. Having fused with the cell membrane, HIV releases its genetic material (viral RNA) and enzymes into the CD4+ T-cell. Here viral RNA is reversely transcribed into HIV DNA, which is compatible with genetic material of the host cell [reverse transcription (RT)]. To perform the reverse transcription of RNA into DNA, HIV carries its own enzyme called reverse transcriptase. Viral DNA is transported to the cell’s nucleus and incorporated into the DNA of the host cell (integration). This process is made possible by the enzyme integrase. The individual components of HIV are then produced within the CD4+ T-cell. Differential Equations and Dynamical Systems 1 3 The individual components of HIV are then assembled together to make new HIV viruses. This process depends on the enzyme protease. The newly matured HIV particles are released from the CD4+ T-cell. These are ready to infect other CD4+ T-cells and begin the replication process all over again. The process of HIV life cycle described above is illustrated in Fig.  1. Antiretroviral drugs blocking the enzyme reverse transcriptase (called Reverse Transcriptase Inhibitors) will be able to prevent the production of new viruses [18, 37]. Our basic mathematical model in this work is the model proposed in [42] which describes the effect of Reverse Transcriptase Inhibitor (RTI) on the dynamics of HIV infection. In this model the infected class of CD4+ T-cells is subdivided into two subclasses: pre-RT class and post-RT class. Pre-RT class consists of the infected CD4+ T-cells in which reverse transcription is not completed, and post-RT class consists of those infected CD4+ T-cells where the reverse transcription is completed such that they are capable to produce virus. Throughout the paper we denote by ΩT = [0,T] the time domain for T > 0 , where T is the final observation time. The mathematical model which we use in this note, is: with initial conditions (1) ⎧⎪⎨⎪⎪⎪⎩ du1 dt = f1(u(t), (t)) = s − ku1(t)u4(t) − u1(t) + ( (t) + b)u2(t), du2 dt = f2(u(t), (t)) = ku1(t)u4(t) − ( 1 + + b)u2(t), du3 dt = f3(u(t), (t)) = (1 − (t)) u2(t) − u3(t), du4 dt = f4(u(t), (t)) = N u3(t) − cu4(t), Fig. 1 The HIV life cycle Differential Equations and Dynamical Systems 1 3 In system (1) the functions are defined as follows: • u1(t) – uninfected target cells population • u2(t) – infected target cells before Reverse Transcription (pre-RT class) • u3(t) – infected target cells in which Reverse Transcription is completed and they are capable of producing virus (post-RT class) • u4(t) is the virus population function. The parameters which are used in system (1) are described in Table  1 and are taken from the specialized literature [32, 33]. The initial data (2) is chosen such that they satisfy two steady states which we discuss in “Existence and Lyapunov Stability of Steady State”. Figure  2 illustrates the effect of Reverse Transcriptase Inhibitor (RTI) on the dynamics of HIV infection for the mathematical model (1). The system (1) can be presented in the following compact form: where we have denoted all involved functions as (2) u1(0) = u 0 1 = 300 mm, u2(0) = u 0 2 = 10 mm, u3(0) = u 0 3 = 10 mm, u4(0) = u 0 4 = 10 mm. (3) { du dt = f (u(t), (t)) t ∈ [0, T], u(0) = u, Fig. 2 The effect of reverse transcriptase inhibitor (RTI) on the dynamics of HIV infection described by the model problem (1) Differential Equations and Dynamical Systems 1 3 Here, (⋅) denotes transposition operator. In the model (1) we assume that f ∈ C1(ΩT ) is Lipschitz continuous and the function (t) ∈ C(ΩT ) represents the unknown drug efficacy which belongs to the set of admissible functions M : The control parameter is dosage of the reverse transcriptase inhibitor. This parameter protects the cells and prevents infection. In this work we assume that all parameters in system (1) are constants except the control parameter which depends on time, e.,e., = (t), t ∈ [0, T] . This means that the control parameter = (t) tells us which dosage of the reverse transcriptase inhibitor should be given to the concrete patient at any time moment t for t ∈ [0, T] . Thus, personalized determination of this parameter for every individual is of vital importance for treatment of HIV. Existence and Lyapunov Stability of Steady State Let us now assume that the parameter in system (1) is constant. That is, (t) ≡ c ∈ (0, 1) . Setting du dt in (1) to zero and solving for u1 , u2 , u3 and u4 we can see that there are two possible steady states: an infected and an uninfected one [42]. The uninfected steady state is given by (4) u = u(t) = (u1(t), u2(t), u3(t), u4(t)) T , u = (u1(0), u2(0), u3(0), u4(0)) T , du dt = ( u1 t , u2 t , u3 t , u4 t )T , f (u(t), (t)) = (f1, f2, f3, f4) T (u(t), (t)) = (f1(u1,... , u4, (t)),... , f4(u1,... , u4, (t))) T . (5) M = { (t) ∶ (t) ∈ [0, 1] inΩT , (t) = 0 outside of ΩT}. Table 1 Parameters dataset Parameter Value Units Description s 10 mm day Inflow rate of T cells 0.01 day Natural death rate of T cells k 2.4E−5 mmday Interaction-infection rate of T cells 1 0.015 day Death rate of infected cells 0.4 day Transition rate from pre-RT infected T cells class to post-RT class b 0.05 day Reverting rate of infected cells return to uninfected class 0.26 day Death rate of actively infected cells c 2.4 day Clearance rate of virus N 1000 vir/cell Total number of viral particles produced by an infected cell [0, 1] NA Control, dosage of the reverse transcriptase inhibitor (drug efficacy) Differential Equations and Dynamical Systems 1 3 and the infected steady state is achieved when In [42] was shown that the infected steady state can exist only when is less than the following critical value For our system of parameters, presented in Table 1, this critical value is crit ≈ 0.88375 . Whenever ≥ crit only the uninfected steady state can exist. Plugging in the values of Table 1 into (7) for η < ηcrit , or (6) if ≥ crit , we obtain the numerical values for solutions (u1, u2, u3, u4) of (1) presented in the Table 2. Stability of Solutions Let us define (6) ⎧⎨⎪⎩ u1 = s , u2 = 0, u3 = 0, u4 = 0, (7) ⎧⎪⎨⎪⎪⎩ u1 = ( 1+ +b)c N k(1− ) , u2 = s− u1 1+ (1− ) , u3 = (1− )u2 , u4 = N u3 c . (8) crit = 1 − c( 1 + + b) N ks . (9) ⎧⎪⎨⎪⎪⎩ m = min{ , 1}, Ξ = s m , Φ ∶= Φ( ) = s(1− ) m , Ψ ∶= Ψ( ) = N s(1− ) mc , Table 2 Stable steady states for different values of , while keeping the other parameters fixed u 1 u 2 u 3 u 4 0.0 116 21 33 3549 0.1 129 23 32 3483 0.2 145 26 31 3402 0.3 166 28 30 3298 0.4 194 32 29 3162 0.5 233 36 27 2975 0.6 291 41 25 2702 0.7 388 45 21 2269 0.8 581 44 14 1469 0.9 1000 0 0 0 1.


Introduction
Mathematical modeling of the immune processes is an essential part of the research in immunology [15,16,22]. Despite the emergence of a great amount of new high-performance methods for experimental analysis of the immunity, the results of mathematical modeling are relevant, in particular, in clinical practice in order to work out optimal, individually customized strategies for treatment of the pathological process (bacterial/viral infections or tumor growth).
It is well known that physiological parameters vary between individual patients and thus, personalized treatment approaches require the development of robust and efficient 1 3 parameter estimation methods to assimilate individual data with mathematical models [6]. The problems of parameter identification in mathematical models are often non-linear and ill-posed, and thus challenging to solve numerically [45,46]. The computational algorithms for parameter identification that we work on, are based on an adaptive time-mesh refinement [13] for coefficient inverse problems (CIPs). The main idea of our work is adoptation of results for space mesh refinement developed in [13] for solution of hyperbolic CIPs, to the parameter identification problem (PIP) for reconstruction of parameters on the time mesh. More precisely, first, we determine candidate parameter at known initial (coarse) time partition. Then we refine time-mesh locally only at a such time intervals where a posteriori error indicator is large and compute new time-dependent control functions on a new time mesh until the error in the computed residual is reduced to the desired accuracy. The adaptive finite element method has shown that it significantly improves reconstruction of parameters when solving the coefficient inverse problems for hyperbolic PDE [5,7,8,12,13].
We note that the main goal of our work is to present mathematical framework of a posteriori error estimation for solution of PIP's and to show usefulness of the time-adaptive error control for determination of parameters in PIP which we demonstrate on the example of the model of HIV infection. More than 35 years have passed since the discovery of the etiological agent of AIDS-human immunodeficiency virus (HIV), nevertheless, the problem of the spread of HIV infection, treatment and quality of life for people living with HIV still remains relevant: the number of newly infected does not decrease. Appearing of highly active antiretroviral therapy (HAART) in 1996-1997 has led to a significant improvement in the quality of life of patients, has caused a clear decrease in AIDS-related diseases and mortality. HAART provides treatment protocols with using combinations (two or more) of antiviral drugs which affect both the different stages of viral replication and prevent HIV from entering the host cell.
We took the simplified model of HIV infection proposed in [42] as an illustrative example to show effectiveness and robustness of identification time-dependent parameter drug efficacy using a local time-mesh refinement algorithm. This work is a continuation of the works by authors [9,10] where was introduced the time-adaptive finite element method for parameter identification problems. Compared to other optimal control algorithms for solution of PIP, see, for example, works [1,2,26] and references therein, our time-adaptive algorithms are based on rigorous finite element a posteriori error analysis for the error in the Tikhonov's functional or for the error in the reconstructed parameter. The same approach can be applied to the solution of any other PIP and particularly, for more complicated models of HIV infection which involve more unknown functions and parameters [3, 27, 34-36, 40, 41, 48]. However, these models are much complicated compared with studied here model of [42], and can be considered as topic for a future research.
As was mentioned above, in [10] was studied the optimal control problem of reconstruction of drug efficacy in the model of HIV infection when measurements of all functions in time of the model ODE system were used what, actually, is not realistic problem. Moreover, numerical simulations were not presented in [10]. In the current work we fill this gap and study a more realistic case when instead of measuring all four functions in ODE system, we take measurements only of the virus population function. New a posteriori error estimate between regularized and computed drug efficacy is derived. Based on this estimate, a time adaptive algorithm is formulated and numerically tested on the optimal determination of drug efficacy in time domain from noisy measurements of virus population function. Extended numerical studies for different noise levels in the virus population function and for different values of time for initial observation of this function are presented in

The Mathematical Model and Its Biological Description
The main cellular targets of HIV are the immune system cells that have CD4 receptors at their surface, called CD4+ T-cells. HIV differs from other viruses by a high mutation level, the mutation rate is 10 −5 -10 −4 per nucleotide during one replication cycle [30]. High genetic variability allows HIV to skillfully avoid humoral and cellular defense factors and the effects of drugs. The constant presence of large reservoirs of latently infected cells is one of the important feature of the pathogenesis of HIV infection. Another paradoxical feature of HIV is that activation of the immune system does not lead to a suppression of virus multiplication, but rather to activation of latently infected cells, which start to produce new viral particles. These factors are the major obstacles for antiviral therapy and the development of efficient vaccines [18,19]. According to the latest data on HIV (UNAIDS, 2020) [47], there are currently 38 million people living with HIV globally and 1.7 million people became newly infected with HIV in 2019.
The HIV life cycle starts with attachment of the viral envelope protein gp120 to the cell surface via interaction with CD4 receptor. In general case for living organisms the genetic information goes from the storage in DNA through messenger RNA (mRNA) to protein synthesis in the ribosomes. The process of converting the genetic information from DNA to mRNA is called transcription [25]. In the case of retroviruses, such as HIV, HIV's genetic information is encoded in form of RNA. Having fused with the cell membrane, HIV releases its genetic material (viral RNA) and enzymes into the CD4+ T-cell. Here viral RNA is reversely transcribed into HIV DNA, which is compatible with genetic material of the host cell [reverse transcription (RT)]. To perform the reverse transcription of RNA into DNA, HIV carries its own enzyme called reverse transcriptase. Viral DNA is transported to the cell's nucleus and incorporated into the DNA of the host cell (integration). This process is made possible by the enzyme integrase. The individual components of HIV are then produced within the CD4+ T-cell.
The individual components of HIV are then assembled together to make new HIV viruses. This process depends on the enzyme protease. The newly matured HIV particles are released from the CD4+ T-cell. These are ready to infect other CD4+ T-cells and begin the replication process all over again. The process of HIV life cycle described above is illustrated in Fig. 1. Antiretroviral drugs blocking the enzyme reverse transcriptase (called Reverse Transcriptase Inhibitors) will be able to prevent the production of new viruses [18,37]. Our basic mathematical model in this work is the model proposed in [42] which describes the effect of Reverse Transcriptase Inhibitor (RTI) on the dynamics of HIV infection. In this model the infected class of CD4+ T-cells is subdivided into two subclasses: pre-RT class and post-RT class. Pre-RT class consists of the infected CD4+ T-cells in which reverse transcription is not completed, and post-RT class consists of those infected CD4+ T-cells where the reverse transcription is completed such that they are capable to produce virus.
Throughout the paper we denote by Ω T = [0, T] the time domain for T > 0 , where T is the final observation time. The mathematical model which we use in this note, is: In system (1) the functions are defined as follows: • u 1 (t) -uninfected target cells population • u 2 (t) -infected target cells before Reverse Transcription (pre-RT class) • u 3 (t) -infected target cells in which Reverse Transcription is completed and they are capable of producing virus (post-RT class) • u 4 (t) is the virus population function.
The parameters which are used in system (1) are described in Table 1 and are taken from the specialized literature [32,33]. The initial data (2) is chosen such that they satisfy two steady states which we discuss in "Existence and Lyapunov Stability of Steady State". Figure 2 illustrates the effect of Reverse Transcriptase Inhibitor (RTI) on the dynamics of HIV infection for the mathematical model (1).
The system (1) can be presented in the following compact form: where we have denoted all involved functions as du dt
In the model (1) we assume that f ∈ C 1 (Ω T ) is Lipschitz continuous and the function (t) ∈ C(Ω T ) represents the unknown drug efficacy which belongs to the set of admissible functions M : The control parameter is dosage of the reverse transcriptase inhibitor. This parameter protects the cells and prevents infection. In this work we assume that all parameters in system (1) are constants except the control parameter which depends on time, e.,e., = (t), t ∈ [0, T] . This means that the control parameter = (t) tells us which dosage of the reverse transcriptase inhibitor should be given to the concrete patient at any time moment t for t ∈ [0, T] . Thus, personalized determination of this parameter for every individual is of vital importance for treatment of HIV.

Existence and Lyapunov Stability of Steady State
Let us now assume that the parameter in system (1) is constant. That is, (t) ≡ c ∈ (0, 1) . Setting du dt in (1) to zero and solving for u 1 , u 2 , u 3 and u 4 we can see that there are two possible steady states: an infected and an uninfected one [42].

Stability of Solutions
Let us define where , 1 , s etc. are the parameters of (1). Consider the set It can be proven [42] that if u(0) ∈ Γ( ) , then the solution trajectories of (1) will stay inside Γ( ) for all t ∈ Ω T .

Remark 1
It is not required that is constant. As long as ∈ M , we may allow (t) to vary with time.
For our parameters presented in Table 1, these bounds are quantitatively defined as (Table 3) shows upper limits for u 3 and u 4 for different values of . It can furthermore be proven that if and only if ≥ crit the uninfected state is globally asymptotically Lyapunov stable. On the other hand, if the steady state exists, then it is locally asymptotically Lyapunov stable whenever the following condition is satisfied [42] where We can calculate that, when is constant and the other parameter values are chosen as in Table 1, then the infected steady state is locally asymptotically stable for all values of such that is less than the critical value, crit ≈ 0.88 . Figure 3 illustrates this statement.
Thus, if is constant, and less than the critical value crit ≈ 0.88 , it suffices to know the solution of (1) at steady state to deduce .
Although it is often a reasonable assumption that the drug efficacy is constant for a given individual, viruses mutate readily, which can alter the efficacy of a RT-inhibitor. Thus, it is interesting to know how to determine (t) when it is not constant. So let us for the remainder of this note consider the case when (t) is not necessarily constant.

Well-Posedness of the Forward Problem
Let D = Ω T × Γ( ) be the bounded domain. Let functions u(t), f(t, u(t)) are defined as in (4). Further, let f(t, u(t)) be a continuous function for t in Ω T and Lipschitz continuous function for u(t) in D. Then f(t, u(t)) is clearly Lipschitz continuous on the compact set Γ( ) × Ω T . In other words, ∃L = const. ∶ ∀t ∈ Ω T , ∀u 1 (t), u 2 (t) ∈ D,

3
Thus, using the Picard-Lindelöf theorem (Theorem 2.2 in [43]) one can prove that, for given initial condition u(0), the model problem (1) has a unique solution. Furthermore, the solution depends continuously on data of the problem (1) in the following sense (Theorem 2.8 in [43]):

Inverse Problems and Ill-Posedness
Since the parameter identification problem can be considered as an inverse ill-posed problem it is clear that this problem is difficult to solve properly. In this section we show how parameter identification problem can be solved accurately via construction of proper Tikhonov regularization functional.
Let us consider the following problem: Let B 1 and B 2 be Banach spaces. Let G ⊆ B 1 be an open set in B 1 and F ∶ G → B 2 an operator. Let y ∈ B 2 be given, and suppose we want to find x ∈ G such that Problems of this sort, when you want to identify x in (17), given observations, y, are called inverse problems. A special class of inverse problems are called parameter identification problems (PIP), i.e. x is some parameter of a differential equation, and F(x) is the solution of the differential equation, with this parameter. (17) is said to be well-posed by Hadamard if it satisfies the following conditions [45,46]:

Definition 1 Problem
3. Stability: For each y such that a unique solution of (17) exists, the solution x = x(y) is a continuous function of y. 1 (17) is said to be ill-posed if it is not well-posed.

Definition 2 Problem
PIP and other inverse problems are often ill-posed. Ill-posedness means that it is difficult to solve (17) numerically, since measurement errors, or even errors induced by finiteprecision computer arithmetic, can have disastrous consequences. Let y * denote noiseless observations, > 0 be the noise level, and B [y * ] = {y ∶ ||y − y * || B 2 ≤ } . The solution to the slightly perturbed equation F(x) = y (with y ∈ B [y * ] ) could be entirely different from the solution to F(x) = y * . Perhaps a solution to F(x) = y does not even exist. No matter how small is. A generally ill-posed problem (17) can be well-posed if we consider the restriction of F in (17) to certain subsets of its domain. In this case is introduced the following definition.
(17) F(x) = y. 1 We will call a problem well-conditioned if it is well-posed and such that a small change in data results in a small change in the solution: even if x(y) is continuous it might still be very sensitive to changes in y.

3
Conditional well-posedness Let B 1 and B 2 be Banach spaces. Suppose G ⊂ B 1 is the closure of an open subset in B 1 . Let F ∶ G → B 2 be a continuous operator. Assume that y * ∈ F(G) is our ideal noiseless data, and pick a noise level > 0 . Suppose we want to solve where y ∈ B [y * ] . This problem is called conditionally well-posed on G if it satisfies the following conditions [45,46]: 1. Existence: It is a priori known 2 that there exists an ideal solution x * = x * (y * ) ∈ G for an ideal noiseless data y * . 2. Uniqueness:

Definition 3
The set G in Definition 3 is called the correctness set of the problem (18).
Continuity of the inverse operator F −1 can be guaranteed if the domain of F is compact. Hence, any compact set with nonempty interior such that F is one-to-one is a correctness set. This suggests a method to solve (18) by choosing a suitable correctness set, G, and then finding a x ∈ G such that ||F(x) − y || is as small as possible. The Tikhonov's theorem offers a such method.
Theorem 3 (Tikhonov [44]) Let B 1 and B 2 be Banach spaces, and U ⊂ B 1 a compact set. Let F ∶ U → B 2 be a continuous one-to-one operator and V = F(U) .
For a proof of this fundamental theorem see, for example [13,44].

Quasi-Solution
Let H 1 and H 2 be Hilbert spaces, 3 and assume that F ∶ G → H 2 is a continuous mapping defined on a compact correctness set, G ⊂ H 1 . Let > 0 and assume, as before, that we want to solve with y defined as before. We know that a solution exists for perfect data y * , but in general (19) has no solution, since y ∉ F(G) (implying that we are dealing with an ill-posed problem). Our goal in this, and the following subsection is to sketch how to construct a family of approximate solutions {x } in G that converges to x * as → 0 . Let us define Since F is continuous, it takes compact sets to compact sets, thus F(G) is compact in B 2 . And since F(G) is a compact subset of a Hilbert space, and therefore closed, a minimum of (20) exists (and if F(G) also happens to be convex, this minimum is unique). Any x ∈ G , unique or not, that minimizes J y in (20) is called a quasi-solution to (19).
Since the inverse mapping, F −1 , is continuous by the Theorem 3, and is defined on a compact metric space, it admits a modulus of continuity, F −1. 4 From Theorem 1.5 in [13] it follows that, given y ∈ B 2 , then for any quasi-solution x ∈ min x∈G J y (x) the following error estimate holds: where F −1 (z) is the modulus of continuity of the inverse operator F −1 . Thus x → x * as → 0 . Hence, we can take a sequence of quasi-solutions to be our desired family. However, sometimes the set of all plausible solutions to (19) does not form a compact set, and in these cases, F need not be continuous on the set of all plausible solutions-such problems are called essentially ill-posed. Thus, it is not obvious how to choose a suitable correctness set. And even if all the plausible solutions form a compact correctness set, the minimum of (20) may not be unique: there may be local minima or regions where the gradient of the functional is very small, where a minimization algorithm could get trapped. In the next subsection, we will discuss how a stable solution to essentially ill-posed problems could be obtained in practice.

The Tikhonov Functional
The Tikhonov functional makes sure that when minimizing (20), we will stay in the neighborhood of some point, x 0 , which is a priori known to be close to the true solution, x * . A general Tikhonov functional is given below The first term is essentially the same as in (20), the second term is the regularization term and ∶= ( ) is the regularization parameter. The regularization parameter can be chosen as, for instance, where ∈ (0, 1) , see details in [4].
In general, the Tikhonov functional (23) might not actually attain its infimum; we can only guarantee the existence of the minimizing sequence, {x k } . However, without loss of generality, we can assume that G is the closure of an open and bounded set containing the initial guess, x 0 , the (bounded) minimizing sequence, {x k } , and the exact solution, x * .
Hence, if we consider finite dimensional Hilbert spaces, the Tikhonov functional, defined on G, would have a minimum, since closed and bounded sets on finite dimensional Hilbert spaces are compact, and functionals defined on compact sets attain their infimum according to the Weierstrass' extreme value theorem. Suppose now that G is convex and that (23) is Fréchet differentiable, 5 with a Fréchet derivative that is uniformly bounded and Lipschitz continuous. Then one can prove, see [13,14], that for given noise level and regularization parameter (23) is locally strongly convex in a neighborhood of its minimum and that x * is also contained in this neighborhood if ||x * − x 0 || is small enough. Assume that we have a single noise and our goal is to minimize (23). Let ( ) is chosen as in (24), then it can be proven, see [28], that there exists a 0 such that in particular it follows that if (23) attains a minimum, any x k ∈ min x∈G J(x) would be a better approximation to x * than x 0 , if the noise level is small enough.
Thus, to sum up, under reasonable assumptions discussed above, a minimum of (23) exists and is a better approximation than the starting guess, x 0 . And under reasonable assumptions, and if the initial guess is good enough, there is only one unique point that is the global minimum, and we do not need to worry about local minima. These facts explain why Tikhonov regularization is so useful for solving ill-posed problems. To find the zero of the Fréchet derivative, one can use common minimization techniques, such as the conjugate gradient method (CGM) or the method of steepest descent. Obviously, the minimum of the Tikhonov functional will not be exactly the same as the quasi-solution if the noise level and regularization parameter are constants. On the other hand, by letting the regularization parameter decrease for each iteration of the minimization algorithm we will have a minimum of the Tikhonov functional that approaches the quasi-solution. In [4] it was suggested that can be updated as where p ∈ (0, 1] and k = 0, 1, 2, …. From what have been discussed, it is clear that a good first guess is essential for successful identification of the desired parameter. Of course, we do not in general have any idea at all what the solution to (17) might be, and therefore, in general, we need to devise some kind of globally convergent algorithm to solve ill-posed PIP.

The Parameter Identification Problem
To formulate the parameter identification problem we assume that all parameters in system (1) are known except the control parameter (t) which describes efficacy of the drug. The typical values of parameters {s, , k, 1 , , b, , c, N} in (1) are taken from [42] and they are described in Table 1.

3
Parameter Identification Problem (PIP). Assume that conditions (5) hold and parameters {s , , k, 1 , , b, , c, N} in system (1) are known. Assume further that the function (t) ∈ M is unknown inside the domain Ω T . The PIP is: determine (t) for t ∈ Ω T , under the condition that the virus population function g(t) is known Here, the function g(t) presents observations of the function u 4 (t) inside the observation Note, that we solve the PIP on the time interval [0, T] and assume that observations of g(t) can even be on the more narrow interval [T 1 , T 2 ] ⊂ [0, T] . "Numerical Results" section show that reconstruction of the parameter (t) is not very good on the time interval where observations are not available and thus, observations of the virus population function u 4 (t) should be taken as early as possible from the date when the virus started to be reproduced in the body of the host.

Optimization Method
Let H be a Hilbert space of functions defined in Ω T . To determine (t) , t ∈ [0, T] in PIP we construct the Tikhonov functional (23) in the following form: Here, the solution u 4 (t) of the system (1) with parameter (t) , g(t) is the observed virus population function, 0 is the initial guess for the parameter (t) and ∈ (0, 1) is the regularization parameter, z (t), ∈ (0, 1) is smoothness function which can be defined similarly to [10]. Our goal now is to minimize the Tikhonov functional (29) with respect to the function (t) ∈ H.
Let us introduce following spaces needed for further analysis where all functions are real valued.
To derive the Fréchet derivative of the Lagrangian (31) we assume that functions v = ( , u, ) can be varied independently of each other in the sense that Thus, we consider L(v +v) − L(v) , single out the linear part with respect to v of the obtained expression and neglect all nonlinear terms. The optimality condition (33) means that for all v ∈ U we have i.e., every component of (34) should be zero out. Thus, the optimality conditions (33) yields (32)

3
The equation (35) corresponds to the forward problem (1)-(2), the equation (36) -to the following adjoint problem which can be rewritten in the compact form as with The adjoint system should be solved backwards in time with already known solution u(t) to the forward problem (1)-(2) and a given measurement function g(t).
For the case when u and are exact solutions of the forward (1)-(2) and adjoint (39) problems, respectively, to the known function , we get from (31) that and thus the Fréchet derivative of the Tikhonov functional can be written as Using (37) in (42), we get the following expression for the Fréchet derivative of the Tikhonov functional Thus, to find the unknown parameter which minimizes the Tikhonov functional (29) we can use the following expression

Finite Element Discretization
For solution of (33) we will use the finite element discretization and consider a partition We define also the piecewise-constant time-mesh function such that For discretization of the state and adjoint problems we define the finite element spaces W u ⊂ H 1 u Ω T and W ⊂ H 1 Ω T for u and , respectively, as For the function (t) we also introduce the finite element space W ⊂ L 2 Ω T consisting of piecewise constant functions We use different finite element spaces since we are working in a finite dimensional space and all norms in finite dimensional spaces are equivalent. Next we denote U = W u × W × W such that U ⊂ U. Now the finite element method for (33) is: find v ∈ U such that Since the forward (1)-(2) and adjoint (36) problems are nonlinear their solutions can be found by Newton's method. For the discretization the variational formulation of the forward problem (1)-(2) for all ū ∈ H 1 u (Ω T ) is: The finite element method for (1)-(2) will be: find u k+1 ∈ H 1 u (Ω T ) such that for all ū ∈ H 1 u (Ω T ) Here, we can compute the Jacobian V � (ũ n ) via definition of V(ũ) in (51) as where I is the identity matrix, f � (ũ n ) is the Jacobian of f (the right hand side of the forward problem (1)) at ũ n and n is the iteration number in Newton's method. The explicit entries in the Jacobian f � (ũ n ) for system (1) are computed as We note that the finite element method (48) will work even in this case, see details in [24]. In a similar way the Newtons's method can be derived for the solution the adjoint problem (39). Since we solve the adjoint problem backwards in time starting from the known (T) = 0 , we discretize time derivative as for the already known k+1 , and write the variational formulation of the adjoint problem for all ̄∈ H 1 (Ω T ) as The finite element method for (39) will be: find k ∈ H 1 (Ω T ) such that for all ̄∈ H 1 (Ω T ) where I is the identity matrix, f � (̃n) is the Jacobian of f (the right hand side of the adjoint problem (39)) at ̃n , and n is the iteration number in Newton's method. The explicit entries in the Jacobian f � (̃n) for the adjoint system (39) are given by Taking into account values of parameters given in Table 1, we observe that detf � (̃n) ≠ 0 as well as detf � (̃n) ≠ 0 . Thus, schemes (53), (59) will converge given the appropriate starting values ũ 1 and ̃1 , correspondingly. For study of convergence of iterative methods we refer to [4].

A Posteriori Error Estimates
We consider the function ∈ C(Ω T ) as a minimizer of the Lagrangian (31), and ∈ W its finite element approximation. Let us assume that we know good approximation to the exact solution * ∈ C(Ω T ) . Let g * (t) be the exact data and the function g (t) represents the error level in these data. We assume that measurements g(t) in (28) are given with some noise level (small) such that Accordingly [14] we assume that and (57) (61) = ( ) = 2 , ∈ (0, 1∕4), ∈ (0, 1)

3
where * is the exact solution of PIP with the exact data g * (t) . Let Assume that for all ∈ V 1 ( * ) the operator has the Fréchet derivative F � ( ) which is bounded and Lipschitz continuous in V 1 ( * ) for D 1 , D 2 = const. > 0

An A Posteriori Error Estimate for the Tikhonov Functional
In the Theorem 1 we derive an a posteriori error estimate for the error in the Tikhonov functional (29) on the finite element time partition J .

Theorem 1
We assume that there exists minimizer ∈ C(Ω T ) of the functional J( ) defined by (29). We assume also that there exists finite element approximation of a minimizer ,

A Posteriori Error Estimate of the Minimizer on Refined Meshes
Theorems 2 and 3 present two a posteriori error estimates for a minimizer of the functional (29). Proof of the next theorem follows from the proof of Theorem 5.1 of [29].
Theorem 2 Let ∈ W be a finite element approximation on the finite element mesh J of the minimizer ∈ L 2 (Ω T ) of the functional (29) with the mesh function (t) . Then there exists a Lipschitz constant D = const. > 0 defined by

3 and interpolation constant C I independent on such that the following a posteriori error estimate for the minimizer holds true
Proof Let be the minimizer of the Tikhonov functional (29). The existence and uniqueness of this minimizer is guaranteed by conditions (62) and follows from Theorem 1.9.1.2 of [13]. By this theorem, the functional (29) is strongly convex on the space L 2 (Ω T ) with the strong convexity constant . This implies that Here, J � ( ), J � ( ) are the Fréchet derivatives of the functional (29) given by (43) for respective .
Since is the minimizer of the Tikhonov functional (29) then Using the splitting where I is an interpolant of , together with the Galerkin orthogonality principle for all , I ∈ W in (79) we obtain We can estimate the right hand side of (82) using (77) as

Substituting above equation into (82) we obtain
Using the interpolation property we obtain a posteriori error estimate for the regularized solution with the interpolation constant C I :

3
We can estimate || || H 1 (Ω T ) in (85) similar to (74). Substituting this estimate into the right hand side of (85) we get In the case when ∈ W terms with jumps in time [ ] disappear and we have a posteriori error estimate Using the interpolation property (84) and further the estimate (74) we obtain following a posteriori error estimate for the regularized solution with the interpolation constant C I :

3
and since ∈ W the above estimate reduces to ◻

Algorithms for Solution of PIP
Here we present two algorithms for solution of PIP: • CGA-usual conjugate gradient algorithm on a coarse time partition, • ACGA-time-adaptive conjugate gradient algorithm which minimized the Tikhonov functional (29) on a locally refined meshes in time.
We denote the nodal value of the gradient at the observation points {t i } by G m (t i ) and compute it accordingly to (43) as The approximate computed solutions u 2 m and 1,3 m are obtained computationally by Newton's method with ∶= m . A sequence { m } m=1,…,M of approximations to is computed as follows with and where d 0 (t i ) = −G 0 (t i ) and G m (t i ) is the gradient vector which is computed by (91) in time moments t i . In (92) the parameter r m is the step-size in the gradient update at the iteration m which is computed as

Numerical Results
In this section we present several numerical results which show performance and effectiveness of the time-adaptive reconstruction of unknown parameter (t), t ∈ [0, T] in PIP using ACGA algorithm. Numerical tests are performed in Matlab R2019b using the developed code for solution of the studied problem available for download at [31]. Numerical results of reconstruction of function (t) using usual conjugate gradient Algorithm 1 on the nonrefined timemeshes are presented in [25]. We note that observations of all u i , i = 1, 2, 3, 4 functions in system (1) were used in [25]. The goal of numerical tests of this note is to determine the unknown function (t) from observation of the virus population function u 4 (t) in (1) on the interval In all numerical tests assumed that parameter (t) satisfy conditions (5) and is unknown in the system (1), but all other parameters {s , , k, 1 , , b, , c, N} of this system are known and their values are chosen as in the Table 1. The observation interval [T 1 , T 2 ] is such that T 2 = T = 300 , but T 1 is taken differently in different tests since observations of the virus population function u 4 (t) can be taken after the first 3 − 9 weeks since the virus started to be reproduced in the body of host.
For generation of data u 4 (t) = g(t) the problem (1)-(2) was solved numerically with exact values of the test model function (t) . For solution of problem (1)-(2) was used Newton's method presented in "The Parameter Identification Problem". Next, the random noise was added to the observed solution u 4 (t) as where ∈ [0, 1] is nose level and ∈ [−1, 1] is random number.
In Algorithms 1, 2 it is of vital importance to take initial guess 0 such that it satisfy condition (62) which means that 0 is located in the close neighborhood of the exact solution. This condition is fulfilled in our PIP since we can compute explicitly values of the parameter (t) on the initial non-refined time mesh using, for example, the third equation of system (1) as We used following discretised version of this equation to get initial guess 0 (95) u 4 (t) = u 4 (t)(1 + ),

3
Here, u 3 k+1 , u 3 k , u 2 k are known computed approximations of functions u 3 , u 2 at time iterations k + 1 and k, respectively. In our computations we take values of u 3 k+1 , u 3 k , u 2 k using solution of the problem (1)- (2) with exact values of the test model function (t) and then adding the noise as where ∈ [0, 1] is nose level and ∈ [−1, 1] is random number. We note that denominator in (97) is not approaching zero because = 0.4 and u 2 (t) > 0 ∀t ∈ [0, T] . To get reasonable approximation 0 for the initial guess 0 in Algorithm 2 we assume that noisy functions u 3 , u 2 are known on the initial non-refined mesh, apply (97) and then use polynomial fitting to obtained noisy data 0 in order smooth them. Finally, the condition (5) was applied for the computed 0 in order to ensure that 0 (t) belongs to the set of admissible parameters M . Second order discretization of the first time derivative in (97) is also possible. We note that numerical differentiation of noisy data is an ill-posed problem and it is discussed in detail in Section 4 of [25].
All tests are performed with tolerance Θ = 10 −7 in ACGA algorithm and 1 = 0.1 in (8.94). The value of 1 is chosen such that it allows local refinements and avoids refinement of the very large time region in the time mesh. All tests are performed for different T 1 = 25, 50, 100 for the time interval [T 1 , T 2 ] = [T 1 , 300] which corresponds to the fact that HIV virus can be detected in the first 3-9 weeks after infection.
Relative errors in the reconstructed parameters (t) presented in the Tables are measured in L 2 -norm and are computed as (98) u i (t) = u i (t)(1 + ), i = 2, 3,

3
Complete description of all numerical tests with reconstruction results are presented in the recent work [11].
In this test we present the reconstruction results of the smooth model function  Table 1. Right figures of Fig. 4 show the reconstruction results of the function (100) for a noise level = 10% in the data u 4 (t) for starting observed time T 1 = 50 . Right figures of Fig. 5 shows the reconstruction results of this function for a noise level = 40% in the data u 4 (t) and for T 1 = 100. Table 1 and Figs. 4 and 5 confirm that with local time-mesh refinements the reconstruction of the drug efficacy function is significantly improved compared to the reconstruction of obtained on the initial non-refined time-mesh.

3
In this test we present numerical reconstruction results of the constant model function (t) = 0.7 from noisy observations of the virus population function u 4 (t) at the observation interval [T 1 , T 2 ] . We again took T 1 = 25, 50, 100 , but number of observation points were 20 at the time interval [T 1 , T 2 ] = [T 1 , 300] . We generate initial time partition J with equidistant time step = 300∕19 . The reconstruction results of the model function (t) = 0.7 for noise levels = 5%, 10%, 20%, 40% in data u 4 (t) are presented in Table 2. Right figures of Figs. 6, 7 show the reconstruction results of the function (t) = 0.7 for noise level = 40% in the data u 4 (t) for T 1 = 50 and T 1 = 100 , respectively. We again observe from the results of Table 2 and Figs. 6 and 7 that with local timemesh refinements the reconstruction of the drug efficacy is significantly improved compared to the reconstruction of obtained on the initial non-refined time-mesh even if we add large noise = 40% to the observed data u 4 (t).

Conclusion
The finite element time-adaptive optimization method for determination of the drug efficacy in a mathematical model of HIV infection with drug therapy is presented. Timeadaptive optimization means that first the time-dependent drug efficacy is determined at a known coarse time partition using several known values of observed function (usually, we used 15-20 observations). Then the time-mesh is locally refined at points where the computed residual |R( )| attains its maximal values and the drug efficacy is computed The proposed new time-adaptive method can eventually be used by clinicians to determine the drug-response for each treated individual. The exact knowledge of the personal drug efficacy can aid in the determination of the most suitable drug as well as the most optimal dose for each person, in the long run resulting in a personalized treatment with maximum efficacy and minimum adverse drug reactions.
The proposed time-adaptive method can be adopted to solve multi-parameter identification problems for a bread class of problems stated by the system of ODE.