Variational finite element method to study the absorption rate of drug at various compartments through transdermal drug delivery system

Peer review under responsibility of Alexandria University Faculty of Medicine.

Available online 25 November 2014

Abstract

Background

The delivery of drugs through dermal layers known as transdermal drug delivery is one of the important contributions in medical practice. It represents an alternative to oral drug delivery and is poised to provide a substitute to hypodermic injections too. The diffusion of drug from the source to the target site is therefore an important issue to address.

Aim

This paper is an attempt to establish a mathematical model for the diffusion of drugs through the transdermal drug delivery system. The model identifies the pattern of drug diffusion in human body and its effective absorption rates at various compartments of skin and sub-cutaneous tissues.

Methods

The finite element method has been used to obtain the solution of the mass diffusion equation with appropriate boundary conditions. The tissue absorption rate of drug has been taken as the decreasing function of drug concentration from the skin surface towards the target site. The concentration at nodal points has been calculated which in turn determines the drug absorption at various layers.

Conclusion

The drug concentration at the nodal points of different dermal layers has been computed and the graphs were plotted between drug concentration and thickness of dermal layers using MATLAB software. It has been observed that due to dense network of connective tissues in dermal and sub-dermal parts, the drug absorption is maximum as compared to cutaneous tissues.

Keywords:

• Transdermal drug delivery
• Diffusion equation
• Finite element method

1 Introduction

In medical sciences, various attempts have been made to search new drugs and to find the improved ways of drug delivery in the body. Oral drug mechanism delivery like pills has been considered as one of the most appropriate methods of drug administration for decades worldwide. There are many drugs which cannot be taken orally so they are injected into the body by making use of hypodermic needles. However, the hypodermic injections have many disadvantages such as the presence of pain, the appearance of infections and the requirement for medical expertise to complete the process.¹ For the administration of drug and its absorption to target the desired compartment of the biological tissue; it is one of the challenging problems to control the adverse effects of the drug delivery through in-vivo tissues. Thus, it is desirable to address this issue by using an appropriate method of drug administration through dermal regions of the human body with minimal side effects. Transdermal drug delivery (TDD) is an alternate route of drug administration to pills and injections. This method operates by delivery of drugs into the human body across the skin and subcutaneous layers. TDD system has the ability to overcome the problems associated with the traditional delivery of drugs.

Transdermal drug delivery system has been in existence since last thirty years. Transdermal therapeutic systems are defined as self contained, continuous dosage forms which when applied to the skin deliver the drug through it at a controlled rate to the systemic circulation. TDD is a feasible administration route for potent and low molecular weight therapeutic agents. The choice of therapeutic agent is determined by a number of factors which includes the physio-co-chemical properties of the drug and its interactions with the membrane. This is due to the barrier function of the skin represented by the outer layer- the stratum corneum, which generally allows diffusion of only small molecular weight solutes. To overcome this diffusion limitation at the stratum corneum, various methods have been developed to enhance drug delivery across stratum corneum which includes chemical enhancers or physical enhancer technique e.g. iontophoresis^† and ultrasound.

The transport of drug through the device and skin follows the Fickian diffusion process by means of a simple homogeneous membrane. Kalia and Guy² have established the mathematical transport models describing the release of drug from various delivery systems and formulations based on the solution of Fick’s second law. In these models, skin was simply treated as a boundary condition and the main attention was focussed on the drug transport, also analytical expressions were found for the amount of drug released from the device. Qaliaf et al.³ developed a mathematical model and carried numerical simulations to describe the pharmacokinetics of drug penetration into the skin from microneedle array. Missel⁴ used FEM for modelling of diffusion and partitioning in biological systems. One-dimensional mathematical models were developed by Lee et al.^5,6 which considers a single component drug penetration via both intercellular and transcellular pathways of the skin. The skin was divided into stratum corneum and viable epidermal layers and diffusion of drug from the porous device into the multiple pathway model of skin was calculated by finite difference approximations. Saxena, Pardasani and their co-workers^7–9 had carried out a notable research for the study of heat distribution in biological tissues. Khanday¹⁰ studied the temperature distribution in multi-layered human skin and sub-cutaneous (SST) by assuming the thermal conductivity as a function of temperature. Khanday et al.^11–14 studied the diffusion of heat and mass in biological tissues particularly dermal regions and human head. Sharma and Saxena¹⁵ used FEM to study the drug distribution in TDD systems in steady state case by making use of quadratic shape function. Most of the models have taken into consideration only the drug transport in the devices and did not discuss the flow of drug in the skin and sub-cutaneous tissues. In this paper, a variational finite element approach with linear shape function is used to obtain the solution of diffusion equation taking into account the role of diffusion constant and the rate constants. The method has been used because of its applicability over the irregular geometries. Since the human skin has irregular geometry, therefore this method guarantees a reasonable outcome as compared to the other numerical methods.

2 Formulation of model

The transport of drugs through transdermal system is governed by the differential equation given in Crank¹⁶(1) $\frac{\partial }{\partial x}\left(D\frac{\partial C}{\partial x}\right)-A(C)-B=\frac{\partial C}{\partial t}$(1) where $C\text{,}D\text{and}t$ denote the drug concentration, mass diffusivity and time. Also, $A(C)$ is the absorption rate of drug by the tissues, which is a function of concentration and B is the rate of drug intake by the blood.

Consider the domain consisting of three layers viz., epidermis, dermis and subcutaneous region with outer and inner boundaries as $a_0$ and $a_3$. The concentration of the drug at the outer skin surface (epidermis) is known as $C_0$ and it has negligible concentration beyond the subcutaneous tissue as the drug is absorbed through the subcutaneous tissue to reach the target site. Hence at the innermost boundary the concentration is $C_3=0$.

By the law of mass action, the absorption rate of drug by the tissue decreases with changing time, so we can say that absorption rate is decreasing function of drug concentration. Hence, we assume(2) $A_i(C^{(i)})=\exp (-k_i{C}^{(i)})\text{,}(i=1\text{,}2\text{,}3)$(2) where $k_i(\ll 1)$ represents the rate constants. The drug distribution has a mild barrier at the interfaces due to non-homogeneous character of tissues in various layers of dermal regions. So,(3) $C^{(i)}={\rho }_i{C}^{(i+1)}{|}_{x=a_i(i=1\text{,}2)}$(3) Here $x\text{and}{\rho }_i(i=1\text{,}2)$ represents the depth below the skin surface and the skin partition coefficients for the drug at the respective interfaces respectively. ${\rho }_i\approx 1$ for most of the drugs¹⁷.

The solution $C(x)$ of the problem can be approximated by polynomial functions of x of higher degrees. Since the distance between the skin layers is of minute length, the higher order shape functions may not contribute a significant change in the drug diffusion, therefore the concentration of drug in each region is approximated by a linear polynomial as(4) $C^{(i)}={\varphi }_i+{\phi }_{i}x\text{,}(i=1\text{,}2\text{,}3)$(4)

3 Mathematical method

The solution of such type of problem has been carried out by various researchers. The numerical solution of the boundary value problem is important for various situations due to the non-availability of the analytical and exact solutions. But, the finite element method is treated to be one of the most reliable methods in numerical systems. The method has advantage over other methods due to position dependent properties of parameters and their flexibility. Moreover, the method is applicable to understand the feasible diffusion in irregular geometries. The variational integral(5) $I=\int F(C\text{,}{C}^{\prime }\text{,}x)\mathit{dx}$(5) in optimum form is equivalent to Euler-Lagrange differential equation given by Myer¹⁸(6) $\frac{\partial F}{\partial C}-\frac{d}{\mathit{dx}}\left(\frac{\partial F}{\partial C^{\prime }}\right)=0\text{;}\text{where}{C}^{\prime }=\frac{\partial C}{\partial x}$(6) On comparing Eq. (1(1) $\frac{\partial }{\partial x}\left(D\frac{\partial C}{\partial x}\right)-A(C)-B=\frac{\partial C}{\partial t}$(1) ) with Euler-Lagrange Eq. (6(6) $\frac{\partial F}{\partial C}-\frac{d}{\mathit{dx}}\left(\frac{\partial F}{\partial C^{\prime }}\right)=0\text{;}\text{where}{C}^{\prime }=\frac{\partial C}{\partial x}$(6) ), we obtain the variational integral(7) $I={\int }_0^{a_3}\left[\frac{D}{2}{\left(\frac{\partial C}{\partial x}\right)}^2+(1+B)C-\frac{k}{2}{C}^2+\frac{1}{2}\frac{\partial C^2}{\partial t}\right]\mathit{dx}$(7) The variational integral on each of the sub domains are given as(8) $I_1=\frac{1}{2}\left[D_1\frac{{(C_1-C_0)}^2}{a_1}+(1+B_1)(C_1+C_0)a_1-\frac{k_1{a}_1}{3}\left(C_0^2+C_1^2+C_0{C}_1\right)+\frac{a_1}{3}\left(2C_0\dot{C_0}+2C_1\dot{C_1}+C_1\dot{C_0}+C_0\dot{C_1}\right)\right]$(8) (9) $I_i=\frac{1}{2}\left\{\frac{D_i}{a_i-a_{i-1}}\left(C_i-\frac{C_{i-1}}{{\rho }_{i-1}}\right)+(1+B_i)\left(\frac{C_{i-1}}{{\rho }_{i-1}}+C_i\right)\left(a_i-a_{i-1}\right)-\frac{k_i}{3}\left(a_i-a_{i-1}\right)\left[\frac{C_{i-1}{C}_i}{{\rho }_{i-1}}+C_i^2+{\left(\frac{C_{i-1}}{{\rho }_{i-1}}\right)}^2\right]+\frac{a_i-a_{i-1}}{3}\left[\frac{\dot{C_{i-1}}{C}_i}{{\rho }_{i-1}}+\frac{C_{i-1}\dot{C_i}}{{\rho }_{i-1}}+2C_i\dot{C_i}+\frac{C_{i-1}\dot{C_{i-1}}}{{\rho }_{i-1}}\right]\right\}(i=2\text{,}3)$(9)

Now, these integrals $I_i$ are assembled to obtain $I={\sum }_{i=1}^3{I}_i$ and I is optimised with respect to the nodal values $C_i(i=1\text{,}2)$ to get the system of equations(10) ${\alpha }_1{C}_1+{\alpha }_2{C}_2+{\beta }_1\dot{C_1}+{\beta }_2\dot{C_2}={\gamma }_1$(10) (11) ${\alpha }_2{C}_1+{\alpha }_3{C}_2+{\beta }_2\dot{C_1}+{\beta }_3\dot{C_2}={\gamma }_2$(11) where the coefficients are given in the Appendix.

The Eqs. (10(10) ${\alpha }_1{C}_1+{\alpha }_2{C}_2+{\beta }_1\dot{C_1}+{\beta }_2\dot{C_2}={\gamma }_1$(10) ) and (11(11) ${\alpha }_2{C}_1+{\alpha }_3{C}_2+{\beta }_2\dot{C_1}+{\beta }_3\dot{C_2}={\gamma }_2$(11) ) in the matrix form can be written as(12) $\begin{array}{cc}\hfill & \mathit{AC}+B\dot{C}=Y\hfill \\ \hfill & \text{where}A=\left(\begin{array}{cc}\hfill {\alpha }_1\hfill & \hfill {\alpha }_2\hfill \\ \hfill {\alpha }_2\hfill & \hfill {\alpha }_3\hfill \end{array}\right)\text{;}B=\left(\begin{array}{cc}\hfill {\beta }_1\hfill & \hfill {\beta }_2\hfill \\ \hfill {\beta }_2\hfill & \hfill {\beta }_3\hfill \end{array}\right)\hfill \\ \hfill & \text{and}Y=\left[\begin{array}{c}\hfill {\gamma }_1\hfill \\ \hfill {\gamma }_2\hfill \end{array}\right]\text{;}C=\left[\begin{array}{c}\hfill C_1\hfill \\ \hfill C_2\hfill \end{array}\right]\hfill \end{array}$(12)

Case:I (Steady State)

Solving equation $\mathit{AC}=Y$, we get the values of the concentration of the drug at the nodes $a_1\text{and}{a}_2$ as $C_1(a_1\text{,}0)\text{and}{C}_2(a_2\text{,}0)$ which in turn on substitution to Eq. (4(4) $C^{(i)}={\varphi }_i+{\phi }_{i}x\text{,}(i=1\text{,}2\text{,}3)$(4) ) are used to determine the values of $C^{(i)}(i=1\text{,}2\text{,}3)$.

Case:II (Unsteady State)

Applying the Laplace transform to the Eq. (12(12) $\begin{array}{cc}\hfill & \mathit{AC}+B\dot{C}=Y\hfill \\ \hfill & \text{where}A=\left(\begin{array}{cc}\hfill {\alpha }_1\hfill & \hfill {\alpha }_2\hfill \\ \hfill {\alpha }_2\hfill & \hfill {\alpha }_3\hfill \end{array}\right)\text{;}B=\left(\begin{array}{cc}\hfill {\beta }_1\hfill & \hfill {\beta }_2\hfill \\ \hfill {\beta }_2\hfill & \hfill {\beta }_3\hfill \end{array}\right)\hfill \\ \hfill & \text{and}Y=\left[\begin{array}{c}\hfill {\gamma }_1\hfill \\ \hfill {\gamma }_2\hfill \end{array}\right]\text{;}C=\left[\begin{array}{c}\hfill C_1\hfill \\ \hfill C_2\hfill \end{array}\right]\hfill \end{array}$(12) ), we have$\begin{array}{cc}\hfill & \mathcal{L}\left(\mathit{AC}+B\dot{C}\right)=\mathcal{L}(Y)\hfill \\ \hfill & \text{or,}A\mathcal{L}(C)+B\left(s\mathcal{L}(C)-C(x\text{,}0)\right)=\mathcal{L}(Y)\hfill \\ \hfill & \text{i.e.,}(A+\mathit{Bs})\mathcal{L}(C)=\mathcal{L}(Y)+\mathit{BC}(x\text{,}0)\hfill \end{array}$(13) $\text{which gives}({\alpha }_1+{\beta }_{1}s){\widehat{C}}_1+({\alpha }_2+{\beta }_{2}s){\widehat{C}}_2=\frac{{\gamma }_1}{s}+{\beta }_1{C}_1(a_1\text{,}0)+{\beta }_2{C}_2(a_2\text{,}0)$(13) (14) $\text{and}({\alpha }_2+{\beta }_{2}s){\widehat{C}}_1+({\alpha }_3+{\beta }_{3}s){\widehat{C}}_2=\frac{{\gamma }_2}{s}+{\beta }_2{C}_1(a_1\text{,}0)+{\beta }_3{C}_2(a_2\text{,}0)$(14) where  ̂ represents Laplace transform.

On solving Eqs. (13(13) $\text{which gives}({\alpha }_1+{\beta }_{1}s){\widehat{C}}_1+({\alpha }_2+{\beta }_{2}s){\widehat{C}}_2=\frac{{\gamma }_1}{s}+{\beta }_1{C}_1(a_1\text{,}0)+{\beta }_2{C}_2(a_2\text{,}0)$(13) ) and (14(14) $\text{and}({\alpha }_2+{\beta }_{2}s){\widehat{C}}_1+({\alpha }_3+{\beta }_{3}s){\widehat{C}}_2=\frac{{\gamma }_2}{s}+{\beta }_2{C}_1(a_1\text{,}0)+{\beta }_3{C}_2(a_2\text{,}0)$(14) ), we have(15) ${\widehat{C}}_1=\frac{({\alpha }_3+s{\beta }_3)(m_1+\frac{{\gamma }_1}{s})-({\alpha }_2+s{\beta }_2)(m_2+\frac{{\gamma }_2}{s})}{(s-\alpha )(s-\beta )}$(15) (16) $\text{and}{\widehat{C}}_2=\frac{({\alpha }_1+s{\beta }_1)(m_2+\frac{{\gamma }_2}{s})-({\alpha }_2+s{\beta }_2)(m_1+\frac{{\gamma }_1}{s})}{(s-\alpha )(s-\beta )}$(16)

Applying the inverse Laplace transform and making use of Heaviside Expansion theorem,^‡ we get the value of $C_i(i=1\text{,}2)$ which on substitution to Eq. (4(4) $C^{(i)}={\varphi }_i+{\phi }_{i}x\text{,}(i=1\text{,}2\text{,}3)$(4) ) is used to determine the value of $C^{(i)}(i=1\text{,}2\text{,}3)$ at any instant of time. The overall drug concentration can be formulated by assembling the linear variations of each region. This indicates polygonal variations of concentration with discontinuities of flux at the interfaces. In order to obtain the better approximations, the Lagrange’s interpolation method has been used. The method determines the drug concentration with continuous flux between the interfaces and within the region.(17) $C(x)=H_0{x}^3-H_1{x}^2+H_{2}x-H_3$(17)

4 Numerical calculations

To solve Eq. (12(12) $\begin{array}{cc}\hfill & \mathit{AC}+B\dot{C}=Y\hfill \\ \hfill & \text{where}A=\left(\begin{array}{cc}\hfill {\alpha }_1\hfill & \hfill {\alpha }_2\hfill \\ \hfill {\alpha }_2\hfill & \hfill {\alpha }_3\hfill \end{array}\right)\text{;}B=\left(\begin{array}{cc}\hfill {\beta }_1\hfill & \hfill {\beta }_2\hfill \\ \hfill {\beta }_2\hfill & \hfill {\beta }_3\hfill \end{array}\right)\hfill \\ \hfill & \text{and}Y=\left[\begin{array}{c}\hfill {\gamma }_1\hfill \\ \hfill {\gamma }_2\hfill \end{array}\right]\text{;}C=\left[\begin{array}{c}\hfill C_1\hfill \\ \hfill C_2\hfill \end{array}\right]\hfill \end{array}$(12) ), we make use of the values given in Sharma and Saxena¹⁵ of the constants given in the Table 1.

Table 1 Numerical values.

	${\rho }_1$	${\rho }_2$	$B_2(\mathrm{\mu }{\text{m}}^{-3})$	$B_3$	$k_1({\text{sec}}^{-1})$	$k_2$	$k_3$
Set-I	0.8	0.9	0.002	0.02	0.4	0.5	0.55
Set-II	0.98	1	0.1	0.2	0.5	0.7	0.8

We can assign different values to the constants $a_i(i=0\text{,}1\text{,}2\text{,}3)$ and $D_i(i=1\text{,}2\text{,}3)$ depending on the sample of the skin under study. The set of values considered here are$a_0=0\mathrm{\mu }\text{m,}{a}_1=0.3\mathrm{\mu }\text{m,}{a}_2=0.6\mathrm{\mu }\text{m,}{a}_3=1\mathrm{\mu }\text{m}$$D_1=0.00068\mathrm{\mu }{\text{m}}^2{\text{s}}^{-1}\text{,}{D}_2=0.0205\mathrm{\mu }{\text{m}}^2{\text{s}}^{-1}\text{,}{D}_3=0.002\mathrm{\mu }{\text{m}}^2{\text{s}}^{-1}$

The numerical calculations have been carried out for 5 grams of drug. The curves have been plotted for $C^{(i)}(i=1\text{,}2\text{,}3)$ versus thickness of skin layers at time $t=0\text{,}5\text{and}10\text{sec}$. The profiles of drug concentration in two different set of values are shown in Figs. 1 and 2.

Figure 1 Concentration profile for Set-I at $t=0\text{,}5\text{and}10\text{sec}$.

Figure 2 Concentration profile for Set-II at $t=0\text{,}5\text{and}10\text{sec}$.

5 Discussion and conclusion

The present work is an attempt to estimate the drug concentration in the human dermal layers. Initially, we have calculated the drug concentration at the nodal points of different dermal layers and graphs have been plotted using MATLAB software between concentration and thickness of dermal layers. From the graphs given in Figs. 1 and 2, the drug concentration profiles at different time intervals have been estimated. It has been observed that the curves for $C^{(i)}(i=1\text{,}2\text{,}3)$ fall more quickly in Set-II and reach the equilibrium state earlier than the Set-I. Analysing the concentration profiles, it seems that the concentration decreases with increase in the partition coefficient, absorption coefficients and rate constants. The curves given in figures reflect the amount of drug absorption in the region of papillary and reticular layers due to dense network of capillary bed.

Lagrange’s interpolation formula has been used and the curve describing the continuous distribution of drug inside the dermal layers has been obtained and is given in Figs. 3 and 4. The results obtained in this paper appear to be realistic and are stronger than the results obtained by Sharma and Saxena¹⁵. They had taken the rate of absorption of drug by the tissues and rate of drug intake by blood as constants while as in our model, keeping the physiological aspects of the dermal layers into consideration, the absorption rate has been taken as an exponentially decreasing function of the concentration. Also they have studied the drug distribution in TDD systems for steady state case while the present work is time dependent problem and therefore their model is a special case of our model. Our work shows a close resemblance and gives modified outcome than the work carried out by Sharma and Saxena¹⁵ for steady state case. Moreover, our model is time dependent and the time plays a key role in drug diffusion.

Figure 3 Modified drug concentration based on Lagrange’s interpolation method for Set-I.

Figure 4 Modified drug concentration based on Lagrange’s interpolation method for Set-II.

The developed model helps us in predicting the amount of the drug concentration in different tissues at different time intervals. Also, this model may be helpful in analysing the role of the diffusion constants and the absorption coefficients in drug delivery system. The study may be very helpful in the medical sciences in various situations.

Conflict of interest

None declared.

Acknowledgements

The authors wish to thank the DST, Govt. of India for providing financial support under DST- INSPIRE programme.

Notes

Peer review under responsibility of Alexandria University Faculty of Medicine.

Available online 25 November 2014

† A physical process in which ion flow diffusively in a medium driven by an applied electric field.

‡ Let $F(s)\text{and}G(s)$ be two polynomials in s where the degree of $F(s)$ is lower than that of $G(s)$ and if $G(s)$ has n distinct roots ${\alpha }_i(i=1\text{,}2\text{,}\dots \text{,}n)$, then${\mathcal{L}}^{-1}\left[\frac{F(s)}{G(s)}\text{;}t\right]={\displaystyle \sum _{i=1}^n}\frac{F({\alpha }_i)}{G^{\prime }({\alpha }_i)}{e}^{{\alpha }_{i}t}$

Appendix A

$\begin{array}{cc}\hfill & h_i=a_{i+1}-a_i\text{,}(i=0\text{,}1\text{,}2.)\text{;}{g}_i=\frac{D_i}{h_i}\text{;}{r}_i=k_i{h}_i\text{,}(i=1\text{,}2\text{,}3)\hfill \\ \hfill & q_j=\frac{g_j}{{\rho }_{j-1}^2}\text{;}{n}_j=\frac{r_j}{{\rho }_{j-1}^2}\text{;}{p}_j=\frac{h_j}{{\rho }_{j-1}^2}\text{,}(j=2\text{,}3)\hfill \\ \hfill & {\alpha }_1=2\left(g_1+q_2\right)-\frac{2}{3}\left(r_1+n_2\right)\text{,}{\alpha }_2=\frac{-1}{{\rho }_1}\left(2g_2+\frac{r_2}{3}\right)\hfill \\ \hfill & {\alpha }_3=2\left(g_2+q_3\right)-\frac{2}{3}\left(r_2+n_3\right)\text{;}{\beta }_1=\frac{2}{3}\left(h_1+p_2\right)\text{;}{\beta }_2=\frac{{\rho }_1{p}_2}{3{\rho }_1}\hfill \\ \hfill & {\beta }_3=2\left({\rho }_1{\beta }_2+\frac{p_3}{3}\right)\text{;}{\gamma }_1=\left(2g_1+\frac{r_1}{3}\right)\hfill \\ \hfill & {\gamma }_1=2g_1+\frac{p_3}{3}\text{;}{\gamma }_2=\left(\frac{n_3}{3}+2q_3\right){\rho }_2{C}_3-(1+B_2)h_2-(1+B_3)p_3{\rho }_2\hfill \\ \hfill & m_i={\beta }_i{C}_1(a_1\text{,}0)+{\beta }_{i+1}{C}_2(a_2\text{,}0)\text{,}(i=1\text{,}2)\hfill \\ \hfill & H_0={\sum }_{i=0}^3\frac{C_i}{{\mathrm{\Delta }}_i}\text{;}{\mathrm{\Delta }}_i={\mathrm{\Pi }}_{{\displaystyle \genfrac{}{}{0ex}{}{i\ne j}{j=0}}}^3(a_i-a_j)\text{,}(i=0\text{,}1\text{,}2\text{,}3)\hfill \\ \hfill & H_1={\sum }_{i=0}^3{H}_{i1}\frac{C_i}{{\mathrm{\Delta }}_i}\text{;}{H}_{i1}={\sum }_{{\displaystyle \genfrac{}{}{0ex}{}{j=0}{j\ne i}}}^3{a}_j\text{,}(i=0\text{,}1\text{,}2\text{,}3)\hfill \\ \hfill & H_2={\sum }_{i=0}^3{H}_{i2}\frac{C_i}{{\mathrm{\Delta }}_i}\text{;}{H}_{02}=a_1{a}_2+a_1{a}_3+a_2{a}_3\text{;}{H}_{12}=a_0{a}_2+a_0{a}_3+a_2{a}_3\hfill \\ \hfill & H_{22}=a_0{a}_1+a_0{a}_3+a_1{a}_3\text{;}{H}_{32}=a_0{a}_1+a_0{a}_2+a_2{a}_3\hfill \\ \hfill & H_3={\sum }_{i=0}^3{H}_{i3}\frac{C_i}{{\mathrm{\Delta }}_i}\text{;}{H}_{i3}={\mathrm{\Pi }}_{{\displaystyle \genfrac{}{}{0ex}{}{j=0}{i\ne j}}}^3{a}_j\text{,}(i=0\text{,}1\text{,}2\text{,}3)\hfill \end{array}$