Continuum Modeling with Functional Lennard–Jones Parameters for DNA‐Graphene Interactions

Carbon nanostructures are of particular interest as platforms for molecular storage and adsorption. In this paper, the adsorption of a single stranded DNA molecule onto a graphene sheet is considered. Even though DNA molecules are complicated heterogeneous structures comprising several types of atoms, it is found that the repeated patterns within the DNA molecules enable the use of a continuum approach to model the DNA‐graphene sheet interaction. Here, a model is proposed such that the heterogeneity across the DNA molecule is captured by interaction functions, which are used to replace the attractive and repulsive constants in the Lennard‐Jones potential. Result from this new model shows better agreement to molecular dynamics simulations compared to the traditional continuum approach where atoms on the DNA are averaged evenly across the molecule. Finally, the paper comments on the model, its parameters, and suggests ways for improvement.


Introduction
Due to their high surface area, carbon materials have been explored as platforms for molecular adsorption for uses in a variety of applications. Gas storage, [1][2][3][4] pollution capture, [5][6][7] biomedical uses, [8][9][10] gas sensing, [11][12][13] and drug encapsulation [14][15][16] are some of the applications that make use of the non-covalent interactions between carbon nanostructures and stored materials. Since carbon compounds are typically nonpolar, the non-covalent interactions between carbon nanostructures and other molecules are dominated by van der Waals forces. To model these interactions and obtain the potential energy for the system, the 6-12 Lennard-Jones potential is commonly used due to its simple form and computational efficiency. [17] DOI: 10.1002/adts.202200896 The 6-12 Lennard-Jones potential for two atoms at a distance apart is given by where A and B are the attractive and repulsive constants, respectively. As shown in Equation (1), the Lennard-Jones potential is an interatomic potential, so for intermolecular interactions, one method to obtain the total potential energy is to sum over the pairwise interaction between atoms i and j of the two interacting molecules, namely This method is referred to as a discrete approach. Another method to obtain the total interaction energy, which we refer to as a continuum approach, is to replace the above double summations with double surface integrals, where here is a distance between typical surface elements dS 1 and dS 2 , 1 and 2 are the surface densities of atoms on the two molecules. This approach was first proposed by Girifalco [18] to determine the potential energy between interacting fullerenes. Over the past decade, extensive development of this approach has been made to determine interaction energies for a variety of molecular structures with various shapes, sizes, and chemical compositions. [19][20][21][22] The key advantage of a continuum model lies in the resultant analytical expression of the potential energy, which can be readily used to determine characteristics of an interaction, such as critical distance, energy minima, and forces at an arbitrary distance and configuration. However, the continuum model is not straightforwardly used for irregular shaped structures.
In the standard continuum approach mentioned above, the coefficients (A and B) in the Lennard-Jones potential are constants, which is suitable for modeling interactions between mononuclear molecules since they can be thought of as homogeneous surfaces with uniform atomic distribution. [23,24] Examples include interactions between two carbon nanostructures, such as fullerene-fullerene, fullerene-nanotube, fullerenenanocone, nanotube-nanotube, and graphene-graphene. [25] For interactions involving heterogeneous structures or heteronuclear molecules, the standard continuum approach can be adapted by either averaging contributions of the constituent atoms as a homogeneous approximation, [26,27] or dividing the molecules into multiple surfaces of the same atomic type. [28,29] Even though these extensions are relatively simple to implement, the homogeneous approximation produces results that poorly agree with molecular dynamics simulations, and the multi-surface method requires multiple integral computations which can be cumbersome for large molecules. Recently, a new development for the continuum approach has been introduced to capture heterogeneous structure where the constants A and B are replaced by functional coefficients known as interaction functions, A( ) and B( ). [30] These functions represent the change in the interaction profile depending on the location of different types of atoms across the heterogeneous molecule. Although excellent agreements to molecular dynamics simulations are obtained, the integration involved becomes more complicated than having A and B as constants. Importantly, care must be taken to determine the interaction functions so that the integration is tractable.
We comment that the method of averaging contribution from all constituent atoms is the simplest to upscale for studying large heterogeneous molecules, whereas the multi-surface method is not as straightforward due to the summation of all interactions from all the surfaces. As such, the multi-surface method can get computationally intensive for large heterogeneous structures. The method with the interaction functions is relatively easy to upscale compared to the multi-surface method as there is only one surface to handle.
In refs. [30][31][32], we use the continuum approach with the interaction functions to model the interactions of carbon structures with coronene and methane, which are small heterogeneous structures comprising only two types of atoms. In this paper, we investigate potential use of interaction functions to model larger heterogeneous structures comprising more than two types of atoms.
In particular, we consider a heterogeneous molecule with clear patterns or repeated structures due to the organization of their atoms, such as a periodic structure in polymers. One example of such structures is a single-stranded DNA (ssDNA) molecule, and so we model its interaction with a graphene sheet. Based on the patterns within a ssDNA molecule, we propose a technique where we combine the method of averaging over regions of homogeneous surfaces and using the interaction functions to describe the interaction profile across the regions within a heterogeneous molecule.
We note that we particularly choose to study the interaction between a ssDNA molecule and a graphene sheet in order to compare our results to that of Alshehri et al. [26] who adopt a continuum approach with homogenous approximation, assuming atoms on the ssDNA molecule are uniformly smeared across DNA's helical surface.
We also note a number of proposed applications involving DNA-graphene interactions, including molecular logic gates, [33] material synthesis, [34] and electronic biosensing. [35] Thus, understanding adsorption mechanisms of DNA onto graphene sheet is crucial for the development of such applications.
In the following section, we describe the newly combined approach taken to obtain analytical expressions for the interaction energy. Note that full details for evaluating integral expressions are given in Appendices A and B, while detail for molecular dynamics simulations used to benchmark analytical solution is given in Appendix C. In Section 3, we investigate the effects of the parameters on the interaction energy, as well as discuss improvements that can be made when implementing this modeling technique to other molecules in the future.

Modeling Approach
In this section, we model the interaction between a singlestranded B-DNA molecule and a graphene sheet (see Figure 1). We assume that a graphene sheet is lying flat on the xy-plane and that the DNA is a helix of length c and radius b. Following Alshehri et al., [26] a general point on the surface of the helix has coordinates ( c
Adv. Theory Simul. 2023, 6, 2200896  While Alshehri et al. [26] average the contribution of atoms on the DNA to obtain constants A and B, this paper considers the atomic structure of the DNA and thus, we use the Lennard-Jones potential with functional interaction coefficients to model the DNA.
Here, the potential energy between a ssDNA molecule and a graphene sheet is evaluated by the expression where g is the atomic density of the graphene sheet, and K n is determined by the computation of the integral where 1 is the atomic density of the ssDNA molecule, is the distance between a point on the surface of the ssDNA and the graphene sheet, dS DNA and dS g are the surface elements of the DNA helix and graphene sheet respectively, and f n for n = 3, 6 are the interaction functions replacing the standard constants A and B. The physical constants of DNA and graphene used in this model are found in Table 1.
Since a graphene sheet is a homogeneous surface of carbon atoms, the choice of interaction function only depends on the structure of the ssDNA molecule. The ssDNA molecule consists of two main regimes along the length of the strand, namely the phosphate-deoxyribose backbone, comprising alternating molecules of phosphate and deoxyribose in a periodic fashion, and the inner region comprising the four (sometimes five) types of nucleobases (see Figure 2a). A single turn of ss-DNA typically has either 10 or 11 nucleobases, often written as 10.5. [36] This allows us to assume the interaction function that is periodic in with a frequency of 11. While there is a variability in the distribution of the atoms throughout the DNA as shown in Figure 2a, the regularity in the distribution of the constituent molecular groups is clear, as shown in Figure 2b. Thus, it seems intuitive to model each molecular group on the DNA as a homogeneous region. So, for example, rather than attempting to capture the heterogeneity within a phosphate group, we treat the entire phosphate as a homogeneous surface with average contribution of phosphorous and oxygen. Additionally, since we do not consider a particular DNA, specific sequence of nucleobases is not relevant in the model. As such, we assume that all four nucleobases have an equal likelihood to be on the DNA and thus we approximate an averaged nucleobase as a homogeneous surface over which its constituent atoms are smeared and the contributions of all four bases are averaged.
Based on these assumptions, there are three possible candidates for the interaction function: a) a piecewise interaction function that varies in (i.e., f n,1 ( ) for t ∈ [0, ] and f n,2 ( ) for t ∈ [ , 1]); b) a piecewise function that varies in t (i.e., f n,1 (t) for ∈ [ (2i)2 11 , (2i+1) 2 11 ] and f n,2 (t) for ∈ [ (2i+1) 2 11 , (2i+2) 2 11 ] for i = 0, 1, ⋯, 5, not including the final interval as there is an odd number of strips); and finally c) an interaction function that varies continuously in both variables f n ( , t). Visually, (a) and (b) correspond to the images in Figure 3, where the first is two continuous surfaces separated along t = , the second image has 11 continuous surfaces stacked vertically along the helix, and (c) would be a fully continuous combination of (a) and (b). This division into regions also naturally leads to using individual densities ( i ) for each region rather than a single density for the whole DNA strand.
According to (a), we consider interaction functions of the form cos(m ) + in order to capture the alternating nature of the two regimes, phosphate and deoxyribose in the outer band and nucleobases and empty space in the inner band. The resulting integral under this choice of the interaction function is given by where the values of the coefficients i and i (i = 1, 2) are determined from the interaction coefficients (A and B) of the con-stituent molecules interacting with graphene, and m i is determined by the number of base pairs in one turn. According to (b), we consider interaction functions that are sigmoidal in shape since we need to capture the transition from nucleobase to deoxyribose and empty space to phosphate across the strips. Interaction functions that are sigmoidal in shape have previously been used to capture transition between two types of atoms [30] so we assume similar behavior to hold when transitioning between two regions within a molecule. In particular we use functions of the form arctan( (t − )) + , which implies that the resulting integral for evaluation is of the form where the coefficients i and i are determined in a similar manner to the first choice above, 1,1 and 1,2 are atomic densities for each of the horizontal strips, and i is determined by the slope of the transition from the inner to outer region. According to (c), we consider an interaction function that has combined properties from both (a) and (b), that is, we expect a sigmoidal profile when fixing and a sinusoidal shape when fixing t. This behavior can be achieved by taking the cosine function above and replacing the coefficients and with arctan functions, namely ( 1 arctan( 1 (t − )) + 1 ) cos(m ) + ( 2 arctan( 2 (t − )) + 2 ), where the coefficients of the arctan functions are determined from the coefficients of the original cosine functions from the first choice. Note that this only works if both cosine functions have matching frequencies m, otherwise the frequency will need to be modified by some function www.advancedsciencenews.com www.advtheorysimul.com of t that smoothly transitions from m 1 to m 2 . With this choice of interaction functions, the integral K n of the form Of all these three choices of interaction functions and the corresponding integrals, we find that only Equation (7) can be evaluated analytically, namely where and where 2 F 1 denotes the standard hypergeometric function and F 1 denotes an Appell's hypergeometric function. Detailed derivation of these expressions is given in Appendix A. For Equation (8), we obtain a semi-analytical solution where one integration is remained to be evaluated numerically (see Appendix B). Due to the high computational time to compute the semi-analytical solution numerically, we do not adopt this case here. We also do not progress with Equation (9) since we are only interested in analytical solution for K n . High computational time is also expected for Equation (9).

Determination of Coefficients 1,2 and 1,2 Used in Equation (10)
In order to adopt (10) to determine the potential energy between a DNA and a graphene sheet, we first need to determine the coefficients i and i (i = 1, 2), which are the interaction coefficients of the constituent molecular groups interacting with a graphene sheet. We note that two sets of the interaction coefficients are considered. The first set is derived from the carbon-carbon values in Table 2, while the second set uses the graphene-graphene values.
In refs. [30,31] methane-nanotube and coronene-graphene interactions were considered. In these studies, the heterogeneous molecule (methane/coronene) can be modeled as two regions (inner region of carbon atoms and outer region of hydrogen atoms). The interaction functions are then used to describe the interaction profile of a carbon surface (nanotube/graphene) interacting with the two regions. In refs. [30][31][32] the sigmoidal  Molecule functions were used to represent a smooth transition from the inner carbon core to the outer hydrogen region. As shown in refs. [30][31][32] the coefficients for the interaction functions are relatively simple to determine since they can be restricted to match the Lennard-Jones constants A C-C (or B C-C ) at the carbon region and A C-H (or B C-H ) at the hydrogen region. For our choice of interaction functions presented in Equation (7) for modeling a ssDNA molecule, we determine i and i (i = 1, 2) such that when n = 3, 1 + 1 = A nucleobase-x and 2 + 2 = A deoxyribose-x and when n = 6, 1 + 1 = B nucleobase-x and 2 + 2 = B deoxyribose-x . Note that x represents interaction involving either carbon or graphene. Similarly, we prescribe that 1 − 1 = 0 (when n = 3) or 0 (when n = 6) and 2 − 2 = A phosphate-x (when n = 3) or B phosphate-x (when n = 6).
We derive the Lennard-Jones constants of each region by employing the homogeneous smearing method using the parameters given in Table 2. Using the standard combination rules for Lennard-Jones parameters, we generate two sets of interatomic interaction constants A i-C , B i-C , and A i-graphene , B i-graphene , where i is any other type of molecular groups. For example, the attractive constant for carbon and phosphate (PO 4 ) is calculated from Using the values from Table 2, we derive the interaction coefficients for the constituent molecules of DNA interacting with either carbon or graphene as given in Table 3.
Using the values shown in Table 3, we obtain 1,2 and 1,2 for both sets of coefficients as presented in Table 4. In the following section, we plot the energy obtained from Equation (10) using these two set of coefficients. The results are compared with molecular dynamics simulations and the homogeneous model where atoms in the DNA are smeared over the helical surface.

Results and Discussion
In this section, the interaction between a ssDNA and a flat graphene sheet is considered. Particular attention is made to the Table 4. Numerical values for 1,2 and 1,2 for the interaction function (7) derived from the values in Table 3. Note the set of coefficients (i) and (ii) are based on the interaction involving carbon-carbon and graphenegraphene, respectively.  DNA structure in order to use a continuum approach to study such interaction. In Alshehri et al., [26] the DNA is modeled as a homogeneous helical surface where the effect of constituent atoms is averaged throughout the entire structure. As shown in Figure 4, the potential energy obtained from the homogeneous model does not agree well with the results from molecular dynamics (MD) simulations (see Appendix C for details of the MD study). We compare our new results from the expression derived from Equation (10) using the two sets of coefficients found in Table 4, denoted as models (i) and (ii) respectively, with the results from the homogeneous model described in Alshehri et al. [26] and the numerical results from the MD study. For the interaction coefficients used in the homogeneous model, Alshehri et al. [26] considered a specific DNA sequence, 5 ′ -CCACTAGTGG-3 ′ , so the attractive and repulsive constants are obtained by taking the arithmetic mean of all the individual  atomic interaction coefficients (e.g., A DNA−C = 1 320 (97A C−C + 10A P−C + 114A H−C + 61A O−C + 38A N−C )). We note that the total number of atoms in a DNA is determined from the number of atoms in each molecular group in its bonded state, meaning that there are less hydrogen atoms than the standalone molecules. In general, for every base pair there is a corresponding deoxyribose molecule (C 5 OH 7 ) and a phosphate molecule (PO 4  O as the chemical composition of an averaged nucleobase when determining the interaction coefficients. It can be seen that model (i) performs even more poorly compared to the homogeneous model. For model (ii), the result is closer to the MD simulations. This indicates that using coefficients based on graphene-graphene interaction is more appropriate than using coefficients based on C-C interaction for modeling interactions involving graphene sheet.
Next, we explore the impact of other parameters in the model on the interaction energy profile. We first investigate the two shape parameters of the helix, namely the radius b and the length c. From Figure 4, we find the location of the minimum energy ( min ) and the corresponding minimum energy (E min ) as shown in Table 5. The values of min from the three models are found to be in good agreement with the MD results, despite varying widely in the values of E min .
We show in Figure 5 that the location min is strongly affected by the radius of the ssDNA (b). In Figure 5, we demonstrate that varying the values of b by 0.2 Å results in the same amount of displacement in min , as well as a slight change in E min due to helical surface being further away from the graphene sheet.
On the other hand, it is shown in Figure 6 that an increase in the length of the helix, c, over 10 Å is required to produce a decrease in E min of a similar magnitude to that produced by increasing the radius, b by 0.2 Å.
These results demonstrate the significance of DNA's radius over its length in the model. However, these findings indicate that the size of the helix largely impacts the location of the minimum energy ( min ), but only minimally affects the value of minimum energy (E min ).
To further investigate the cause of the discrepancies in the energy profiles shown in Figure 4, we consider interaction coefficients, A and B, used for interactions involving phosphate, deoxyribose, and nucleobase. We comment that the coefficients A and B calculated from averaging the contribution of constituent atoms within a molecule do not provide a result that is in good agreement with the MD simulations. Due to the non-linear effect that the distance has on the Lennard-Jones potential, the configuration of the atoms within the molecule needs to be accounted for. In Figure 7, two homogeneous sphere models for phosphate interacting with graphene are compared against simulation results. One has coefficients determined from averaging the contribution of each constituent atoms (i.e., A = 1 5 A PC + 4 5 A OC ) and the other has coefficients determined by fitting a standard sphere-plane continuum model [25] with the externally obtained well depth ( ) and the van der Waals radius ( ) (i.e., E min = − and E ′ (2 1 6 ) = 0). From the energy profiles, it can be seen that the average atomic contribution approach leads to interaction coefficients A phosphate-carbon and B phosphate-carbon that do not produce result that is in good agreement to simulation results. This is contrasted with the coefficients obtained from fitting a sphere-plane model, which matches the simulation results more closely. Even though there is a big difference between the two approaches, the fitted coefficients approach comes with some drawbacks. These include the needs for a means of accounting for the molecule's shape for a homogeneous model and for information about the well depth and the location of the well depth for the interaction considered. Thus, this method would only be suitable for determining coefficients of smaller molecular units as part of a larger molecule. As a result, we do not pursue the fitted coefficients approach for DNA molecule studied here as it requires further assumptions for the shapes of nucleobases and deoxyribose and the determination of and for these molecular groups and graphene. Another possibility for the discrepancy between the continuum model and the simulation results is the assumption made to the structure of the DNA molecule. Figure 8a shows how the helical-surface model, with the right parameters, captures the overall shape of the ssDNA molecule but not the volume. As the nucleobases are orthogonal to the backbone and the backbone itself is not flat, a model comprising two helical cylinders capturing nucleobases and backbone atoms may be more suitable for modeling a ssDNA strand. An example of this modeling approach is shown in Figure 8b.

Concluding Remarks
The Lennard-Jones potential is widely used for determining interaction energy between non-bonded molecules. This potential together with a continuum approach are adopted when dealing with regular shaped molecules since the resultant energy can be found explicitly as a function of distance between two interacting molecules. Challenges arise when the interaction involves irregular shaped structures or heterogeneous molecules, which comprise more than one types of atoms. For heterogeneous molecules, the assumptions underlying traditional continuum approach that atoms are evenly smeared throughout the molecule and the effects of atoms are the same regardless of its location on the molecule seem invalid. To address this issue, we propose a new continuum approach for a heterogeneous structure, where the Lennard-Jones coefficients are replaced by interaction functions. [30][31][32] While this approach has been shown to work well with methane and coronene, we note their relatively simple structures, which can be modeled as a sphere and a disk, respectively. Furthermore, their atomic arrangement of carbon atoms at the central region and hydrogen atoms at the perimeter enables the interaction functions in the form of a sigmoidal profile, which makes the integrals in the interaction energy traceable yielding analytical outcomes. For larger and more complex heterogeneous molecules, there has yet to be a continuum method that takes into account both atomic arrangement and structure of a molecule. As such, this paper can be viewed as an introduction to a continuum approach using the Lennard-Jones potential with interaction functions to model complicated heterogeneous molecules, such as polymers and DNAs.
A continuum model for heterogeneous molecules is proposed here for the modeling of single-stranded DNA interacting with a graphene sheet. By recognizing the patterns of atoms on the DNA molecule, we are able to develop a technique combining the method of averaging atomic contributions over regions of surfaces and using the interaction functions to describe the interaction profile across the regions within a heterogeneous molecule.
This new approach appears to improve the accuracy of the model when benchmarking with molecular dynamics simulations. Despite some improvement shown here, the model can be further benefited by a systematic method for determining the interaction coefficients between the constituent molecules (namely phosphate, deoxyribose, and the nucleobases) and the graphene sheet. Different interaction functions may also be explored to better capture the heterogeneity of the DNA molecule. Additionally, the model can also be further improved by using a volume-based model to better capture all atoms of the DNA.
These suggested improvements to the continuum modeling for DNAs provide areas for future research directions. Finally, we note that the approach presented in this paper can be extended to model polymers or other heterogeneous molecules with repeated patterns and structures.

Appendix A: Computation of the Integral K n in Equation (7)
A ssDNA helix has coordinates ( c 2 , bt sin , bt cos + ), where c is the length of the DNA strand, b is the radius of the DNA strand, is the distance from the helix axis to graphene sheet. Here, we assume that > b since the DNA molecule cannot intersect the graphene sheet. For a graphene sheet lying flat on the xy-plane, its coordinates are given by (x, y, 0). Typical surface elements for graphene and DNA are given respec- From the coordinates of graphene and the DNA, distance between two typical surface element is given by 2 (bt cos + ) 2 . Thus, the integral K n over both surfaces can be written as Next, we write 2 = fl 2 1 + ( c 2 − x) 2 , and make the substitution x − c 2 = fl 1 tan . This changes the limits of the integration from (−∞, ∞) to [− 2 , 2 ] and the differential dx = fl 1 sec 2 d . With the substitution the in-tegral K n becomes Note that ∫ 2 0 cos 2n xdx = (2n−1)!! (2m)!! 2 and n!! is the double factorial. Now, we again write fl 2 1 = fl 2 2 + (y − bt sin ) 2 and substitute y − bt sin = fl 2 tan , so since ∫ 2 0 cos 2n+1 xdx = (2n)!! (2m+1)!! www.advancedsciencenews.com

www.advtheorysimul.com
The rest of the analytical calculation depends on the choice of f n ( , t). Dividing the DNA strand into two regimes along the t-axis, we have the nucleobases in the interval [0, ] and the phosphate-deoxyribose backbone in the interval [ , 1]. Accordingly, f n is set to be two different functions of in these intervals, namely f n,1 = 1 cos(m 1 ) + 1 in [0, ] and f n,2 = 2 cos(m 2 ) + 2 in [ , 1]. This also splits K n into two double integrals, 1 cos (m 1 ) + 1 (bt cos + ) 2n−2 d dt The integral in is effectively the same for both regions, so we only need to solve the integral of the form After applying a double-angle trigonometric identity, the integral I 2 can be readily computed as a hypergeometric function using a special case shown in equation 3.681.1 in ref. [39], which we then expand into an infinite sum in order to facilitate the integration with respect to t: Using a quadratic transformation of the 2 F 1 hypergeometric function, both the coefficient and the argument of the hypergeometric function in I 2 are simplified to The integral I 1 is more interesting due to the cosine term. By making the assumption that m is an integer renders it tractable as we can make use of De Moivre's formula: where i is the imaginary unit. Thus, we have Combining the results for I 1 and I 2 , the integral I can be written as a function of t with parameters , , and m, namely Thus, the expression for K n is given by which means only four integrals remain, namely In some of the above integrals, namely J 1 and J 3 , we approximate √ 1 + 4b 2 2 c 2 t 2 by its Taylor series about a point t 0 to render them more tractable. The general series is given by so, for = 4b 2 2 c 2 and r = 1 2 , we have As a general rule, we select t 0 to be the midpoint of the interval over which we are integrating. This means t 0 = 2 for J 1 and t 0 = 1− 2 for J 3 . We first evaluate J 2 since it is the most straightforward integral compared to the others. From equation 3.254.1 in ref. [39] we have ∫ u 0 x −1 (u − x) −1 (x 2 + 2 ) dx = 2 u + −1 B( , ) 3 F 2 (− , 2 , +1 2 ; + 2 , + +1 2 ; −u 2 2 ), where p F q is a generalized hypergeometric function given by p F q (a 1 , ⋯, a p ; b 1 , ⋯, b q ; z) = ∑ ∞ i=0 (a 1 ) i ⋯(a p ) i Thus, where the reduction in order of the 3 F 2 to a 2 F 1 occurs as the . As such, Next, we consider J 1 . From equation 3.194.1 in ref. [39] we have ∫ u 0 x a−1 (1+bx) dx = u a a 2 F 1 ( , a; 1 + a; −bu), so that www.advancedsciencenews.com www.advtheorysimul.com Last, we evaluate J 3 by using a linear substitution, t = (1 − )u + , and rewrite the integral in the form of an Appell's F 1 function, that is, = B( , )F 1 ( , , , +  ; u, v). Appell's F 1 function is a generalized 2 F 1 function with two variables instead of one. Thus, Using the results of J 1 , J 2 , J 3 , and J 4 , we obtain analytical expression for K n as presented in Equation (10).