An individualized digital twin of a patient for transdermal fentanyl therapy for chronic pain management

Fentanyl transdermal therapy is a suitable treatment for moderate-to-severe cancer-related pain. The inter-individual variability of the patients leads to different therapy responses. This study aims to determine the effect of physiological features on the achieved pain relief. Therefore, a set of virtual patients was developed by using Markov chain Monte Carlo (MCMC) based on actual patient data. The members of this virtual population differ by age, weight, gender, and height. Tailored digital twins were developed using these correlated, individualized parameters to propose a personalized therapy for each patient. It was shown that patients of different ages, weights, and gender have significantly different fentanyl blood uptake, plasma fentanyl concentration, pain relief, and ventilation rate. In the digital twins, we included the virtual patients’ response to the treatment, namely, pain relief. Therefore, the digital twin was able to adjust the therapy in silico to have more efficient pain relief. By implementing digital-twin-assisted therapy, the average pain intensity decreased by 16% compared to conventional therapy. The median time without pain increased by 23 h over 72 h. Therefore, the digital twin can be successfully used in individual control of transdermal therapy to reach higher pain relief and maintain steady pain relief. Graphical Abstract (Created with BioRender.com)


Introduction
Fentanyl is a strong synthetic opioid that is 50 to 100 times more potent than morphine and 30 to 50 times more potent than heroin [1][2][3]. It has high lipophilicity, allowing it to quickly pass the blood-brain barrier [4]. Fentanyl transdermal patches are being extensively used for severe pain management, including cancer-induced pain [4,5]. Fentanyl therapy has a narrow therapeutic range [6]. Therefore, keeping the fentanyl concentration in plasma in this window is crucial to reach the therapeutic effect while avoiding its toxicity. Clinical results of fentanyl transdermal therapy show that a high inter-individual variability occurs in response to the therapy [7]. As such, there is a different uptake and response for each patient. Therefore, the therapeutic transdermal dose must be adapted for every patient. Due to the different physiological features of patients, their cause of pain, the intensity of the pain, and concomitant diseases, every patient needs a different dose of fentanyl. However, it is challenging to find the optimum therapy for each patient by trial and error. Such trial and error is discomforting and dangerous for the patient and requires a higher time commitment from the medical doctor.
By implementing the experimental and clinical data, and involving physics, physics-based modeling of the drug uptake can predict the patient's plasma concentration and clinical response to different therapies. Via this approach, we are able to examine different approaches to modify them based on the patient's needs in a shorter time without putting the patient in danger or discomfort. The aim of applying in silico study is providing additional and key information that is not possible or too complicated to gain by in vivo and in vitro studies. Numerous studies were conducted in order to develop a mathematical or computational model to predict the outcome of therapies for the patient [8][9][10][11][12][13]. These studies focus on different administration routes for medicinal products, such as oral, pulmonary, parenteral, and transdermal routes. These models for transdermal drug delivery are divided into two groups: (1) The studies focused on the drug's transport from the patch through the skin at different scales [14][15][16][17][18][19][20][21]. (2) The studies are focused on pharmacokinetics and pharmacodynamics modeling of fentanyl transport in the body [22][23][24][25]. However, simulations that cover the whole process of transdermal drug uptake, distribution and metabolization, clinical effect, and even patient feedback are rarer [26]. In order to reach this aim, in our previous studies, we developed a physics-based digital twin by including virtual patient feedback. The digital twin is a virtual presentation of the real-world patient that contains all his/hers organs and processes, and it is connected to the real world by sensor data or patient feedback. In this digital twin, we took to account the drug uptake by skin, pharmacokinetics, and pharmacodynamics model. Additionally, the feedback from the virtual patient was included in the digital twin. Every 24 h, the digital twin received virtual pain intensity feedback generated by a standard randomized number with 1 VAS unit standard deviation [26]. This digital twin was used to simulate transdermal fentanyl therapy for virtual patients of different ages. However, in our previous study, other patient physiology features besides age, such as weight and gender, were not considered, although based on the clinical studies, they strongly impact therapy outcomes.
In this study, a virtual population containing 3000 virtual cancer patients was generated based on sample data. This sample data was obtained from the available data in the literature [27], which included information on age, weight, height, gender, cancer diagnosis, and prior treatment based on the WHO steps of 20 cancer patients. The members of this population differ from each other by age, gender, weight, and height. Considering other impactful physiological features such as liver and kidney conditions will help to tailor the twin to the patients; however, we did not include them in this study due to the lack of sufficient data. In order to evaluate the outcome of fentanyl transdermal therapy for each patient, a physics-based digital twin was developed. The model parameters of this digital twin were modified based on patient physiology and considering inter-individual variability. The physics-based digital twin contains three main parts: (I) drug uptake model to predict the absorption of the drug from the patch through the skin to reach the blood circulation, (II) pharmacokinetics model to predict the fentanyl concentration in plasma, and (III) pharmacodynamics model to predict the two clinical effects of fentanyl, which are pain relief (therapeutic effect) and reduction in ventilation rate (adverse effect). This tailored digital twin proposes a therapy for each individual to relieve the pain efficiently and avoids hypoventilation based on the calculated pain level compared to the targeted pain intensity.

Drug uptake model
The drug uptake model simulates the drug diffusion in the transdermal patch and skin layers (stratum corneum, viable epidermis, and upper part of dermis). The thickness of the skin layer was modified for each patient, as this was dependent on patient physiology, such as gender, weight, and age. The geometric dimensions for the patch were based on Duragesic ® fentanyl patches. The Duragesic ® fentanyl patches provide a different flux of fentanyl to the body based on their surface area. The surface area of these patches varies between 5.25 and 42 cm 2 . In this model, the chosen section of skin and the patch have a square shape. We considered the skin block has the same length and width as the transdermal patch. We assumed that there was perfect contact between the patch and the skin. The implemented geometry in this study is shown in Fig. 1.

Computational system configuration
The governing equation The drug penetration from the patch through the skin was modeled by transient diffusion which is described by Fick's second law (Eq. 1).
where c i and D i are concentration [ng/ml] and diffusion coefficient [m 2 /s] of fentanyl in domain i . At the interface of each two layer, partitioning should be considered. The partition coefficient at each interface is mentioned in Eq. (2).
where K i∕j is the partition coefficient between the two domains i and j . The partition coefficient is the concentration ratio of one solute at the interface of two immiscible phases. Due to the partition coefficient, the concentration of fentanyl at the two sides of the interface of the layers might be different. This discontinuity is inconvenient to solve computationally. To solve this problem, the drug's potential i , as described in Eq. (3), it was used instead of the concentration, as it is continuous throughout the layers.
where K i is drug capacity in domain i, and it is related to the partition coefficient as: K i∕j = K i ∕K j (as a result of the continuity of i over the whole domain).
Fentanyl penetrates through the epidermis and reaches the dermis. In the dermis layer, the drug will be taken up by capillary vessels at different lengths of the dermis. Besides diffusion, other physical mechanisms, such as advection, will play a role. In this study, we only modeled the diffusion process. Additionally, the capillaries will absorb the fentanyl in the upper part of the dermis, and all of the drug or part of it will be absorbed before reaching its bottom. To obtain the same drug uptake through the skin as the clinical studies, we used an equivalent dermis thickness, which is smaller than the real dermis thickness. The details of this assumption and modification and validation of drug uptake model are mentioned in our previous studies [26,28,29].
Boundary and initial conditions At the top boundaries ( Fig. 1, b1) and for the peripheral boundaries ( Fig. 1, b2 (3) c i = K i × i and b4), we assumed no flux of fentanyl. In this study, the considered patch contains 12.6 mg of fentanyl ( Fig. 1, d1). At the bottom boundary ( Fig. 1, b3), the concentration of fentanyl is taken equal to the concentration of fentanyl in blood. This concentration is obtained from the pharmacokinetic model. The initial concentration of fentanyl in the skin layers ( Fig. 1, d2, d3, and d4) is equal to zero, as there is no drug in the skin when applying the patch. This assumption implies that the patch is always applied at a new location.

Pharmacokinetics modeling
In this study, a physiologically based pharmacokinetic model was used to predict the concentration of fentanyl in the blood plasma. Due to the complexity of the human body, different organs were lumped together into PK compartments [26]. We considered five compartments: (I) the central compartment, which contains blood and lungs (Eq. 4); (II) the rapid equilibrated compartment, which includes the brain, heart, skin, and kidneys (Eq. 5); (III) the slow equilibrated compartment, which contains muscle, fat tissue, and the carcass (Eq. 6); (IV) gastrointestinal compartment which includes gut, spleen, and pancreas (Eq. 7); and (V) hepatic compartment where the metabolization happens and which only includes the liver (Eq. 8). The concentration of fentanyl in these compartments is evaluated by the following ordinary differential equations: The parameters of these equations are mentioned in Table 2. In this set of equations, c p , c r , c s , c g , and c l are the concentration of fentanyl in the central, rapid equilibrated, slow equilibrated, gastrointestinal, and hepatic compartments, respectively. f u is the unbound fraction of fentanyl in the central compartment. k ij , k met , and k re are the first-order equilibrium rate constants for inter-compartmental clearance, metabolization, and renal clearance, respectively. In our previous study, the pharmacokinetic model was validated [26].
Implemented geometry of skin and patch for drug uptake model; boundaries b1, b2, and b4 have no flux; however, fentanyl is absorbed by the blood circulation system in the b3 boundary. In this model, the patch's thickness is 50.8 µm, the stratum corneum's thickness varies in the range of 12-27 µm, the thickness of the viable epidermis is between 15 and 51 µm, and the thickness of the upper part dermis is between 153 and 374 µm

Pharmacodynamics model
There is a delay between the concentration of fentanyl in the blood plasma and the corresponding clinical effects on pain relief and breathing rate [30]. To account for this delay, a virtual compartment is defined, called the effect compartment. The concentration of fentanyl in the effect compartment related to the central compartment is calculated in Eq. (9).
As fentanyl is a small lipid-soluble drug, it will rapidly penetrate through the blood-brain barrier from the blood circulation system and bind to opioid receptors in the central nervous system [31]. Through this mechanism, fentanyl produces pain relief. We used the visual analog scale (VAS), which is a pain scale between 0 (indicates no pain) and 10 (indicates the worst possible pain), to evaluate pain relief. Additionally, fentanyl will bind to the opioid receptor on chemoreceptors in the brain stem, which are responsible for detecting the carbon dioxide level in the blood and stimulating breathing. This binding will reduce the sensitivity of carbon dioxide levels and may lead to stopping breathing in extreme cases. A sigmoidal model models the reduction of the ventilation rate and VAS pain score as a function of fentanyl concentration in the effect compartment (Eq. 10).
where E i , E 0 i , and E max i are the intensity of the pharmacologic, baseline, and maximum reachable effect, respectively. EC 50,i is the concentration related to the response halfway between the baseline and maximum possible effect, and i is the Hill coefficient. In our previous study, the parameters for the pharmacodynamics model for pain relief were derived based on experimental data [26].

Sample data of cancer patients
Ten women and ten men with age 40 to 68 with different tumor sites were studied by Zech et al. [27]. All the involved patients used fentanyl transdermal therapy for chronic pain management. The patient's physiological features are mentioned in Fig. 2, including age, gender, weight, and height. These parameters are cross-correlated. To explore the correlation between these physiology features, the correlation test is implemented on this set of data, and its result is shown in Fig. 2. This result shows a strong correlation (correlation coefficient (C.C.) = 0.87) between gender and height. However, it should be clarified that the value for gender is considered Boolean in this set of data. The value for the male gender is one, and for the female, gender is zero. Therefore, the positive correlation coefficient of 0.87 represents the higher height between male genders, as is expected. The correlation test result also represents a weak correlation between gender and weight (C.C. = 0.43) and weight and height (C.C. = 0.38).

Virtual population
In this study, we developed a virtual population of patients from sample data to examine the performance of digital twins for patients with a wider range of physiological features. Members of this population differ from each other based on their age, gender, weight, and height. In reality, these physiological features of patients are correlated, as shown in Fig. 2. Therefore, our virtual patients should have correlated features in the same way. We used the Markov chain Monte Carlo (MCMC) method to produce the virtual population with realistic feature combinations for each patient. MCMC predicts the posterior distribution by sampling from the joint distribution of likelihood and prior distribution. This method uses the sample demographic data of cancer patients (N = 20) as a likelihood function. The prior distribution for age, weight, and height is considered normal, and gender is considered a uniform distribution. By using the Gibbs sampling method, the posterior distribution (N burn = 20,000 and N keep = 3000) was evaluated. The demographic characterization of virtual patients of the generated population (N = 3000) was obtained from this joint distribution. In Fig. 3, the generated virtual population is shown. By comparing the correlation coefficients between predicted features for the virtual population and sample data, we find that the MCMC method was able to mimic the correlation between sample data. However, the distribution of parameters in the sample data and population data are visually different due to the limited number of cases in the sample data. The virtual population members are chosen based on sample data to reach a population with physiological features compatible with sample data of real patients. The summary of sample data and population data is shown in Table 1.

Estimation of model parameters
The patient's physiology changes the skin properties, volume of the organs, blood flow, glomerular filtration rate, metabolization, and opioid tolerance of the patient. By using the patient's individual physiological features, we can predict the input parameters for the model. It should be noted that despite these equations being obtained from experimental data, as a result of inter-individual variability, different values for these parameters correspond to patients with the same age, gender, weight, and height. To account for this inter-individual variability, the calculated value was modified by the exponential of a normal random number with mean = 0 and SD = 0.1. By including this random factor, 95% of the evaluated parameters are in the range of ± 20% of the evaluated values. The overall dependency of the model parameter on patient physiology is mentioned in Table 2. As mentioned earlier, if the gender of the patient is male, the value for the gender parameter is equal to 1 and for the female is 0.
Based on Table 2, skin layer thickness, compartment volumes, rate constants, clearance rate, and the concentration of half maximum effects for pain relief are directly calculated based on physiological features. The rest of the model parameters in this study are considered independent of the patient's physiological features and are in Table 3. In real life, these parameters might depend on the patient's physiology. However, for the four features considered in this study, based on our best knowledge, there is no proposed experimental model to calculate these parameters based on these four physiological features. By considering both parameters, dependent and independent of patient physiology, all input parameters change between the patients.

Spatial and temporal discretization
The grid sensitivity analysis by using Richardson extrapolation was done over the skin layers and fentanyl patch based on the flux out of the dermis, in which spatial discretization error was considered 0.1%. As a result of this grid sensitivity analysis, about 900 quadrilateral grids were built. However, this number varies from case to case as the thickness of the skin's layer differs between cases. In some cases, almost 100 quadrilateral grids were sufficient. In order to increase the numerical accuracy, the accumulation of grids near the interfaces is more. To define the time steps, sensitivity analysis based on maximum fentanyl concentration was done, which defined a maximum time step of 6 h; however, to reach a higher temporal resolution, time steps were assigned as 1 h.

Numerical implementation and simulation
COMSOL Multiphysics version 5.4 was used in this study to solve the diffusion process of fentanyl from the patch through the skin in the mechanistic model, the distribution of fentanyl in the human body in the pharmacokinetics model, and the drug's effect in the pharmacodynamics model. The MUMPS (MUltifrontal Massively Parallel sparse direct Solver) was chosen as the solver scheme in this set of simulations. The partial differential equation (PDE) interface was implemented to solve the diffusion process of fentanyl in the mechanistic model. To take to account the distribution, elimination, and metabolization of fentanyl, the boundary ordinary differential equation (ODE) was used. In the pharmacodynamics model, the concentration of fentanyl in the effect compartment was calculated by the PDE interface, and the boundary probe calculated the drug's effect. To apply the change of the patch and skin location during the therapy, the event interface was used. The population generation was done in RStudio by using the "mixAK" package. Analyzing the sample data, calculating the posterior distribution, generating the virtual patients' characteristics, calculating the model parameters, and analyzing the result of digital twins are done in RStudio.

Proposed therapy by the digital twin
The tailored physics-based digital twin predicts the outcome of fentanyl transdermal therapy for the patient at any moment of the treatment. Based on the summary of product characteristics (SmPC), which we will call "Conventional therapy" in this study, the fentanyl transdermal patch needs to be changed every 72 h. After applying the patch, fentanyl accumulates in a depot in the skin, and it is gradually released into the blood circulation, and after reaching its peak in blood, the concentration in the blood decreases as the amount of drug in the patch decreases. As a result of the lower concentration, the analgesic effect decreases too. Patients who do not reach sufficient analgesic effect within 72-h intervals can switch to 48-h intervals [43]. This prescription is similar for every patient as there is a limited number of alternative approaches [26]. However, by using the predicting capability of the digital twin, we can propose a dynamic application time for each patch based on patient needs. In this regard, every 8 h, the digital twin considers the condition for changing the patch. If the pain intensity is above target (in this case, VAS pain score equal to 3), the patient needs to change the patch. The reason for putting a limit of 8 h is to avoid changing the patch too frequently, which may lead to an intensive increase in concentration in plasma and drug waste. Additionally, there is a 1 to 2 h of time lag between the application of the fentanyl transdermal patch and reaching the blood circulation [44]. Based on our previous study, the time to reach half maximum concentration is about 8 h [26]. In this way, the digital twin gives enough time to evaluate the patient's pain relief. The general approach of digital-twin-assisted therapy is shown in Fig. 4.

Conventional therapy
The conventional therapy was applied to all virtual patients over 3 days. Three thousand virtual patients implemented a Duragesic ® fentanyl patch with a nominal flux of 75 µg h −1 over 72 h based on the SmPC. The outcome of therapy was evaluated based on three output parameters: fentanyl concentration in plasma, pain intensity, and ventilation rate. This result is shown in Fig. 5. In Fig. 5a-c, the outcome based on patient age is represented, which shows no distinct trend in the effect of age on fentanyl concentration in plasma, pain intensity, or ventilation. This implies inter-individual variability, and the effect of characteristic features of the patient overshadows the effect of age. However, our previous study showed that the age of the patient drastically affects the therapy's outcome [26]. By comparing the results of these two studies, we can conclude that the impact of age on fentanyl therapy outcomes is considerable when the other physiological features are similar between patients. However, the impact of other physiological features, such as weight, overshadows the impact of age. Therefore, due to the high impact of other characteristics of the patient, with solely the patient's age, we cannot predict the therapy outcome. Figure 5d-f show the outcome of therapy versus weight. Based on the result, by increasing weight, the concentration in plasma drops, which leads to lower pain relief and simultaneously lower side effects, which is the reduction in ventilation rate. The reduction of fentanyl concentration in plasma by increasing weight is due to the larger volume of each compartment; the drug will be diluted. Therefore, weight has a key impact on pain relief, which is the main aim of the therapy. In the next step, we studied the effect of height on the outcome of the therapy. This result is demonstrated in Fig. 5g-i, which shows that the concentration in plasma reduces by increasing height. As a result of the reduction in plasma concentration, pain relief reduces as well. However, based on Figs. 2 and 3, higher height is correlated to higher weight, and the reduction in plasma concentration is due to higher weight, not necessarily higher height. Here in these three graphs, we see a gap between genders due to the height gap between the genders. Nevertheless, standard fentanyl therapy clearly is not gender-neutral. This deviation is partially due to different height distributions between genders; however, besides the height impact, the gender of the patient directly impacts the skin layer thickness and the volume of PK compartments. Therefore, by considering similar age, weight, and height for the two genders, we still will end up with different fentanyl therapy outcomes.

Digital-twin-assisted therapy
As mentioned in "Proposed therapy by the digital twin" section, digital twins predict the pain intensity and how the treatment influences it. Based on the intensity of the pain, the digital twin can propose changing the patch if the pain intensity is higher than the target (VAS pain score 3). The digital twin is able to calculate the pain intensity at any moment and update the therapy based on that. However, based on our previous study, it takes about 8 h for a fentanyl patch to reach half of the maximum contraction in plasma. Therefore, we put a limit of 8 h for the checking points to update the therapy. At first, for patient number one (age = 58 years, male, weight = 74.2 [kg], and height = 1.73[m]), we explored the effect of different checking points at every 2, 4, 6, 8, 12, 24, 36, and 72 h. We did not consider frequencies shorter than 2 h, as the time lag of fentanyl in reaching the blood circulation is about 2 h [44]. In Fig. 6a, the profile of pain relief for different checking points is shown, and based on Fig. 6b, except for the checking point every 24 and 72 h, the average pain intensity is lower than 3. The checking point at every 36 h has better pain relief than 24 h because this specific patient has higher pain intensity around 36 h after therapy. Additionally, based on Fig. 6c, for all the checking points duration every 24 and 72 h, the time without pain is more than 48 h (more than 2/3 of the therapy duration). As shown in Fig. 6d, except for the checking point every 72 h, all other cases pass the plasma concentration threshold; however, the criteria for breathing rate are met by all the cases (shown in Fig. 6e). In Fig. 6f, the number of the needed patch by the first patient, by considering different checking points. Considering the pain relief performance while avoiding excessive  The concentration of fentanyl in plasma for all the virtual patients during conventional therapy versus time is shown in Fig. 7a. As shown in the graph for a large number of patients, the fentanyl concentration in plasma reaches the upper threshold during the therapy. It should be noted that the threshold of 2 ng ml −1 [45] is an average for all patients. This number in other studies is 3 ng ml −1 [46]. This threshold is

242.6μm
Modified equation [33][34][35] PK model Compartment volume -Modified equation [36,37] First-order equilibrium rate constant  ×exp (θ k ) defined based on the fentanyl concentration, which leads to toxicity for the patient. By considering inter-individual variability between patients and different conditions in different studies, the measured fentanyl concentration corresponding to toxicity might vary. As the concentration of fentanyl in plasma after reaching the maximum reduces, the patient's pain intensity starts to increase as well, as shown in Fig. 7b. More than 15% of the patients experience minute ventilation below 4 L/min during the treatment, as shown in Fig. 7c. In the next step, we implemented digital-twin-assisted therapy for each patient. The outcome of digital-twin-assisted therapy is shown in Fig. 7d-f. By comparing the fentanyl concentration in plasma for digital-twin-assisted therapy (as it is shown in Fig. 7d) with conventional therapy, we realized the average concentration in plasma increased by 11.5%, while the standard deviation decreased. This increase in concentration was necessary to keep the pain intensity below target. However, via this method, the fluctuation in concentration reduces. The result in Fig. 7e shows that the digital twin successfully kept the pain intensity for 98.8% of the patients below 3 for more than half of the duration of the treatment. However, for conventional therapy, only 57.1% of the patients have pain intensity below 3 VAS for more than half of the treatment duration. By implementing digital-twin-assisted therapy, the pain intensity decreased by 16% compared to conventional therapy, and the standard deviation decreased. A lower standard deviation corresponds to a lower deviation for pain intensity between patients. The higher fentanyl concentration in plasma reduced the average minute ventilation by 15% compared to conventional therapy. Besides reducing the overall pain intensity, the aim of digital-twin-assisted therapy is to keep the pain below the target. In this study, we defined a parameter for the duration that the patient experiences pain lower than VAS pain score 3, called the time without pain. By implementing digital-twin-assisted therapy, the median time without pain for the population increases by 23 h compared to conventional therapy, which is shown in Fig. 8a. The interquartile range (IQR) of time without pain for conventional therapy is 23 h, while it is 7 h (70% less) for digital-twin-assisted therapy. This shows that the variation of therapy's outcome for digital-twin-assisted therapy is less than conventional therapy. In the next step, we studied the effect of patient physiology in time without pain. In Fig. 8b and f, the effect of gender on the duration that the patient has lower pain intensity is shown. During conventional therapy, the difference between the median time without pain for the male patients and the female patients is 5 h, while this number is 40% less (3 h) during digital-twin-assisted therapy. This shows with proposed therapy by digital twin, the effect of gender on the outcome of therapy will reduce. The time without pain corresponds to the patient's age, weight, and height both for conventional and digital-twin-assisted therapy, as shown in Fig. 8c-e and g-i. The correlation coefficient between time without pain and patient age through both therapies is 0.1. However, the correlation coefficient of time without pain through conventional therapy and digitaltwin-assisted therapy and weight is − 0.41 and − 0.31, respectively. This implies that age does not affect the outcome of therapy, similar to conventional therapy. However, during digital-twin-assisted therapy, the outcome of therapy is less dependent on the patient's weight compared to conventional therapy. This is evaluated for the studied range of weight for virtual patients (42-107 kg), and lower or higher weights might change the results. The correlation coefficient between time without pain and patient's heights through conventional therapy and digital-twin-assisted therapy is − 0.15 and − 0.23. However, as mentioned earlier, the patient's height affects the model parameters via its effect on the BMI. Therefore, calculating the correlation coefficient for BMI will be − 0.30 and − 0.14 for conventional therapy and digital-twin-assisted therapy, respectively. By considering the result in Figs. 7 and 8, we conclude that digital-twinassisted therapy successfully increased the pain relief for the  Average weight and age of the patients, based on the number of patches they used in 72 h during digital-twin-assisted therapy; the blue area corresponds to the patients that they needed to change the patch more frequently patient and decreased the dependency of therapy outcome on patient physiology. With digital-twin-assisted therapy, we are able to provide the same efficiency of pain relief for men and women, old and young patients, and patients with higher or lower BMI.
Through digital-twin-assisted therapy, each patient will use an individualized number of patches for 72 h based on their pain intensity experience. Figure 9 presents the number of used patches during 72 h based on patient weight. Based on this result, patients with higher weights, on average, will need to change the patch more frequently compared to lower weights to achieve constant pain relief. As also mentioned earlier, age does not influence the number of needed patches. By analyzing the data, we realize that the majority of the patients who need to change the patch considerably and more frequently are male (shown in the blue circle). This is due to the fact that, on average, male patients are taller; therefore, they tend to have a higher weight. Therefore, they proportionally need more fentanyl to keep the pain intensity below the target.
In the next step, we analyzed the patients that needed to change the patch more frequently. We separated the part of the population that they need more than one patch per day. We compare the patient physiology distribution between the population and this subgroup in a density graph. The result in Fig. 10 shows that the two groups' age distribution is not significantly different (p value = 0.45). However, for weight distribution, there is a shift in the mean of the distribution (by + 6 kg), and with a p value of 2.2 * 10 −5 , these two distributions are different. For height distribution, the mean of the distribution is shifted by 3 cm, and the two distributions are significantly different (p value = 0.004). The ratio of females to males in the population is 1.02; however, in the subgroup, this ratio is 0.29. Based on this result, patients of the male gender and with higher body mass need to change the fentanyl patch more frequently to reach the target pain relief.

Outlook
This study analyzed the patient characteristics that affect fentanyl transdermal therapy. These characteristics were age, weight, gender, and height. Based on patient physiological features and published experimental studies, model parameters for drug uptake, pharmacokinetics, and pharmacodynamics model were calculated. In the next step, the tailored digital twin for each patient proposed a modification of the therapy based on the pain intensity. Digital-twin-assisted therapy successfully reduced the pain intensity and the variation in therapy outcomes for virtual patients. However, there are approaches to increase the compatibility to the realworld and higher efficiency and safety in the therapy. These approaches are mentioned below: -In order to reduce the error in the calculation and prediction of the digital twin, the model parameters should be calculated based on the outcome of fentanyl therapy on the patient, such as fentanyl concentration in plasma and pain intensity. -Liver enzyme activity and kidney filtration rate greatly impact the plasma concentration level as they determine the elimination and metabolism rate of fentanyl [47]. These physiological features must be considered in tailoring the digital twin for the patient. -In this study, only pain intensity was considered the main factor in modifying the therapy. However, in order to reach a safer therapy, it is important to consider adverse effects such as hypoventilation in decision-making. -Besides respiratory depression, fentanyl transdermal therapy can cause other side effects such as nausea, skin dryness, muscle spasms, confusion, and anxiety [43]. Taking to account these side effects and trying to avoid them while having efficient pain relief will lead to more tolerable therapy for the patient.

Conclusions
In this study, we tailored a physics-based digital twin of transdermal fentanyl therapy for 3000 virtual patients. By implementing a set of sample data and Markov chain Monte Carlo (MCMC), these virtual patients that were different from each other based on age, gender, weight, and height were developed. The model parameters for each Fig. 10 a Age distribution, b weight distribution, c height distribution, d gender percentage for the population and the subgroup of patients needing to use more than one patch a day patient were calculated based on these physiological features and applied inter-individual variability. Conventional SmPC-based therapy of fentanyl transdermal patch, which is implementing fentanyl patch and replacing it every 72 h, was studied. Based on the outcome of conventional therapy, the maximum fentanyl concentration varied in the range of 1.3-3.2 ng ml −1 . Additionally, the minimum VAS pain score varied in the range of 0-5, and the minimum minute ventilation rate was in the range of 2.5-10 L min −1 for the virtual patients. This result shows that the outcome of therapy has a high variation between patients as some patients receive excessive amounts of fentanyl, while some of the virtual patients could not reach the target pain relief. In order to reach the target pain relief for each patient, we designed a new therapy by digital twin. In this digitaltwin-assisted therapy, the patient's pain intensity will be checked every 8 h, and if it is higher than the target, the patient implements a new patch of the same size. As a result of modifying the treatment based on the patient feedback, the average pain intensity decreased by 16%, and the standard deviation among the population decreased by 54%. Based on the result, the digital twin successfully kept the pain intensity for 98.8% of the patients below 3 for more than half of the treatment duration. However, this number is only 57.1% for conventional therapy. Digital-twin-assisted therapy increased the median time without pain (VAS pain score under 3) by 23 h in 3 days and reduced its interquartile range (IQR) by 70%. Based on this result, digital-twinassisted therapy reduces the overall pain intensity for the whole virtual population and reduces the variation of the outcome of therapy between the patients. Therefore, the tailored digital twin can predict the patient's needs and conditions, propose a therapy, monitor the endpoint of therapy, and update when the patient needs it. In this way, all patients suitable for fentanyl transdermal therapy can receive suitable pain relief.