Determination and biological application of a time dependent thermal parameter and sensitivity analysis for a conduction problem with superficial evaporation
Natalia N. Salva Maria S. Herrera Andrea Monti Hughes Claudio Padra Gustavo A. Santa Cruz
physics.bio-ph


Horizontal patterns from finite speed directional quenching
     Rafael Monteiro

University of Minnesota, School of Mathematics,  
206 Church St. S.E., Minneapolis, MN 55455, USA
    December 30, 2023
=======================================================================================================================




A boundary value problem, which could represent a transcendent temperature conduction problem with evaporation in a part of the boundary, was studied to determine unknown thermophysical parameters, which can be constants or time dependent functions. The goal of this paper was elucidate which parameters may be determined using only the measured superficial temperature in part of the boundary of the domain. We formulated a nonlinear inverse problem to determine the unknown parameters and a sensitivity analysis was also performed. In particular, we introduced a new way of computing a sensitivity analysis of a parameter which is variable in time. We applied the proposed method to model tissue temperature changes under transient conditions in a biological problem: the hamster cheek pouch. In this case, the time dependent unknown parameter can be associated to the loss of heat due to water evaporation at the superficial layer of the pouch. Finally, we performed the sensitivity analysis to determine the most 
sensible parameters to variations of the superficial experimental data in the hamster cheek pouch.




 


Keywords:  Boundary value problem, Determination of parameters, Sensitivity analysis, Hamster cheek pouch.










§ INTRODUCTION


Diffusion problems have been studied in many different areas, such as heat conduction  <cit.>, sediment transport in a river  <cit.>, the study of the weather <cit.>, or in drying food technology  <cit.>. The heat transport in a certain object may be modeled through a convection-diffusion differential equation. The boundary conditions usually considered are: a constant temperature (Dirichlet condition), a fixed heat flux (Neumann condition) or a mixed condition, where the heat flux depends on the temperature at the boundary.

In <cit.> some thermal inverse problems were considered in a mathematical tumor model. They estimated simultaneously unknown thermophysical and geometrical parameters, through an evolutionary algorithm. All the parameters were constants and no evaporation was considered. 
In <cit.> a space-dependent convection parameter was determined through a nonlinear least squares technique, using several temperature measurements at different times in different locations of the domain.
In <cit.>, we observed that exophytic tumors and irradiated tissues, principally those with ulcers, exhibit mass moisture transfer in the tissue-air interphase.

In this work,  we resolved a partial differential problem considering a time-dependent extra term in the mixed condition, which can be associated with a loss of heat due to superficial evaporation. Having measurements of temperature in a section of a particular domain during time, a nonlinear inverse problem was formulated to determine some thermophysical parameters such as the superficial heat loss or the heat conductivity. A sensitivity analysis was also  performed. Finally the parameter determination and the sensitivity analysis were validated using our previous study on tissue temperature responses in the hamster cheek pouch <cit.>. 




§ DESCRIPTION OF THE CONDUCTION PROBLEM

Let us consider a one dimensional spatial domain Ω=[0,X_max], and the following boundary value problem:

 

    c_1 ∂ T ∂ t (x,t)=c_2  ∂^2 T ∂ x^2 (x,t) +c_3 T(x,t)+c_4,    ∀   (x,t) ∈Ω× [0,t_max] 
        c_5 ∂ T∂ x(x,t)= c_6 T(x,t)+c_7+f(t),    ∀   (x,t) ∈Γ_u × [0,t_max] 
        T(x,t)=g(t),    ∀  (x,t) ∈Γ_b × [0,t_max] 
        T(x,0) =  h(x)    ∀   x ∈Ω





where  c_i are real constants, f,g:[0,t_max] →ℝ and h:Ω→ℝ are real functions, and Γ_u={X_max},  Γ_b={0}[We chose one spatial dimension for simplicity. This procedure can be easily extended to more dimensions without loss of generality.].

In general, the parameters c_i,   i=1..7, f,g and h are fundamental in using this type of differential equation problem with boundary conditions given by Eqs.  (<ref>), (<ref>) and (<ref>). Note that Eq. (<ref>) includes a time-dependent parameter, f(t), which is defined only in the boundary Γ_u. In particular, c_i  are often approximated, and sometimes their values are just guessed. If there is an experimental measure of the temperature T(x,t)  on the boundary Γ_u, an inverse problem can be formulated defining the following cost function:

    J(T)=∫_0^t_max∫_Γ_uT(x,t)-T^* (x,t)^2 dx dt  ,

where T^*(x,t) is the measured superficial temperature. Therefore, some parameters can be determined minimizing the cost function J. 

The goal of this paper is to find which parameters c_i  may be determined using only the measured superficial temperature. First we proposed a method to identify the function f=f(t) such that the solution T of the problem (<ref>) minimizes the cost function (Section <ref>). Then, we introduced a new way of computing a sensitivity analysis of a variable parameter (Section <ref>). As an example, in Section <ref>, we applied this problem to model tissue temperature changes, under transient conditions, in the hamster cheek pouch. In this biological problem, the function f(t) can be associated to the loss of heat due to water evaporation at the superficial layer of the cheek pouch.  Finally, we performed the sensitivity analysis to determine the most sensible parameters to variations of the superficial experimental data (Section <ref>).   




§ DETERMINATION OF THE FUNCTION F
 

The goal of this section is to determine the function f:[0,t_max] →ℝ, assuming that the rest of the parameters are known constant.  
To solve the direct problem (<ref>) we used the Finite Element Method (FEM), meshed the spatial-time[The mesh used is a triangular uniform mesh, where the nodes are obtained by the Cartesian product of a discretization in time and a discretization in space. Using polynomials of degree one and this particular mesh makes the FEM equivalent to a finite difference scheme.] domain Ω× [0,t_max] and obtained the triangulation 𝔗={𝒯_i,i=1...M}. 
We define the following finite element spaces:


    S_𝔗={ v ∈ H^1(Ω× [0,t_max])/ v ∈𝒫_1(𝒯), ∀ 𝒯∈𝔗}, 𝒫_1



    V={ v ∈ S_𝔗/ v=g,   Γ_b × [0,t_max]  ∧  v=h Ω×{0}} ,


    V_0={ v ∈ S_𝔗/ v=0,   Γ_b × [0,t_max]  ∧  v=0 Ω×{0}} .



The weak formulation of  problem (<ref>) is the following:


    (VP):  T ∈ V/    a(T,η)=l(η),  ∀η∈ V_0


where,


    a(u,v)=Ω× [0,t_max]∫ c_1  ∂ u/∂ t  v+ c_2  ∂ u/∂ x∂ v/∂ x - c_3   u v   dx dt   -  Γ_u × [0,t_max]∫c_2   c_6/c_5 u v   dx dt 
     l(v)=Ω× [0,t_max]∫ c_4   v   dx dt  + Γ_u × [0,t_max]∫c_2/c_5  (c_7+f(t)) v  dx dt  .


Let {(x_i,t_i), i=1,..,N} be the nodes of the triangulation 𝔗, and  I ⊂{1,...,M} be the subset of index whose nodes are in the boundary  Γ_u ×  [0,t_max]. We define the function f as follows,


    f(t)=∑_i ∈ Iβ_i η_i(X_max,t),


where {η_i }_i=1..N is the base of the space V. 

The inverse problem of determining f consists in obtaining the values of {β_i }_i ∈ I such that the solution T of the weak formulation problem is close to the measured temperatures. Let 
T be a solution of (VP). The cost function can be approximated by:


    J(T)=Γ_u × [0,t_max]∫ (T-T^*)^2   dx dt ≈∑_j=1^N_t (T(X_max,t_j)-T^*(t_j))^2 .
 

We defined the following optimization problem: 


 
   Find B={β_i }_i ∈ I such that minimizes the cost function j({β_i }_i ∈ I)=J(T_B), where T_B is a solution of (VP) using the values of B={β_i}_i ∈ I to define f(t). 

 


We used the Lagrange Method to determine the variation of j by each β_i, and obtaining a direction of descent. Suppose that for r ∈ I-{i} the values of β_r are fixed, and we vary only the parameter β_i. We define the following function:

    ℒ(u,w,β_i)=J(u)+a(u,w)-l_β_i(w)

where a(·,·) y l(·) is defined in (<ref>) and (<ref>), and  the subscript β_i means that we are using in f(t) the fixed values of {β_r}_r ∈ I-{i} and the value of β_i (which may vary).

Note that if u_β_i is a solution of (VP) then  ∀ w ∈ V_0:  ℒ(u_β_i,w,β_i)=J(u_β_i)=j(β_i), and therefore their derivatives are equal.

Using the chain rule we obtained the derivative of the Lagrangian ℒ respect to each β_i:

    δℒ(u_β_i,w_β_i,β_i)/δβ_i= Γ_u × [0,t_max]∫η_i w_β_i  dx dt

where w_β_i∈ V(Ω) is the adjoint state, which is defined as the unique solution to  the following variational problem:

    (VPA) {[ a(v,w_β_i)=-2 Γ_u × [0,t_max]∫ (u_β_i-T^*) v   dx dt,                                               ∀ v ∈ V;                                               w_β_i=0                                       Γ_b × [0,t_max] ].


To find the optimal values of {β_i} we used the gradient descent method.
This method is based on the observation that if the function j(β ) is defined and differentiable in a neighborhood of a point β^0={β_1^0, ⋯, β_N^0}, then j(β ) decreases fastest if one goes from β^0 in the direction of the negative gradient of j at β^0 , e.g. -∇ j(β^0 )={∂ j(β^0 )/∂β_1, ⋯, ∂ j(β^0 )/∂β_N}. It follows that for γ small enough, the value of j in β^1 =β^0 -γ∇ j(β^0 ) is smaller than j(β^0 ).  We get a sequence β^0,β^1,β^2,… such that  β^n+1=β^n-γ _n∇ 
j(β^n) and  j(β^n+1)≤ j(β^n), for all n ∈ℕ_0. The convergence of this sequence to a local minimum of j depends on the properties of j, (for example, j convex and ∇ j Lipschitz).

The value of γ_n >0 is different in every step. For each n, we searched for the greatest γ_n ≥ 0 such that it minimizes j(β^n-γ∇ j(β^n)). This is achieved combining a dichotomy method and a method that approximates j by a parabola, for more information about this procedure, see page 51 of  <cit.>. 




§ SENSITIVITY ANALYSIS


In <cit.> Blackwell and Dowding analyzed the usage of sensitivity parameters in connection with the estimation of thermal properties in the heat conduction equation. They stated that although parametric investigations may be done, sensitivity parameters are rarely computed (see also <cit.>). Sensitivity parameters help to understand the parametric dependence of an
experiment and shape our experience and intuition for future cases .

In inverse problems, the sensitivity parameter is the partial derivative of the output function (in our case the temperature) with respect to a parameter being determined, which is ∂ T/ ∂ p for a parameter p. Due to the general interest on the comparison of magnitudes for different parameters, a scaled (sometimes called “modified”) sensitivity parameter is used:

    S_p := p ∂ T/∂ p .


Note that equation (<ref>) has units of temperature for all parameters, therefore magnitudes for various parameters can be directly compared. 
Small sensitivity parameters, or general insensitivity, are beneficial when the parameters are not well quantified, such as materials with no characterized thermal properties. Then the parameter is not influential in the thermal response. To estimate a parameter, however, the measured response has to be sensible to that parameter. In this case, the scaled sensitivity parameters are desired to be larger in magnitude (compared to the representative temperature) and linearly independent (having different shapes). The more sensitive the temperature is, the more valuable the temperature measurements are. In a similar way, the estimation of multiple parameters requires that the sensitivity, or the effect on temperature of each parameter, is different or independent of one another for each parameter. If two parameters have similar effects on temperature, their individual influence is difficult to distinguish.

In this work we analyzed the sensitivity of the superficial temperature, for which we defined the modified  sensitivity parameter as follows:

    S_p (t) := lim_Δ→ 0 p  T_p+Δ(X_max,t) -T_p(X_max,t)/Δ ,


where X_max indicates the superficial boundary, and T_α is the solution of the differential problem (<ref>), using the value α in the parameter p.

We used the method of finite differences <cit.> to determine the sensitivity parameters. First we solved the direct problem  (<ref>) using the values of the parameters  p=(p_1,p_2,...,p_i,..,p_n), and then we solved the direct problem again but with the following values: p_=(p_1,p_2,...,p_i+ p_i,..,p_n). Finally we approximated the sensitivity parameters as follows:


    S_p_i (t) ≈  p_i  T_p_ (X_max,t) -T_p (X_max,t)/ p_i .


We remark that as every partial derivative, it is dependent not only on time, but also depends strongly in the values assumed by the parameters p_i,   i=1...n.

This previous analysis works if the parameters are constant. In the case of f which is variable in time, we introduced a new way of defining the sensitivity parameter of a variable parameter:



    {[ p=(p_1,p_2,...,f,..,p_n),    p_=(p_1,p_2,...,(1+).f,..,p_n);                     S_f ( t)≈ T_p_ (X_max,t) -T_p (X_max,t) ].


where (1+).f represents the product of a real number and a function.


Note that in order to compute S_p_i we need to work in a concrete problem, and therefore we will perform the sensitivity analysis for the biological application, where f will be determined (Section <ref>).









§ BIOLOGICAL APPLICATION

It was previously demonstrated that the hamster cheek pouch is useful for the  study of tissue temperature affected by tissue superficial humidity <cit.>. The hamster cheek pouch is widely used as a model of oral cancer and mucositis, an adverse side effect induced by several cancer therapies <cit.>. Our group is focused on the study of BNCT (Boron Neutron Capture Therapy), a binary treatment modality that can selectively target neoplastic tissue <cit.>. Particularly, we study BNCT therapeutic effect on tumors and BNCT induced mucositis in the hamster cheek pouch with a non-invasive complementary method called Dynamic Infrared Imaging (DIRI). This method is based on the observation of temperature changes under transient conditions associated with mass moisture transfer in the tissue-air interface of the pouch. 
In our previous studies, we described different temperature changes for normal and tumor tissue, and also for non-irradiated and irradiated pouches <cit.>. However, the study of the mass moisture transfer as a function of time was not quantified.   




 §.§ Dynamic Infrared imaging (DIRI) studies in the hamster cheek pouch

Dynamic Infrared imaging (DIRI) is based on the acquisition of thermal images during transient processes, caused by sudden and sustained changes in surface temperature due to the application of a thermal stimulus (provocation test) that forces the neurovascular system to respond in order to maintain local and body temperature within normal parameters  <cit.>. Other authors followed this concept, including our group <cit.>, in different clinical research studies using thermography <cit.>.  DIRI provides a non-invasively supplementary in vivo information potentially useful to characterize normal and pathological tissues and their response to cancer therapy. 


The biological model, experimental setup and procedures of the DIRI studies can be found in <cit.>. Briefly, a total of 61 hamsters were examined under DIRI protocol. Following an acclimatization period in the room, the animals were anesthetized and the pouch was everted using a plastic pipette held by hand. Thermal responses were measured using a FLIR T420 infrared camera, before, during and after the provocation test, namely, Transient Equilibrium Phase (TEP), Provocation Test (PT) and Recovery Phase (RP), respectively. The PT consisted of a mild air current applied at ambient temperature for about 120 seconds. The purpose of an air stimulus is to eliminate the initial moisture condition of the tissue so that, in the RP, we can focus on its thermal behavior and evaporation process that occur in response of the PT. In TEP and RP no air was applied, leaving the pouch exposed to ambient conditions without perturbations during approximately 280 and 400 seconds, respectively. 

Figure <ref>(a) shows the normal hamster cheek pouch tissue. The measured temperature values during time were extracted from the thermal image (Figure <ref>(b)) and averaged in a user-defined region of interest (ROI) used to delineate the  normal tissue. 






 §.§ Determination of the superficial temperature and superficial heat loss

In this study, the differential equation model (<ref>) was applied to model the hamster cheek pouch temperature, imposing a loss of heat in the superficial tissue, due to the water evaporation. For this simulation, we used the experimental temperature data as a function of time in a given ROI  <cit.>.

The heat transfer modeling in organs  has proposed numerous equations, studied by Pennes since 1948 <cit.>. He suggested that the heat transfer rate between blood and tissue was proportional to the product of the volumetric perfusion rate and the difference between the arterial blood temperature and the local tissue temperature. Therefore the temperature of a tissue depends on the rate of blood perfusion, the metabolic activity and the heat conduction between the tissue and the environment. 


Taking into account the suggestion made by Pennes <cit.>, Eq. (<ref>) of problem (<ref>) takes the form:

















    ρ c ∂ T (x,t)/∂ t =k Δ_x T +ω_b ρ_b c_b (T_b-T)+q_m,    (x,t) ∈Ω× [0,t_max]


where  ρ (ρ_b) represents the tissue (blood) density, c (c_b) is the tissue (blood) specific heat, k is the thermal conductivity, ω_b is the blood perfusion coefficient, q_m is the metabolic heat source and T_b is the constant blood temperature. The boundary condition Eqs. (<ref>), (<ref>) and  (<ref>) are:





    {[                             -k ∂ T∂ x=h (T-T_amb)+L(t),             Γ_u × [0,t_max];                                                  T=T_D,             Γ_b × [0,t_max];                                        T(x,0) =  F_0(x)             ∀ t ∈ [0,t_max];                             ].


Where,  Γ_b={0} represents the inner boundary and let Γ_u={X_max}  represents the superficial boundary. Here, f(t) = L(t) (W/m^2) represents the superficial heat loss due to water evaporation, h is the heat transfer coefficient  between the tissue and the air and T_amb is the ambient temperature. We used the following initial temperature, that assures the continuity of the temperature values between the two boundaries  Γ_b and Γ_u:


    F_0(x)=(T_0-T_b)X_max x+T_b ,

 where T_0 is the initial superficial temperature measured. We used a linear function for simplicity.

Since thermal responses in the hamster cheek pouch were assessed before, during and after the application of a thermal stimulus (TEP, PT and  RP) <cit.>, these three different phases were modeled considering that the heat transfer coefficient h assumed different values in each stage:


    h=h(t)= 
    		     h_1     0<t≤ 120
    
    		     h_2     120<t ≤ 200
    
    		     h_3      200<t ≤ 400


The values used for h_1,h_2 and h_3 and others parameters are shown in Table <ref>.







To our knowledge, there are no published data related to the thermophysical properties (such as thermal conductivity, diusivity, etc.) of the hamster cheek pouch. In <cit.>, Poppendiek et al. suggested that tissues may be considered accurately for thermal analysis as being composed of water, protein and fat. Thus, for the thermal properties needed to compute Eq. (<ref>) in our biological application, we only considered water and protein to calculate the thermal properties, shown in Table <ref>. The other parameters not related to the composition of the tissue were obtained from <cit.>.
























 

 
 
 







Figure <ref> shows the normal hamster cheek pouch experimental and calculated thermal response as a function of time. In Figure <ref>(a) it can be seen the good approximation between the calculated  superficial temperature and the experimental data.  In particular, Figure <ref>(b) depicts the convective coefficient obtained using  an initial constant coefficient L(t) = 1000 W/m^2, and 15 steps of the gradient descent method described in Section <ref>.

In Figure <ref>(b),  we observed that at the beginning of each stage the superficial heat loss had oscillations. These oscillations can be seen also in the derivative of the functional (see Eq. (<ref>)), and therefore it seems that they are intrinsic of this differential problem. Besides, during the provocation test, the superficial heat loss decreases, due to the air stimulus that helps to eliminate the tissue superficial moisture.















 §.§ Sensitivity analysis
 
Table <ref> and  <ref> show those parameters used for the sensitivity analysis. We also used the parameter L(t) in Figure <ref>(b), obtained in the previous section. The initial temperature is in Eq. (<ref>). Table <ref> summarizes the physiological parameters considered, showing their units, and their position in the differential problem[The position in the differential problem, which was taken into account when the sensitivity was computed, changing the parameter p to p + Δ whenever it appeared.], necessary for the sensitivity analysis.  Note that some parameters, such as ω_b and k, appeared in more than one place in the differential problem. Figure <ref> shows different parameters sensitivities, grouped together depending on their order of the sensitivity. A summary of the order of the sensitivity parameters is shown in Table <ref>.





































The first important observation is that the most sensible parameter is the ambient temperature (T_amb) (Figure <ref>(a)). This could be explained by the thinness of the tissue and the Dirichlet condition in the inner surface. Therefore, this parameter should be determined with the smallest error possible.

Secondly, we examined the linear dependence between parameters, studied in <cit.>, which could be established by analyzing the shapes of the sensitivity parameters. If two parameters are linearly dependent, this means that there are infinite possible solutions, which implies that there will be infinite local minimums and therefore the minimization problem will not converge.

Figure <ref>(b) shows that L and k are linearly dependent parameters, because their sensitivity are symmetric with respect to the horizontal axis. Therefore, although they have similar sensitivities, they should not be simultaneously determined. Figure <ref>(d) shows that ω_b and T_b are linearly dependent parameters, because their sensitivities have the same behavior (increasing functions). Moreover, these parameters are linearly dependent with h_1,h_2 and h_3. Therefore to determine a parameter, of this order of sensitivity,  we should choose between ω_b, T_b and {h_1,h_2,h_3}.

Finally the least sensible parameter is the metabolic heat q_m, which may be justified by observing that the total heat source is Q=ω_b ρ_b c_b T_b+q_m, and the predominant term is ω_b ρ_b c_b T_b.







§ CONCLUSIONS


We proposed a method for the determination of time dependent parameters using measured superficial temperatures in a conduction problem with evaporation, and a new way of computing a sensitivity analysis of a variable parameter. We applied this method successfully to a biological problem, modeling tissue temperature changes, under transient conditions, in the hamster cheek pouch. In this study we found a good approximation between the calculated superficial temperature and the experimental data. We performed a sensitivity analysis, which should be done whenever parameters are simultaneously determined. 
Based on temperature measurements and having calculated the loss of heat due to water evaporation at the superficial layer of the pouch, the sensitivity analysis determined which of the studied parameters were the most sensible to variations of the superficial experimental data. We found that ambient temperature should be measured with the smallest error, because it is the most sensible parameter in this problem. Moreover, we noted that the linear dependence between the conductivity and the superficial heat loss is not intuitive,  in contrast with the dependence between the blood perfusion and the blood temperature (an increase in either of them would result in a raise in the superficial temperature). Therefore, we conclude that a choice has to be made between determining the conductivity or the superficial heat loss. 

In previous studies, we observed that tumors and particularly a precancerous tissue bearing ulcers after BNCT had high superficial humidity <cit.>.  Thus, in future studies, the proposed mathematical model will be extended to explore the mass moisture transfer as a function of time in tumors and precancerous tissue in the hamster cheek pouch. 





 §.§ Acknowledgments

This paper has been partially supported by the BNCT project of CNEA, and by CONICET.


ieeetr


