Forecasting sensitive targets of the kynurenine pathway in pancreatic adenocarcinoma using mathematical modeling.

In this study, a new mathematical model was established and validated to forecast and define sensitive targets in the kynurenine pathway (Kynp) in pancreatic adenocarcinoma (PDAC). Using the Panc-1 cell line, genetic profiles of Kynp molecules were tested. qPCR data were implemented in the algorithm programming (fmincon and lsqnonlin function) to estimate 35 parameters of Kynp variables by Matlab 2017b. All tested parameters were defined as non-negative and bounded. Then, based on experimental data, the function of the fmincon equation was employed to estimate the approximate range of each parameter. These calculations were confirmed by qPCR and Western blot. The correlation coefficient (R) between model simulation and experimental data (72 hours, in intervals of 6 hours) of every variable was >0.988. The analysis of reliability and predictive accuracy depending on qPCR and Western blot data showed high predictive accuracy of the model; R was >0.988. Using the model calculations, kynurenine (x3, a6), GPR35 (x4, a8), NF-kβp105 (x7, a16), and NF-kβp65 (x8, a18) were recognized as sensitive targets in the Kynp. These predicted targets were confirmed by testing gene and protein expression responses. Therefore, this study provides new interdisciplinary evidence for Kynp-sensitive targets in the treatment of PDAC.

INTRODUCTION

Pancreatic adenocarcinoma (PDAC) is the 12th most common cancer worldwide. According to The Surveillance, Epidemiology, and End Result program database reports, the 5‐year survival of PDAC patients is around 7.7% with very low improvement. Over 82% of diagnosed patients are at risk of death within a few months because warning signs of PDAC are often not noticeable until the disease has already metastasized to other organs of the body. Statistics alarm that PDAC may become the second leading cause of death among all cancers by 2020. The increase in mortality is due to the lack of effective therapeutic options. The kynurenine pathway (Kynp) plays a significant role in enhancing PDAC growth, immune escape, and immunotolerance. , Plasma levels of kynurenine, kynurenic acid, and quinolinic acid were found significantly higher in PDAC patients, and the level of tryptophan (the substrate of kynurenine) was very low, suggesting high activity of the Kynp in PDAC. Further, Kynp inhibition showed a promising clinical significance among pancreatic cancer patients (NCT02077881, NCT03414229, NCT03432676). The Kynp activates a variety of signaling molecules in cancer cells. It mainly activates G protein–coupled receptor 35 (GPR35) in the PDAC cell that essentially enhances oncogenic signaling and tumorigenesis progression. GPR35 has the potential to stimulate β‐catenin activity in colon, breast, and pancreatic cancers. , , In the PDAC cell, β‐catenin induces IκB kinase β (IKK‐β) and nuclear factor kappa B (NF‐κβ) expression. , NF‐κβ p65 and interleukin 6 (IL‐6) are consecutive molecules of NF‐κβ in the case of Kynp activity. Transforming growth factor beta 1 (TGF‐β1) has been reported as an essential biomarker of the Kynp in PDAC. , , , Therefore, we have delineated the Kynp in PDAC according to the previous literature mentioned above. However, the most sensitive molecules of the Kynp in PDAC have not been defined yet. The exploration of Kynp‐sensitive molecules will substantially contribute to exploring new targets in the treatment of PDAC.

Because the intricacies of sensitive target determination in the whole pathway are challenging, mathematical modelling was established to make the biological changes of Kynp molecules become more quantitative. Different kinds of mathematical models have been established to understand the corresponding mechanism for some diseases. , Hence, the mathematical model has the potential to provide a reasonable explanation for the activity of Kynp molecules in PDAC. In this study, we established a mathematical model to determine sensitive and insensitive molecules of the Kynp in PDAC accurately, which can be considered as critical targets for treating PDAC.

MATERIALS AND METHODS

Cell culture and Kynp stimulation

Human PDAC cell line (Panc‐1 ATCC® CRL‐1469™) was purchased from the ATCC. Panc‐1 cells were plated in Dulbecco's Modified Eagle Medium (DMEM) basic (1×) (Gibco). DMEM was supplemented by 10% fetal bovine serum (FBS) (Gibco) and 1% penicillin/streptomycin (ATLANTA Biologicals). Next, cells were collected and washed twice with phosphate‐buffered saline (PBS). Then, 1 × 106 cells/well were cultured in a 48‐well plastic plate (Corning Costar). Thereafter, cell culture medium was supplemented by 1 µg/mL l‐Tryptophan (Sigma‐Aldrich) and incubated for 72 hours. A total of 20 ng/mL of human purified recombinant indoleamine 2,3 dioxygenase1 (hIDO1) (Prospec) was added every 6 hours to stimulate the Kynp in the tumor cells. This culture system was duplicated and repeated three times.

Gene expression profiling

Culture medium was aspirated carefully into numbered sterile plastic tubes and kept in the refrigerator at 4°C for further analysis. Cells were detached by trypsin and washed with sterile PBS twice. Then, total RNA of tumor cells was extracted by TRIzol standard protocol. Next, the RNA was reverse transcribed to cDNA according to kit instructions (Takara). Later, the template cDNA was used with the primers of the following genes: kynurenine 3‐monooxygenase, GPR35, β‐catenin, IKK‐β, NF‐κβp105, ANK, RelA, IL‐6, and TGF‐β1. Forward and reverse primers listed in Table S1 were designed using NCBI gene IDs and Primer Bank website (https://pga.mgh.harvard.edu/primerbank/). Relative mRNA expression was calculated according to 2−∆∆ct.

Protein expression levels

Enzyme‐Linked Immunosorbent Assay (ELISA)

To validate our model at protein expression level, four molecules of the Kynp have been tested at twelve time points using ELISA kits from 6 hours to 72 hours. Supernatant was aspirated into sterile tubes and then saved in the refrigerator at 4°C. Later, the protein expression levels of kynurenine (Biocompare), β‐catenin (BioVision), NF‐κβp65 (Path‐ScanR), and TGF‐β1 (Multi‐Sciences) in the culture supernatants were tested according to the kit's manufacturers’ instructions.

Western blot

In terms of model verification, tumor cells were cultured by the same protocol as described above. Briefly, eight time points were selected for testing protein expression (6 hours, 10 hours, 20 hours, 30 hours, 40 hours, 50 hours, 60 hours, and 80 hours), and then cells were collected and washed twice with PBS. Total protein was extracted by radioimmunoprecipitation assay (RIPA) buffer mixed with 1% phenyl methyl sulfonyl fluoride (PMSF) proteinase inhibitor for 30 minutes on ice (vortex interval: 10 minutes) and then centrifuged for 10 minutes at 106 g/min. The supernatant was transferred into 1‐mL Eppendorf tubes. Total protein samples were loaded by protein‐loading buffer at 100°C for 10 minutes. Next, 50 ng/mL of total protein was run on the gel and then transblotted on the PVDF membrane. Thereafter, goat milk fixation was performed. Then, protein bands were visualized by human anti‐GPR35 antibody (NOVUSBIO), human anti‐β‐catenin antibody (Proteintech), human anti‐NF‐κβp105 antibody (Cell Signaling Technology), and human anti‐TGF‐β1 antibody (Abcam). Then, they were incubated with anti‐rabbit IgG for 2 hours, washed three time with Tris‐buffered saline buffer, and then were visualized by Odyssey FC (LICOR).

Mathematical model establishment

A mathematical model for the Kynp in PDAC cells was proposed based on previous literature that mention the molecules which actively participate in the Kynp in different tumor types. These molecules were all tested in the Panc‐1 cell line. According to the gene expression of these molecules, a mathematical model was proposed as the following equation system considering concentration and time of expression for every variable: where ,,,,,,,,,,, and represent the relative expression of Tryptophan, IDO1, Kyn, GPR35, β‐catenin, IKK‐β, NF‐κβp105, NF‐κβp65, ANK, RelA, IL‐6, and TGF‐β1, respectively. Parameters (,,,,,,,,,,,,,,,,,,,,,,,) represent production and degradation rate of Tryptophan, IDO1, Kyn, GPR35, β‐catenin, IKK‐β, NF‐κβp105, NF‐κβp65, ANK, RelA, IL‐6, and TGF‐β1, respectively. denotes the activation rate from Tryptophan and IDO1 to Kyn. denotes the activation rate from Kyn to GPR35. denotes the activation rate from GPR35 to β‐catenin. denotes the activation rate from β‐catenin to IKK‐β. denotes the activation rate from IKK‐β to NF‐κβp105. denotes the activation rate from IKK‐β to NF‐κβp65. denotes the activation rate from NF‐κβp105 to ANK. denotes the activation rate from NF‐κβp65 to RelA. and denote the activation rate from ANK and RelA to IL‐6, respectively. denotes the activation rate from IL‐6 to TGF‐β1.

The estimation of model parameters

In order to get the value of the parameters in Equation (1) based on relative mRNA expression levels, an algorithm (fmincon and lsqnonlin function) was programmed to estimate the 35 parameters by Matlab 2017b. In this study, all the parameters were defined to be non‐negative and bounded because each parameter has its own biological significance. Secondly, based on experimental data, the fmincon function was employed to estimate the approximate range of each parameter. The estimated parameters by the fmincon function were regarded as the initial values. Thereafter, further estimation was performed using the lsqnonlin function to achieve the best fitting effect between the simulation curve and the experimental data curve.

Sensitive analysis

In order to determine the sensitive parameters in the Kynp, the extended Fourier amplitude sensitivity test (EFAST) was employed. The total variance was acquired as following: where is the variability associated with the main effect of parameter , and is the variability associated with the interaction between parameters and . The first‐order sensitivity indexes () are derived from the decomposition of the above equation by dividing the importance measures by : where is the interactive sensitivity index between parameters and . Thus, based on the effects of its interaction with other parameters, the total sensitivity indexes () were given by:

A higher or indicates a higher contribution to the response variables of interest, and vice versa. and of thirty‐five parameters in the Kynp were calculated by Matlab, and the codes are presented in Appendix S1. Based on Cannavo, the sensitivity value () larger than 0.5 was defined as a sensitive parameter; otherwise, it was defined as an insensitive parameter.

Biological analysis of sensitive and insensitive factors

In order to confirm the accuracy of sensitivity analysis calculations, experiments were performed using siRNA technique to knock down the expression of sensitive and insensitive parameters, and then TGF‐β1 expression was assessed as an indicator of Kynp activity. Based on the sensitivity analysis results, three siRNA sequences were designed for sensitive and insensitive genes as listed in Table S2. The experiment was divided into three groups: siRNA group, negative control (NC) group, and normal group. For sensitive parameter evaluation, a real‐time polymerase chain reaction (RT‐PCR) experiment was done to test the relative mRNA expression of TGF‐β1 under the condition of knocking down Kyn, GPR35, NF‐kβp105, and NF‐kβp65 gene expression using siRNA. For insensitive parameter evaluation, similarly, the relative expression of TGF‐β1 was tested under the condition of knocking down β‐catenin, IKK‐β, and IL‐6 gene expression using siRNA.

Panc‐1 cell line was cultured in a 48‐well plate in DMEM without antibiotic for 24 hours. Next, cells were transfected by siRNA using the protocol of lipofectamine‐2000 (Thermo Fisher). Later, cells were incubated for 18 hours and then washed, collected, and used for extracting total RNA. The primers listed in Table S1 were used to assess target gene expression. Effective siRNAs have been selected for sensitive analysis test. Kynureninase has been targeted by Kyn‐siRNA, GPR35R by GPR35‐siRNA, NF‐κβp105 by NF‐κβp105‐siRNA, and NF‐κβp65 by NF‐κβp65‐siRNA. Thereafter, transfected cells were collected, washed, and used to test TGF‐β1 expression by RT‐PCR and Western blot using the same protocol as described above. Moreover, insensitive parameters (β‐catenin, IKK‐β, and IL‐6) were assessed by the same method as the one used for sensitive parameters using β‐catenin‐siRNA, IKK‐β‐siRNA, and IL‐6‐siRNA, respectively.

Statistical analysis

All data were presented as the mean ± SD. Differences between the values were evaluated using one‐way analysis of variance (Tukey's post hoc test), except in vivo biocompatibility experiment (Tukey's post hoc test). Data were presented as means with 95% confidence interval. P < .05 was considered statistically significant.

RESULTS

Mathematical model of the Kynp was successfully established based on mRNA expression trends of Kynp molecules

Relative mRNA expression of nine consecutive molecules in the Kynp was assessed for 72 hours, in 6‐hour intervals, as seen in Figure 1A‐I and Table S1. The expression of the tested molecules showed obvious changes under the effect of IDO1 in a time‐dependent manner. It showed significant increase from 6 hours to 72 hours that enhanced the increase of TGF‐β1 in the Panc‐1 cell line. Based on the relative mRNA expression, thirty‐five parameters were estimated by the fmincon function and lsqnonlin function (Matlab 2017b) to establish the model according to expression trends (see the list in Data S1). The model was validated by plotting the experimental findings (red line) and model simulations (blue line) as shown in Figure 1J‐U. We noticed that experimental findings and model stimulation were well fitted. The correlation coefficient (R) was computed between experimental findings and model stimulation. Results as shown in Table S4 illustrated that the R value of all tested variables displayed strong correlations (R > 0.988) suggesting that the model is efficient.

FIGURE 1

Using qPCR the gene relative expression of Kynp parameters in Panc‐1 was analyzed from 6 h to 72 h. A, Expression level of kynurenine (x3). B, Relative expression level of GPR35 (x4). C, Gene relative expression of β‐catenin (x5). D, Gene relative expression of IKKβ (x6). E, Gene relative expression of NF‐κβp105 (x7). F, Gene relative expression of ANK (x9). G, Gene relative expression of RelA (x10). H, Gene relative expression of IL‐6 (x11). I, Gene relative expression of TGF‐β (x12). Then, the mathematical model was established by parameter estimation at the mRNA expression level. J, K, Q, Simulation trends of tryptophan, IDO1, and NF‐κβp65 are presented as blue line from 0 h to 72 h based on the established model. L‐P, R‐U Validation of Kyn, GPR35, β‐catenin, IKK‐β, NF‐κβp105, ANK, RelA, IL‐6, and TGF‐β1 was performed by plotting the experimental findings, shown as red line, and model simulation, shown as blue line 0 h to 72 h (error bars, mean ± SD)

Moreover, in order to test the accuracy and predictability of the model, further experiments were performed at time points 10 hours, 20 hours, 40 hours, and 80 hours, as shown in Figure 2A‐I and Figure S1A‐I. The expression of Kynp molecules confirmed the significant upward changes in a time‐dependent manner. The expression trends of Kynp molecules at time points 10 hours, 20 hours, and 40 hours were used to test the accuracy of the model, while the experimental data at time point 80 hours was used to test predictability of the model. The results, as seen in Figure 2J‐U, showed high accuracy of the model simulation from 0 hours to 72 hours and efficient predictability in the trends of Kynp molecule expression from 72 hours to 80 hours. Further, the calculation of R values between mathematical model simulation and experimental results showed R larger than 0.988, as seen in Table S5, which indicates the high accuracy of the established model in the prediction of the activity of Kynp molecules in PDAC.

FIGURE 2

qPCR analysis of the accuracy and predictability of the Kynp math model was performed at four time points from 10 h to 40 h for an accuracy test and from 72 h to 80 h for the predictability assessment. A‐I, Gene relative expression of the Kynp variables at 10 h, 20 h, 40 h, and 80 h. Then, the verification of the established mathematical model was performed based on the mRNA expression level. J, K, Q, Simulated trends of Tryptophan, IDO1, and NF‐κβp65 are presented as blue line from 0 h to 80 h based on the established model. L‐P, R‐U, Accuracy test of Kyn, GPR35, β‐catenin, IKK‐β, NF‐κβp105, ANK, RelA, IL‐6, and TGF‐β1 was presented as red points at time points 10 h, 20 h, 40 h. Consequently, the predictability was presented by the red points at time points 80 h. Model simulations were correspondingly shown as blue and green line, respectively

Protein expression trends of Kynp molecules confirmed the established model

In order to verify the mathematical model at protein expression level, ELISA and Western blot were performed. The protein expression levels of kynurenine, β‐catenin, NF‐κβp65, and TGF‐β1 were tested for 72 hours in 6‐hour intervals by ELISA kits. The results showed increased protein expression trends of kynurenine, β‐catenin, NF‐κβp56, and TGF‐β1 from 6 hours to 72 hours, as presented in Figure 3A‐D, that interestingly were in line with the genetic profiling results. Western blot analysis, as shown in Figure 3E and Figure S3A‐D, presented the same increasing trend of protein expression of GPR35, β‐catenin, NF‐κβp105, and TGF‐β1 from 6 hours to 80 hours.

FIGURE 3

Mathematical model verification based on the protein expression level. A‐D, ELISA results showed increased protein expression levels of kynurenine, β‐catenin, NF‐kβp65, and TGF‐β1 from 6 h to 72 h. E, Western blot represents protein expression levels of GPR35, β‐catenin, NF‐kβp105, and TGF‐β1 from 6 h to 80 h. Next, the sensitivity analysis was performed using numerical computation. F, The first order sensitive indexes (S) of 35 parameters. G, The total sensitive indexes (S) of 35 parameters calculated by extended Fourier amplitude sensitivity test (EFAST)

The prediction of sensitive molecules in the Kynp based on the mathematical model

The global sensitivity analysis EFAST was performed to evaluate the sensitivity of 35 parameters in the Kynp. The results, as shown in Table 1 and Figure 3F, G, presented the sensitive and insensitive parameters in the Kynp. Sensitive parameters were defined by S > 0.5 (a2 = 0.6516, a6 = 0.7622, a8 = 0.6172, a16 = 0.5870, a18 = 0.5260, b4 = 0.7760, b6 = 0.5710, b7 = 0.5320, b8 = 0.7180, b10 = 0.6950, and b11 = 0.5240). In order to confirm the sensitivity of the selected parameters, gene‐silencing experiments were performed. However, the sensitivity of the degradation parameters (b4, b6, b7, b8, b10, b11) has been ignored because we just focused on the expression level of the activation parameters (a2, a6, a8, a16, a18). Furthermore, the IDO1 produced by cells was mixed with the added IDO1 in the culture; thus, the a2 (IDO1) experimental test was ignored because it would affect the model calculation. Therefore, a6, a8, a16, and a18 were defined as sensitive molecules in the Kynp.

TABLE 1

The results of sensitivity analysis presenting all parameters

Parameter	a1	a2	a3	a4	a5	a6
Si	0.0355	0.0854	0.0364	0.0149	0.0636	0.1514
STi	0.3487	0.6516	0.2603	0.1909	0.4887	0.7622
Parameter	a7	a8	a9	a10	a11	a12
Si	0.0286	0.1465	0.0291	0.0286	0.0333	0.0317
STi	0.2803	0.6172	0.2944	0.1217	0.3292	0.2984
Parameter	a13	a14	a15	a16	a17	a18
Si	0.0474	0.0368	0.0392	0.0764	0.0188	0.0772
STi	0.4130	0.3524	0.3568	0.5872	0.1953	0.5260
Parameter	a19	a20	a21	a22	a23	b1
Si	0.0249	0.0352	0.0469	0.0524	0.0749	0.0424
STi	0.2324	0.2907	0.3934	0.3573	0.4360	0.3276
Parameter	b2	b3	b4	b5	b6	b7
Si	0.0570	0.0483	0.1424	0.0472	0.0761	0.1063
STi	0.4441	0.3952	0.7760	0.3701	0.5709	0.5315
Parameter	b8	b9	b10	b11	b12	b8
Si	0.0982	0.0492	0.1102	0.1176	0.0978	0.0982
STi	0.7178	0.3879	0.6948	0.5240	0.4105	0.7178

Bold values denote sensitive parameters; italic values denote insensitive parameters

Sensitive molecules in the Kynp showed significant effect on the expression of TGF‐β1

As shown in Figure 3G, parameter (a6) S = 0.7622 reflects the high sensitivity of activation rate from Kyn to GPR35. The knockdown of the Kyn gene, as presented in Figure 4A and B and Figure S4A, significantly reduced the expression of TGF‐β1 in comparison with the negative control (P = .0001) at the gene expression level and protein expression level, as seen in Figure 4I, J, K. Further, the activation rate of GPR35 to β‐catenin, NF‐kβp105 to ANK, and NF‐kβp65 to RelA showed sensitive values (a8 = 0.6172, a16 = 0.5870, a18 = 0.5260) that were also successfully silenced using siRNA transfection as shown in Figure 4C, D and Figure S4B for GPR35, Figure4E, F, and Figure S4C for NF‐kβp105 as well as Figure 4G, H and Figure S4D for NF‐kβp65. Similarly, knockdown of GPR35, NF‐kβp105, and NF‐kβp65 significantly reduced the expression level of TGF‐β1 in te PDAC (P < .0001, P < .0005, P < .0012), respectively (gene and protein expression levels are presented in Figure 4L‐N, P‐R, S‐U). It is also noticed that there is a negative relationship between sensitivity value and TGF‐β1 relative expression: Higher sensitivity showed a lower TGF‐β1 relative expression. Also, the comparison of the siRNA group with the normal group presented the same significant relationship between the selected sensitive parameters and TGF‐β1 relative expression (Figure 4I, L, P, S).

FIGURE 4

The biological activity of sensitive parameters of the Kynp in the pancreatic adenocarcinoma (PDAC) cell line was confirmed at mRNA and protein expression levels. Three different siRNAs were designed to test the most effective siRNA for the gene knocking down of calculated sensitive parameters in the Kynp in Panc‐1. A‐B, Protein expression of the knocked down Kyn gene in the Panc‐1 cell line. C‐D, Protein expression of the knocked down GPR35 gene in the Panc‐1 cell line. E‐F, Protein expression of the knocked down NF‐kβp105 gene in the Panc‐1 cell line. G‐H, Protein expression of the knocked down NF‐kβp65 gene in the Panc‐1 cell line. I, L, P, S The effect of the knockdown of Kyn, GPR35, NF‐kβp105, and NF‐kβp65, respectively, on the expression of TGF‐β1 in the Panc‐1 cell line: It showed significant reduction in the expression of TGF‐β1. J, M, Q, T, Protein expression level of TGF‐β1 in the case of kyn, GPR35, NF‐kβp105, and NF‐kβp65 knockdown, respectively. K, N, R, U, Normalized ratio of TGF‐β1 expression in the case of kyn, GPR35, NF‐kβp105, or NF‐kβp65 knockdown, respectively (error bars, mean ± SD; *P < .05, **P < .01, ***P < .001, ****P < .0001, Tukey's post hoc test)

Furthermore, we also verified the significant relationship between sensitive parameters (a6, a8, a16, a18) and TGF‐β1 protein expression in PDAC by Western blot (Figure 4J, M, Q, T). When Kyn was silenced, the TGF‐β1 protein expression ratio was significantly reduced (P < .01) in comparison with the NC and normal group (0.4 ± 0.07, 1.08 ± 0.05, 1.26 ± 0.019, respectively). Moreover, the silencing of GPR35 (P < .01), NF‐kβp105 (P < .001), and NF‐kβp65 (P < .01) significantly reduced the protein expression ratio of TGF‐β1 in the PDAC cell line in comparison with the NC and normal group (0.090 ± 0.01, 0.147 ± 0.01, 0.16 ± 0.02; 0.33 ± 0.02, 1.32 ± 0.04, 1.3 ± 0.07; 0.135 ± 0.04, 0.27 ± 0.04, 0.3 ± 0.008; respectively). Altogether, the experimental evidence interestingly proved the high sensitivity of the calculated sensitive parameters. This evidence strongly indicated the sensitive parameters are attractive, novel therapeutic targets in the Kynp that could contribute to the treatment of PDAC.

To confirm the significant effect of sensitive parameters, insensitive parameters were also tested by qPCR and Western blot to compare their effects on TGF‐β1 expression levels in the PDAC cell line. Three parameters with less than 0.5 (a7 = 0.2803, a9 = 0.2944, a19 = 0.2324) have been selected, which represent the mRNA expression level of β‐catenin, IKK‐β, and IL‐6 in the Kynp, respectively. The experimental analysis of these insensitive parameters showed that when β‐catenin was knocked down, the relative expression of TGF‐β1 showed no statistical differences in comparison with the NC group (P = .6378). Also, knockdown of IKK‐β or IL‐6 molecules presented no statistical differences in the expression of TGF‐β1 in the PDAC cell line (P = .7678 and P = .7099, respectively), as shown in Figure 5J, M, P. The selection of efficient siRNA and the confirmation of the downregulation of insensitive genes in PDAC are presented in Figure 5A‐I. The confirmation of insensitive parameter effects was performed by Western blot that showed no significant differences in the TGF‐β1 protein expression levels when β‐catenin, IKK‐β, or IL‐6 were silenced (P = .2612, P = .6636, and P = .3038, respectively), as presented in the Figure 5K‐L, N‐O, Q‐R. This means the insensitive parameters should be ignored in the therapeutic application trials against the Kynp in PDAC. These results also confirmed the effect of sensitive parameters on the expression level of TGF‐β1.

FIGURE 5

The biological activity of insensitive parameters in the Kynp in the pancreatic adenocarcinoma (PDAC) cell line was confirmed at mRNA and protein expression levels. Three different siRNAs were designed to test the most effective siRNA for the gene knocking down of calculated insensitive parameters in the Kynp in Panc‐1. A‐C, Results of gene relative expression of β‐catenin, Iκκβ, and IL‐6, respectively, in Panc‐1 under siRNA down regulation. D‐I, Performance of the selected siRNA in the knocking down of every tested gene. J, M, P, the TGF‐β1 gene expression results under the effect of β‐catenin, IKK‐β, or IL‐6 knockdown showed no statistical differences in the TGF‐β1 levels in comparison with the negative control group. K, N, Q, Protein expression levels of TGF‐β1 under the effect of β‐catenin, IKK‐β, or IL‐6 knockdown, respectively. L, O, R, Normalized TGF‐β1 expression ratio according to β‐Actin expression

Altogether, the established model importantly showed accuracy in predicting targeted molecules in the Kynp that suggested to be considered in advanced trials of PDAC treatment.

DISCUSSION

Pancreatic adenocarcinoma remains one of the highest‐mortality gastrointestinal tract cancers in the world. Up to date, there is no superior target to reduce PDAC mortality or improve survival rate. Conventional chemotherapy failed to provide any therapeutic improvements. Current research pays much attention to Kynp triggering in pancreatic cancer using IDO inhibitors because it presented a promising immunotherapeutic strategy against PDAC , (ClinicalTrials.gov Identifier: NCT02077881). However, the lack of sensitive targets in the Kynp is challenging in the treatment of PDAC in clinical trials. ,

In this study, we have established a mathematical model for accurately calculating the activity of Kynp molecules and predicting the sensitive parameters in the Kynp in PDAC. Recent statistical and mathematical models of PDAC have considered mortality database or disease risk factors. , However, our mathematical model introduces critical informative outputs describing Kynp molecules’ activity in PDAC based on ordinary differential equations, which is a type of more advanced modeling, as recently reported. , This model describes the dynamic expression trends of every parameter in the Kynp depending on concentration changes, which can not only contribute to deep understanding of Kynp activity but also forecast the trend of molecular activities and quantitatively analyze the molecular mechanisms of the Kynp regulating PDAC, as well as explore novel therapeutic targets.

By tracking the literature, twelve consecutive molecules involving the Kynp were included in the model establishment, starting with tryptophan amino acid, the main substrate of the Kynp; IDO1, the main mediator of Kynp; then GPR35; β‐catenin; , IKK‐β; NF‐κβ; , NF‐κβp65; IL‐6; and finally TGF‐β1, the essential biomarker of the Kynp in PDAC. , , , The correlation efficient of these molecules has been measured by Pearson’s correlation coefficient (R), which measures linear trends. It displayed >0.988 for all selected variables, which indicated the strong relationship between the established model and the expression trends of all variables in the Kynp. Moreover, the analysis of sensitivity explored significant sensitive targets in the Kynp (kynureninase, GPR35, NF‐κβp105, or NF‐κβp65) that can contribute to efficiently inhibiting the Kynp in PDAC. The significant effect of sensitive parameters was confirmed by testing the opposite insensitive parameters that showed no significant effects on the activity of the Kynp in the PDAC cell line. These findings support newly published studies that describe the increase of kynurenine metabolites in the serum as a risk factor for pancreatic cancer, GPR35 as a novel target in pancreatic cancer and other tumors, , and NF‐κβ as the most active player in pancreatic and human cancers. , , Further, our study evidenced that β‐catenin is an insensitive parameter in the Kynp in PDAC, which is in contrast with a published report that suggested that kynurenine activates β‐catenin to enhance breast cancer cell proliferation and metastasis. These contradictory findings could be ascribed to the genetic differences of every tumor cell line. Also, although, β‐catenin has no significant effect on the TGF‐β1 expression levels, it could participate in cancer cell proliferation indirectly. In addition, IKKβ and IL‐6 have been noticed to be insensitive in the Kynp, as they did not show a significant effect on the production of TGF‐β1 in the Panc‐1 cell line. Previous studies showed an unclear role of IL‐6 in PDAC. It could support cancer progression by playing a secondary role, which suggests a low significance of targeting the IL‐6 gene. Indeed, IKKβ could participate with other genes in the rapid proteosomal degradation leading to the release of NF‐κB and enhancing DNA binding, which could be mediated by other genes. Thus, the blocking of IKKβ is useless for inhibiting the activity of the Kynp in PDAC. The limitations of this study can be summarized as the follows: Using one kind of pancreatic cancer cells which is one of the most virulent types could be one of the limitations. However, it could work also with other pancreatic cell types, so this point needs further investigation in the future. Further, the prediction efficiency of the mathematical model in vivo has not been tested. As in vitro evaluation of gene expression is mostly different to in vivo evaluation, in vivo experiments are required in the future for evaluating this model accurately to explore the potential for in vivo prediction efficacy.

In conclusion, the activity of Kynp molecules as shown in Figure 6 have been elucidated by a newly established mathematical model. This is supposed to contribute to tracking, predicting, and modulating the biological activity of the Kynp. It could provide a promising tool to targeting PDAC through the downregulation of sensitive parameters. In addition, the established mathematical model could be applied to other cancers considering the Kynp with appropriate modifications.

FIGURE 6

Kynurenine pathway in the pancreatic adenocarcinoma. The schematic diagram shows the production and effects of all variables in the pathway that stimulate TGF‐β1 production as a marker of tumor cell survival. These variables have a production rate which is presented as X. It is an equation of predicted product during specific time. The model has been built according to the gene profiling of all variables presented in the figure

CONFLICT OF INTEREST

The authors declare that they have no conflict of interest. Our institution did not receive payment or services from a third party (government, commercial, private foundation, etc.) for any aspect of the submitted work. There are no relevant financial activities outside the submitted work.

Appendix S1

Click here for additional data file.