**4. Results**

On Figure 4, one can see two efficiency lines for CO2; upper one (blue) represents the efficiency of the Rankine-like subcritical CO2 power cycle, while the lower one (green) represents its TFC counterpart. The minimal temperature (condenser temperature) were chosen as the triple-point temperature of the CO2 (216.6 K), and the efficiencies were calculated up to the critical point (304.13 K). On Figure 4, the whole curve can be seen, while on Figure 4b, the high-temperature region can be seen, showing the well-developed maximum on the upper curve, 4 degrees below the critical point.

**Figure 4.** The efficiency of a CO2 Rankine-like Power Cycle (upper lines) and Trilateral Flash Cycle (lower lines) with fixed condenser temperature (216.6 K) and with increasing evaporator temperature. (**a**) full temperature range; (**b**) the vicinity of the critical temperature.

On Figure 5, a set of these curves (approximately 20 for each material) can be seen for all materials listed in Table 1. The first sign of this maximum was found earlier in butane [26]. For water (5a) and carbon dioxide (5b) (two typical wet working fluids), the maxima on all ORC curves are very characteristic. For CO2, red dots are marking these maxima for better visualization; these points are omitted for the other materials to avoid them to distort the shape of the curves. The increase of the maximal cycle temperature (with a fixed minimal cycle temperature) causes an increase of the efficiency for both cycles (the Rankine and the flash-types, represented by the upper and lower curves). Reaching higher maximal temperatures, the speed of this increase will be lower and lower, and finally a maximum, then a decreasing part can be seen. It means that with a fixed condenser temperature, the increase of maximal temperature has an initially positive effect of the first-law efficiency, but when it is reaching the vicinity of the critical point, this positive effect disappears, and the further increase would be contra-productive. On the contrary, for the TFC-like cycles, the increase holds up to the critical point, i.e., having a TFC-like cycle with a wet working fluid, the increase of maximal cycle temperature always increases the efficiency. In this way, the two curves are forming an elongated rain-drop form; the pointed end is located at the condenser temperature (with *η* = 0). Increasing the condenser temperature, one can obtain smaller efficiencies; both the upper and lower curves are shifting down, while the location of the maximum (with smaller and smaller value) shift closer to the critical temperature (see the location of the red dots for CO2), but always remains below that value. The rain-drop shape formed by the upper and lower curves remains, but it will be smaller and smaller. For the sake of the better visibility, the location of the maxima are marked by red dots for the CO2 curves (Figure 5b); its movement can be clearly seen.

Hexafluoroethane (R116) is originally an almost isentropic working fluid. The *T-s* diagram has an inverse S-shaped saturated vapor part, which is very narrow, i.e., the entropy value for the vapor phase (except in the vicinity of the critical temperature) is almost constant. In the novel classification system, it is a type ACNMZ. From the mapping of various working fluid types [21], this is an almost wet isentropic fluid. The point M located at 253.5 K, i.e., choosing condenser temperature above this value, R116 acts like a wet working fluid. Point N is located at 220.9 K, i.e., choosing condenser temperature between 220.9 K and 253.5 K, this fluid can mimic a moderately dry working fluid. Efficiency-curves can be seen in Figure 5c. The curves are very similar to the ones seen for water and CO2; the maxima are still well-developed (although these points are closer to the critical point than for the two wet ones), and the joint ORC-TFC curves are also raindrop-shaped. For this material, one can also say that increasing the maximal cycle temperature has a positive effect on the first-law efficiency, but in the vicinity of the critical point, this effect disappears and changes to negative.

**Figure 5.** Temperature dependencies of ORC (or other Rankine-like cycles) and TFC cycles for various working fluids (Water (**a**), CO2 (**b**), R116 (**c**), Butane (**d**), Neopentane (**e**), Dodecane (**f**)); see text for further explanation.

Butane is a moderately dry working fluid; it means that although the *T-s* diagram is not bell-shaped, but tilted, the titling is not very strong. ORC-TFC efficiency curves can be seen in Figure 5d. Although maxima can be still noticed, they are remarkably less developed than for the previous three materials.

For neopentane (Figure 5e), which has a more tilted *T-s* diagram that butane (i.e., it is drier), the rain-drop shape is already very distorted, and even the flat end is almost pointed (although the transition from ORC to TFC is still smooth). Very small, nearly diminished maxima can be seen, located very close to the critical temperature. In this case, one can assume, that increase of the maximal cycle temperature is always good or almost neutral for the efficiency. For the uppermost curve, the location of the maximum is only 1.97 K from the critical temperature, and even above this point, the decrease of efficiency is very small, so, one can consider it almost temperature-independent.

Finally, a very dry (on *T-s* diagram, represented by a very tilted *T-s* curve) material, dodecane (Figure 5f) were analyzed. One can see, that at low evaporation temperatures, ORC efficiency exceed TFC one, as it happened before, but increasing this temperature further, the difference will be smaller and smaller, and at a given temperature (which depend on the condenser temperature, i.e., different for each pair of curves) it disappears at the crossing of the two curves. Passing that crossing temperature, the efficiency of TFC is higher than for ORC; finally, they will join smoothly at the critical point. In this case, the efficiencies of TFC can form a shallow maximum, close to the critical curve, while the efficiency of ORC increases continuously. Therefore one can conclude, that by using very dry working fluid for an ideal (strictly isobaric heat exchanges and strictly isentropic expansion/compression) basic-layout (recuperator and superheater-free) ORC-system, the increase of maximal cycle temperature can always increase cycle efficiency, even in the vicinity of the critical temperature; additionally, small maxima might be expected for the efficiency of TFC.

3-D version of these lines for water (Figure 6a–d) and butane (Figure 6e–h) are also shown Rotating animated gif version of these figures are provided as Supplementary Material. The distortion of the initial rain-drop shape can be clearly seen.

**Figure 6.** *Cont*.

**Figure 6.** 3-D representation of the efficiency-curves for water (**a**–**d**) and butane (**e**–**h**). Rotating version of these figures is provided in the Supplementary Material as Figure S1a (water) and Figure S1b (butane) in animated gif-format.

#### **5. Theoretical Reason for the Existence and Disappearance of the Maximum**

A simple geometrical explanation can be given for the appearance of efficiency maximum in wet working fluids. It can be seen in Figure 7 that the ORC or any other Rankinelike cycle (steam Rankine or CO2 power cycles) can be seen as a superposition of a TFC-like and a Carnot-like cycles, although this superposition is not simply a summing of efficiencies. The TFC-like part is the one at lower entropies, resembling a triangle, while the Carnot-part (located on higher entropies) is the rectangular one. Considering that the Carnot-cycle is the best potential cycle between the given maximal and minimal temperature, adding it to a TFC, the Carnot-like part always has a positive contribution to efficiency. In this way, for wet working fluids, *η*ORC is always bigger, than *η*TFC. One can see, that by increasing the maximal cycle temperature (compare Figure 7a,b), the rectangular part will be smaller

and smaller, almost diminishing in the vicinity of the critical temperature (Figure 7c), and when *T*max = *T*cr, it disappears, and the TFC and ORC cycles will be indistinguishable. Therefore, at *T*max = *T*cr, *η*ORC = *η*TFC, and the transition will be smooth and continuous.

**Figure 7.** Rankine-cycle (black) in a wet working fluid as a superposition (but not the simple sum) of a TFC-like (red) and Carnot-like (green) parts with same condenser temperature with a medium (**a**) high (**b**) and almost critical (**c**) evaporator temperature.

This smoothness and continuity affect the shape of the efficiency-temperature curve. Previously we used *T-η* curves, here, for a short time, the inverse should be used, the *η-T* one (see Figure 8). Being *η*ORC > *η*TFC for wet working fluids at every temperature, except *T*cr, the low-entropy part (left side) represents the TFC, and the high-entropy part (right side) represents the ORC. In this representation, the two curves have to have end-points at the maximal temperature (*T*cr). Continuity causes the equality of the two efficiencies at *T*cr, i.e., the two curves have to join in this point, i.e., instead of two end-points, they will have one common end-point. Due to the smoothness, the derivative d*T*/d*η* will be zero in this point; therefore, this will be a smooth maximum with positive d*T*/d*η* on the left side and negative d*T*/d*η* on the right one. It should be recalled, that the left side represents the efficiency of ORC; being d*T*/d*η* positive on the original, *T*-*η* representation, this part would show negative slope. Therefore, in the vicinity of *T*cr, the smooth and continuous transition between ORC and TFC cycles would cause a part, where efficiency decreases with increasing temperature. Contrary, at lower temperatures, one can see an increasing part; this can be possible only by assuming a maximum on the *η*(*T*) curve of the ORC.

One should remember that a strong point of this reasoning was to handle ORC as a superposition of a TFC-like and a Carnot-like part. While this is true for the wet working fluid (Figure 9a), in a dry one, a third part, resembling a distorted, upside-down, inverse TFC, can also be defined, see Figure 9b–d. This part is also less efficient, that a Carnot-like one. Increasing the "dryness", i.e., having more and more tilted *T-s* diagrams and going from ACZM to AZCM class [23], the weight of the middle Carnot-like part will be smaller and smaller. Finally, this part can completely disappear.

**Figure 8.** Schematic temperature-efficiency (**a**) and its inverse, efficiency-temperature (**b**) curves of an ORC-like and TFC-like system, explaining the necessity of the maximum. Further explanation can be found in the text.

**Figure 9.** Different structure of an ORC cycle (black) with a wet, moderately dry, a medium dry and a very dry working fluid, divided into TFC-like (red), Carnot-like (green) and inverse TFC-like (grey) sub-cycles. The diminishment of the Carnot-like part (green) with increasing dryness can be clearly seen, as well as the loss of triangularity of TFC (red) and the appearance of the inverse TFC-like type (grey). The four figures are representing wet (**a**), moderately dry (**b**), medium dry (**c**) and very dry (**d**) working fluids.

This process can be seen in detail in Figure 9a–c. With wet working fluid (a), the ORC can be seen as a superposition of an almost perfectly triangular TFC (red) and an ideal Carnot (green) cycles, as it has been shown before (Figure 7). For moderately dry one (b), a third part (grey) appears, located at higher entropies, resembling a distorted upsidedown TFC. The appearance of this part makes the Carnot-like part narrower, decreasing

its superficial contribution to net first-law efficiency. For medium dry one (c) the isobaric, but not isothermal part of the heat removal is already so significant, that the Carnot-like part totally disappears and the cycle will be the superposition of an almost triangular TFClike and an upside-down TFC-like part; because the Carnot-part disappeared, one might expect that net efficiency might be smaller than previously. Finally, with a very dry working fluid, the isobaric heat removal part preceding the condensation will be so huge, that it will reach the initial TFC-like part, cutting down the corner and terminating the "triangularity" (Figure 9d).

Because this transition happens with increasing dryness, one can expect to see the similar process by going from wet (water, CO2) through isentropic (R116), moderately dry (butane) and dry (neopentane) working fluids to very dry one (dodecane). Upon this transition, the maximum will be less and less significant, and finally, it can disappear, as it is shown for dodecane.

#### **6. Discussion and Conclusions**

The increase of the maximal cycle temperature (which is equal to the evaporator temperature in ideal, simple, superheater-free systems) is thought to be a good way to increase the first-law efficiency in simple cycles, like steam Rankine, Organic Rankine and CO2 power cycles, although there are signs, that problems might arise when this increase reaches the vicinity of the critical point. Some of the problems related to the properties of the materials in this region; there are significant density fluctuations, causing strong local anomalies; also, the accuracy of the Equation of States can be reduced in this region [27,28]; this can also produce problems during design and in the operational phases. Also, from the technical point of view, high evaporation temperatures, especially for dry working fluids, are also not very favorable [29,30]. For example, choosing the evaporation temperature above the temperature of local entropy maximum (the so-called point M), expansion starts in the wet region, even in dry fluids; therefore, droplets can be formed during the early stage of the expansion, causing troubles in the special "dry" expanders.

Here, we approached the problem of the increase of maximal cycle temperature from a purely thermodynamic side, by determining its effect (with fixed condenser temperature) on the first-law efficiency. Efficiencies for six different materials with increasing dryness (water, CO2, R116, butane, neopentane and dodecane) were evaluated. Although anomalies were expected to be seen for dry working fluids [29,30], surprisingly, a very strange effect was seen for wet working fluids, namely the appearance of an efficiency maximum on the efficiency of the Rankine-like cycle (first seen for butane [26]), located close to the critical temperature. The location of the maximum depends on the condenser temperature as well as on the material. For water, taking the smallest possible condenser temperature (the triple point temperature, 273.16 K), the maximum appears 24.6 K below the critical point, while increasing the condenser temperature, it can approach the critical point. For CO2, R116 and butane, the biggest distance between the location of the maximum from the critical point is in the order of 5–10 K (8.68 K for butane, 4.01 K for CO2 and 4.53 K for R116); for neopentane, it is 1.97 K, and finally, for dodecane, no maximum can be seen. Concerning reduced temperatures, the biggest distances between the entropy maximum and the critical points are 0.038 (water), 0.013 (CO2), 0.015 (R116), 0.020 (butane) and 0.005 (neopentane). Since for most real ORC and RC power plants, the operational temperature is not very near to the critical temperature; therefore, technologically, the existence of the maxima are relevant only for two cases (for water, where the maximal distance between the critical temperature end efficiency maximum is almost 25 K, and for butane, where this difference is almost 9 K). For the rest of the presented working fluids, the relevance is rather theoretical than technical. However, one might assume that the temperature differences will be sufficient to have technological importance for some of the not presented wet working fluids.

The presence of these maxima can be clearly shown for several systems by comparing the efficiency of a Rankine-like cycle with a Trilateral Flash-like Cycle (TFC). Some theoretical explanation is also given to see, that for very dry working fluids (like dodecane), the disappearance of this maximum can always be expected for simple, Rankine-like cycles, using only the basic Rankine-layout (pump, liquid heater with evaporator, expander and condenser), without recuperative unit.

The existence of the efficiency maximum for various working fluids is a disproof of the common misconception that the increase of the maximal cycle temperature always has a positive effect on first law efficiency for simple thermodynamic cycles. This disproof is interesting itself, as a strange thermodynamic fact, but also, it can give a new tool for engineers to find optimal cycle temperature for an ORC (or similar) power plant, by designating a temperature range, where a further increase of the evaporation temperature would be contra-productive.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/1996-1 073/14/2/307/s1, Figure S1a,b: Rotating 3-D representation of the efficiency-curves for water and butane, respectively (in animated gif-format).

**Author Contributions:** Conceptualization, A.R.I.; analysis: A.M.A., preliminary analysis: L.K., writing: A.R.I., A.M.A. and L.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was performed in the frame of the FIEK\_16-1-2016-0007 project, implemented with the support provided from the National Research, Development and Innovation Fund of Hungary, financed under the FIEK\_16 funding scheme. Part of the research reported in this paper and carried out at BME has been supported by the NRDI Fund (TKP2020 NC, Grant No. BME-NC) based on the charter of bolster issued by the NRDI Office under the auspices of the Ministry for Innovation and Technology.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author, [A.R.I.], upon request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Hybrid Solar Systems**

Renewable energy is environmentally friendly; therefore, it is the most viable future option to dispose of some of the environmental problems. Solar is one of the best renewable energy sources, due to availability and high temperature, wherefore solar thermal power, or concentrated solar power (CSP), is an optimal technology to hybridize with other energy technologies for power generation. The hybrid systems using geothermal and solar energy is one of the options with good potentials. This hybrid technology has a promising future, especially in countries that have abundant solar and geothermal energy. The stand-alone geothermal energy can be classified as medium- or low-enthalpy heat source for power plants, therefore, adding solar energy is the best choice, due to increasing the temperature of a geofluid, which increases the first-law efficiency to convert heat to work (and then to electricity) [15,16,31].

In general, a hybrid system can be categorized based on the kind of used energy, such as merging the solar source with one of the renewable energy sources like geothermal, wind, and biomass (these are the so-called high-renewable hybrids); while merging it with natural gas is called medium-renewable hybrids. Finally, merging it with traditionally fueled Brayton or/and Rankine cycle, the systems are called as a low-renewable hybrid.

Reducing CO2 emission is a crucial factor for the power plant. The lowest CO2 emission is with high-renewable hybrids, and the rate increases with medium-renewable, and it is at its highest rate with the low-renewable hybrids; therefore, the increase of the solar part reduces the CO2 emissions [32]. The efficiency of hybrid systems can be increased by adding thermal energy storage that uses the surplus of solar thermal energy at the night time when the temperature is low. Furthermore, the increase of the solar collector field increases the efficiency in all case with or without using thermal energy storage [33]. The second-law efficiency for the hybrid geothermal-solar system is higher than for separate geothermal and solar system at all ambient temperatures [34]. Concerning the increase of maximal temperature above the critical point, the hybrid supercritical ORC (Organic Rankine cycle) outperforms the hybrid subcritical ORC thermodynamically and economically. Therefore, the hybrid supercritical ORC has higher thermal efficiency, and the ability to produce power two to nine times more than a stand-alone subcritical geothermal plant. The geothermal-solar plant is characterized by produce relative constant power during its period operation, concerned to stand-alone solar ones [18,35].

Concerning working fluids, some preliminary result was achieved, showing that that isobutene (also called isobutylene) is a suitable working fluid in the hybrid geothermalsolar system; it has the best performance with the lower GWP, when compared to other working fluids as isobutene, *n*-butane, and 2-butene [36]. The hybrid solar-geothermal is less vulnerable to the seasonal changes in ambient temperatures than stand-alone solar one, as it was observed in an Australian geothermal-solar hybrid system case study, demonstrating that a properly designed hybrid plant can outperform a stand-alone solar thermal plant in terms of the cost of electricity production [37].

As a final issue, the thermal stability of the working fluid and the capacity of the system components might limit the solar superheating in the hybrid geothermal-solar so, to overcome this problem, water can be used as the heat transfer fluid for ORC that operates under 175 ◦C [38].


**4.1. Adenocarcinoma**

12 Dehydrogenases

*3.1.4. ALDH2*

*3.1.5. ALDH7A1*

spleen, muscle, lung and brain.

It's been shown that ALDH1A3 takes part in the development of the eye, nucleus

ALDH1A3 deficiency has been shown to play a critical role in cancer by a number of studies. For instance, in human breast cancer MCF-7 cells, ALDH1A3 expression is downregulated, whereas in cultured human colon cancer cells, *ALDH1A3* is one of two genes that are upregulated by induction of wild type *p*53. In mammary tumor-susceptible BALB/cJ mice that are heterozygous for *p*53, *Aldh1a3* is one of five candidate genes located within a region determined for its linkage to mammary tumorigenesis. In mice resistant to induced mammary tumors, (C57BL/6J), *Aldh1a3* is one of the two upregulated genes. *ALDH1A3* is silenced by methylation in gastric cancer cells, whereas in glioblastoma

ALDH2 is a tetrameric enzyme expressed profusely in lungs and liver; it is also present in organs that obligate high mitochondrial capacity for oxidative ATP generation including heart and brain. Apart from that, ALDH2 is also important in the aldehydic substrate oxidation such as 4-HNE, acrolein, and short-chain, aromatic or polycyclic carbons. To add to its dehydrogenase activity, depending on the substrates, ALDH2 can function as an esterase and reductase. More recent attention has also been focused on ALDH2 in regards to its function in the biotransformation of nitroglycerin, reducing it to 1,2-glyceryl dinitrate for

ALDH7A1 is a homotetramer that's expressed in a large number of tissues; in rat heart, liver and kidney, increased levels of ALDH7A1 are noted, whereas in black seabream fish (sbALDH7A1), ALDH7A1 is significantly formed in the liver and the kidney, excluding the heart. In human fetal tissues, ALDH7A1 has been encountered at elevated levels in the cochlea, eye, ovary, heart and kidney. In contrast, balanced levels are detected in the liver,

Human ALDH7A1's primary role happens in the pipecolic acid pathway of lysine catabolism, in which it catalyzes the oxidation reaction of alpha-aminoadipic semialdehyde (AASA) (*K*m180 μM) to alpha-aminoadipate. *ALDH7A1* mutations form the molecular basis for pyridoxine-dependent epilepsy (PDE), an autosomal recessive disorder characterized by the aggression of tenacious seizures during infancy and early childhood and are avoidable

Remarkably, ADH7A1 expression in the cochlea of the ear, the region dependent on the healthy upkeep of internal hydrostatic pressure, clarifies that mammalian ALDH7A1 might have an accomplishable function in osmotic regulation and in hearing disorders. However, no connection has been revealed yet, including patients with the inner-ear disorder

accumbens and olfactory bulbs, the forebrain, hair follicles and the cerebral cortex.

cells, it is triggered by the antitumor agent IL-13 cytotoxin (Marchitti*, et al.*, 2008).

the production of nitric oxide, which is a critical vasodilator (Chen*, et al.*, 2010).

by daily use of high-dose pyridoxine (Vitamin B6) supplementation.

Ménière's disease, which effects hearing and balance.

Adenocarcinoma is an epithelium cancer that is generated from glandular tissue. Epithelial tissue includes, but is not limited to, the surface layer of skin, glands and a variety of other tissues that line the cavities and organs of the body. Epithelium can be derived from the three germ layers ectoderm, mesoderm and endoderm during embryologic period. Adenocarcinoma classification depends on not only being a part of the gland, but also depends on having the same secretory characteristics. But, this form of carcinoma can occur in some higher mammals, including humans (Fauquier*, et al.*, 2003).

Adenocarcinomas can arise in many tissues of the body due to the ubiquitous nature of glands within the body. While each gland may not be secreting the same substance, as long as there is an exocrine function to the cell, it is considered glandular and its malignant form is therefore named adenocarcinoma. Endocrine gland tumors, such as a VIPoma, an insulinoma, a pheochromocytoma, etc. are typically not referred to as adenocarcinomas, but rather, are often called neuroendocrine tumors. If the glandular tissue is abnormal, but benign, it is called an adenoma. Benign adenomas typically do not invade other tissues and rarely metastasize, whereas malignant adenocarcinomas do both. Colon, urogenital (cervical (Tewari*, et al.*, 2002), prostate, urachus and vagina), breast (Buchholz, 2009), esophagus, pancreas, stomach and throat are several examples of adenocarcinoma (Subramanian & Govindan, 2007).

It is reported that ALDH expression marks pancreatic cancer stem cells. Also, they have mentioned that the enhanced clonogenic growth and migratory properties of ALDHpositive pancreatic cancer cells suggest a key role in the development of metastatic disease that negatively affects the overall survival of patients with pancreatic adenocarcinoma (Rasheed*, et al.*, 2010).

### **4.2. Breast cancer**

From high-grade, absence of hormone receptor expression to positive HER2 status and the basal-like molecular subtype, the expression of ALDH1 is in direct relation with undesired tumor characteristics in breast cancer (Mieog*, et al.*, 2012).

Breast cancer cells with stem-cell-like properties are suggested to be responsible for metastatic spread. Aldehyde dehydrogenase 1 (ALDH1) and cluster of differentiation 44 (CD44) in addition to RhoC GTPase are among the stem cell markers that are expressed by these cells (Chaterjee & van Golen, 2011).

Breast CSCs were initially isolated, established on cell surface marker with CD24/lowCD44 expression. More currently, ''functional'' markers depending on stem cell properties are investigated for their plausible applications in the breast CSCs isolation. By this method, applying the aldefluor assay (Stemcell Technologies), originally designed to isolate viable HSCs and is an enzyme-based assay that recognizes ALDH activity, Ginestier et al. isolated breast CSCs. The assay is thought to precisely recognize ALDH isoform ALDH1A1 activity degree. Besides its application as a prognostic and CSC marker, ALDH activity that is primarily carried out by ALDH1A3 might be functional in breast cancer progression.

Aldehyde Dehydrogenase: Cancer and Stem Cells 15

ALDH1-positive, in contrast to CD44+/CD24-, was tremendously related to sequential paclitaxel- and epirubicin-based chemotherapy resistance, and the expression of ALDH1 increased after neoadjuvant chemotherapy, which stands for an indication of BCSCs, determined by ALDH1, indeed having played a significant role in chemotherapy resistance. This means that ALDH1-positive appears to be a better marker than CD44+/CD24- in identifying BCSCs, at least for the prediction of resistance to chemotherapy (Tanei*, et al.*,

Each year, approximately 171,000 new cases of lung cancer are diagnosed, and 160,000 individuals do not survive from the disease in the United States. This high incidence and mortality makes lung cancer one of the most common cancers and the leading cause of cancer death in men. Lung cancer is still the leading cause of death from malignant diseases worldwide in spite of the advances in surgical treatment and multimodality treatments

Cancer stem cells have attributed resistance of a smaller fraction of cells in the tumor bulk against chemotherapeutics. The isolation of CSCs is important for these reasons and have been isolated using a variety of stem cell markers and phenotypes. CD133 has recently been reported to identify tumor-initiating cells in non-small cell lung cancer (NSCLC). ABCG2 is also a stem cell marker of a variety of tissues and transporter responsible for the multidrugresistance phenotype. However, it was demonstrated that many cells in NSCLC and SCLC cell lines show tumorigenic potential, regardless of ABCG2 and CD133 expression. Recently, ALDH activity has been used for isolation of these kinds of cells. Normal SCs were shown to contain higher levels of ALDH activity than their more differentiated progeny. ALDHpositive cells of tumors have higher proliferation rates, migration and adhesion ability, and metastatic potential than ALDH-negative cells. This may occur because that RA product of ALDHs is thought to participate in cellular differentiation and stem cell self-protection

Epithelial ovarian cancer is the sixth most common cancer in women worldwide and it is still the most lethal gynecologic malignancy (Iorio*, et al.*, 2007)**.** Application of new technologies for detection of ovarian cancer could have an important effect on public health, but to achieve this goal, specific and sensitive molecular markers are essential (Petricoin*, et al.*, 2002)**.** Aldehyde dehydrogenase-1A1 (ALDH1A1) has been a valid marker among several malignant and non-malignant tissues in spite of several stem cell markers to identify CSCs. ALDH plays a role in the biology of TICs as well as being a stem cell marker. Because ALDH1A1 is implicated in chemo resistance pathways, it is questioned that targeting ALDH1A1 can effect cells resistant to chemotherapy and represent a potential target for cancer stem-cell-directed therapy. In a study, ALDH1A1 was investigated in ovarian cancer cell lines and patient samples and examined whether targeting ALDH1A1 sensitizes cells to

2009).

**4.3. Lung cancer**

(Hibi*, et al.*, 1998).

(Serrano*, et al.*, 2011).

**4.4. Ovarian cancer**

Expression of genes and tumor sphere formation in self-renewal and differentiation could be changed by adding chemical RA signaling inducers or inhibitors in breast cancer cell lines (Marcato*, et al.*, 2011).

ALDH1 could work as a marker of breast CSCs better than CD44+/CD24-. Though we could not maintain a conclusion that ALDH1 expression was significantly related with any conventional clinicopathologic attributes, nevertheless, there is a compelling relation between ALDH1-positive breast tumors and resistance to neoadjuvant chemotherapy, because of the pCR rates being obtained, which are lower in ALDH1-positive tumors (9.5%) than ALDH1-negative tumors (32.2%). Moreover, after neoadjuvant chemotherapy, a considerable increase in the proportion of ALDH1-positive tumor cells was observed. These results are an indication of ALDH1-positive tumor cells playing an important role in resistance to chemotherapy. Because of tumor cells being more tumorigenic than CD44+/CD24-, tumor cells of breast CSCs are thought to be richer in ALDH1-positive tumor cells than in CD44+/CD24- tumor cells. As a matter of course, we have shown that ALDH1 positive, in contrast with CD44+/CD24-, is closely associated with colony formation in the collagen gel as well. The subset of ALDH1-positive and CD44+/CD24- tumor cells has been reported to contain the largest proportion of breast cancer stem cells (BCSCs); consequently, it is speculated to have the strongest resistantance to chemotherapy. However, in our current study, pCR rates in the ALDH1-positive and CD44+/CD24- high subset (20%, 2 of 10), are not the lowest among all the subsets consisting of the ALDH1-positive and CD44+/CD24- low subset (0%, 0 of 11), the ALDH1-negative and CD44+/CD24- high subset (34.1%, 15 of 44), and the ALDH1-negative and CD44+/CD24- low subset (30.2%, 13 of 43). Adding CD44/CD24 status to ALDH1 status does not seem to positively improve the prediction of response to chemotherapy. Together, these results direct us to assume that, at least for the prediction of resistance to chemotherapy, ALDH1-positive tumor cells serve as a better marker for BCSCs than CD44+/CD24- tumor cells. Because such tumors contain a higher proportion of CSCs, we suppose that ALDH1- positive tumors are resistant to chemotherapy. However, because ALDH1 has been shown to play an important role in the resistance to chemotherapy in hematopoietic cells, ALDH1-positive tumor cells might be involved in resistance to chemotherapy, regardless of whether they are CSCs or not. In addition to deeper illumination of ALDH1's function in chemotherapy resistance in breast cancers, obtaining a significantly specific marker for BCSCs is necessary to enlighten an authentic role of BCSCs' chemotherapy resistance.

ALDH1-positive, in contrast to CD44+/CD24-, was tremendously related to sequential paclitaxel- and epirubicin-based chemotherapy resistance, and the expression of ALDH1 increased after neoadjuvant chemotherapy, which stands for an indication of BCSCs, determined by ALDH1, indeed having played a significant role in chemotherapy resistance. This means that ALDH1-positive appears to be a better marker than CD44+/CD24- in identifying BCSCs, at least for the prediction of resistance to chemotherapy (Tanei*, et al.*, 2009).


**6. Comparison with other recovery methods**

In this section we compare our version of OMP with an NLS optimization step for the sinusoid frequency and amplitude at each iteration to two common methods for CS recovery: OMP with a linear least squares amplitude estimator at each iteration and convex optimization based on the ell-1 norm of the sparse target vector plus the ell-2 norm of the measurement constraint given by eq. (2). It should be noted that most of the cases presented in the previous sections cannot be solved with OMP/LS or penalized ell-1 norm methods so it is necessary to pick a special case to even perform the comparison. Consider a noise-free signal that consists of 5 unity amplitude sinusoids at 5 different frequencies. We assume N=1024 time samples and an M=30 x N=1024 complex measurement matrix made up of the sum of random reals plus i times different random reals, both sets of reals uniformly distributed between -1 and 1.

#### **6.1 Baseline case OMP-NLS**

We performed 100 different calculations with the frequencies chosen by a pseudo-random number generator. In order to control the number of significant figures, we took the frequencies from rational numbers uniformly distributed between 0 and 1 in steps of 10-6. Table 2 shows the fraction of failed recoveries and the average standard deviation in the value of the recovered frequency as a function of the oversampling ratio.

#### **6.2 OMP with Linear Least Squares**

We performed the same 100 calcuations using conventional OMP in which the NLS step is replaced by LS as in Tropp and Gilbert (2007). Note that the number of failed recoveries is 184 Applications of Digital Signal Processing

Fig. 11. Standard deviation in frequency Vf (red-lower curve) and amplitude Va (green upper curve) for the case with input frequencies {0.3389, 0.3390}, unity amplitude and a 16x1024

In this section we compare our version of OMP with an NLS optimization step for the sinusoid frequency and amplitude at each iteration to two common methods for CS recovery: OMP with a linear least squares amplitude estimator at each iteration and convex optimization based on the ell-1 norm of the sparse target vector plus the ell-2 norm of the measurement constraint given by eq. (2). It should be noted that most of the cases presented in the previous sections cannot be solved with OMP/LS or penalized ell-1 norm methods so it is necessary to pick a special case to even perform the comparison. Consider a noise-free signal that consists of 5 unity amplitude sinusoids at 5 different frequencies. We assume N=1024 time samples and an M=30 x N=1024 complex measurement matrix made up of the sum of random reals plus i times different random reals, both sets of reals uniformly

We performed 100 different calculations with the frequencies chosen by a pseudo-random number generator. In order to control the number of significant figures, we took the frequencies from rational numbers uniformly distributed between 0 and 1 in steps of 10-6. Table 2 shows the fraction of failed recoveries and the average standard deviation in the

We performed the same 100 calcuations using conventional OMP in which the NLS step is replaced by LS as in Tropp and Gilbert (2007). Note that the number of failed recoveries is

value of the recovered frequency as a function of the oversampling ratio.

mixing matrix.

distributed between -1 and 1.

**6.1 Baseline case OMP-NLS**

**6.2 OMP with Linear Least Squares**

**6. Comparison with other recovery methods**

about the same as the baseline OMP-NLS but the frequency error is huge by comparison. This is the natural result of the frequency grid, which is the limit on the OMP resolution. Timing comparisons with our software show that OMP-NLS takes about 50% longer than conventional OMP. We have also windowed the OMP calculations in order to reduce "spectral leakage" and hopefully achieve better performance. Aside from the lowered failure fraction for *Nf* = 2, windowing OMP appears to have no statistically significant effect.

Table 2. Comparing OMP with NLS to OMP and OMP with windowing for 4 values of the overcomplete dictionary Nf = 1,2,4,8. (a) failure fraction, %. (b) rms error in recovered frequencies.

We have also compared windowed OMP to OMP/NLS in the presence of noise. Fig. 12 shows the frequency and amplitude errors, Vf and Va, as a function of the noise standard deviation V for OMP (blue) and OMP-NLS (red) for a signal composed of two sinusoids with *N* = 128, *M* = 20 and *N*f = 4 averaged over 100 trials with randomly chosen input frequencies. Note that the OMP frequency error drops to an asymptote of about 6 x 10-4 and the OMP amplitude error to about 0.23 for V < 0.1 while the OMP-NLS errors continue to drop linearly proportional to V for V < 0.1.


**2.2 Auger electrons**

Auger electrons are generated from suborbital transitions. They are typically very short-range emissions, of the order of 1–1000 nm, depending on their emission energy. Auger electron is intermediate (4–26 keV/μM) of LET [3, 22, 23]. These emissions could be highly cytotoxic if the RPT drug localizes within the cell nucleus [36, 37]. Bromine-77, indium-111, iodine-123, and iodine-125 are the most commonly used Auger electron emitters. In vitro studies had shown highly effective and specific tumor cell killing when they labeled with targeting vehicles that can localize these subcellular-range radiations close to cellular DNA [38, 39]. Human studies using locoregional administration showed promise regarding tumor cell incorporation of the Auger emitters [2]. However, Auger electron-emitter RPT has not been widely adopted yet. Besides, auger electron agents must be incorporated into the DNA, and their unfavorable pharmacokinetics might be the reasons for the lack of efficacy. Technological developments that could overcome those factors will interest RPT development [2].


**2. Theoretical tools**

parameters for all individuals in the population and random effects account for individual differences. Moreover, mixed effects modeling has been shown to be useful in pharmacokinetic/pharmacodynamic studies; particularly for modeling total variation in and between individuals, the models used in pharmacokinetic/ pharmacodynamic analysis are often presented by a system of deterministic differential equations (ordinary, partial, or delay). However, the pharmacological processes in reality are always exposed to incomprehensible effects that are difficult to model, and ignoring these effects can affect the parameter estimation results and their interpretations. So, the introduction of stochastic components to deterministic models is an important tool of analysis [1] and is more appropriate to model the intra-individual variations rather than ODEs. In addition, the extension of ODEs to the SDEs makes it possible to explain the differences between the observations and the predictions by two types of noise: dynamic noise that enters through the dynamics of the system and that can result from its random fluctuations or from the shortages of model and the measurement error which is added in the case of an indirectly observed process, which may be due to a test error or to the existence of a disturbance and represent the uncorrelated part of the residual variability. In the theory, there are rich and developed resources for mixed effects models whether deterministic [2–6] or stochastic and linear or nonlinear models (see many applications of stochastic NLME models in biomedical [7–9] and in pharmacokinetic

*Numerical Modeling and Computer Simulation*

Parameter estimation in mixed effects models with SDEs, known by stochastic differential mixed effects models, is not an obvious procedure except in some cases simpler [14] because it is often difficult to write the likelihood function in its closed form. In this context, we propose a review on estimation methods of SDME models in [13, 15] and an example case that treats a generalized linear mixed models in [16] and also an example to approximate the likelihood function of an NLME with the likelihood of a linear mixed effects model in [17]. Moreover, to strengthen knowledge on estimation methods of SDME, we refer to [18, 19] that propose an example of stochastic mixed effects model with random effects log-normally distributed

In general, it is difficult to obtain an explicit likelihood function because the transition density of the stochastic process is often unknown or that the integral in the likelihood given the random effects cannot be computed analytically, and although the size of the random effects increases, the complexity of the problem increases also rapidly. Therefore, this requires a significant need for approximate methods to compute the transition density in an approximate closed form and for effective numerical integration methods to compute or approximate the integral in the likelihood function, using, for example, the Laplacian and Gaussian quadrature approximation [3, 8, 20] or other approaches [21, 22]. In the literature, several solutions have been proposed to approximate the transition density and have shown their effectiveness despite certain limitations. For example, the transition density could be approximated by the solution of the partial differential equations of Kolmogorov [23]; by the derivation of a Hermite expansion of closed form at the transition density [24–26] (this method has been reviewed and applied for many known stochastic processes for one-dimensional [8] and multidimensional [20] time-homogeneous SDME model); or by simulating the process to Monte-Carlointegrate the transition density [27–29]. These techniques are very useful and can solve the problem, but unfortunately, they involve intense calculations which make

In this work, we focus on two fundamental issues concerning the implementation of SDEs in NLME models. The first is how the transition density of an SDME model can be approximated when it is not known, and the second is about

[10–12] fields).

with a constant diffusion term.

the problem always complicated.

**92**

Consider an N-dimensional continuous and stochastic process *Yt* in the state space *E* ⊂ *N* described by the general first-order nonlinear stochastic differential equations of the Itô type [30]:

$$d\mathbf{Y}\_t^i = \mu\left(\mathbf{Y}\_t^i, t, \theta, b^i\right)dt + \Sigma\left(\mathbf{Y}\_t^i, \theta, b^i\right)dW\_t^i, \quad \mathbf{Y}\_0^i = \mathbf{y}\_0^i, \quad i = 1, \ldots, M,\tag{1}$$

where *Yi t* is defined as the solution of the SDME model Eq. (1) that exists under some conditions that we suppose satisfied [31–33] and represents the observation of individual i from M different experimental units, (*i* ¼ 1, *::*, *M*), at the moment *t*≥*t i* 0, and *Yi* 0 ¼ *Yi t*0 is the initial state of *Yt* for each subject. The process *Yi t* � � *t*≥0, *i* ¼ 1*::M* n o is assumed to verify the same model structure Eq. (1) according to the individual parameters *bi* ; *θ* ∈ Θ ⊂ *p* is a p-dimensional fixed effects parameter which represents the same and common characteristics for all subjects, and *bi* ∈ *B*⊆ *q* are the q-dimensional individual random parameters assumed mutually independent, also called random effects because they are not the same for all individuals; they change between them according to a distribution of density *PB bi* jΨ � � depending on a population parameter Ψ; in the population approach, this parameter vector allows for a data from several subjects to be considered simultaneously. Each component *bi l* may follow a different distribution, ð Þ *l* ¼ 1, … , *q* , and a standard choice for the joint density function *PB bi* jΨ � � of the vector *bi* could be the Gaussian distribution; however, any other distributions may be considered continuous or discrete:

$$b^i \sim i.i.d.\ N(\theta, \phi). \tag{2}$$

The joint density function of the vector *bi* is parameterized by a q-dimensional parameter *ϑ*∈*υ*⊂ *q* and a *q* � *q*-dimensional matrix *ϕ*∈ Φ ⊂ *q*�*q* representing the covariance matrix of *bi* and specifying the parameters of the marginal distributions of the components *bi l* , 1ð Þ ≤*l* ≤*q* ; the components of *ϕ* and *ϑ* represent the population parameters Ψ. For *Yi* 0, it is not necessarily known, and when its components are unknown, they must be considered as random effects since they change between individuals; but in some cases it can be known and assumed equal to a real constant. Also, we assume that the distribution of *Yi* ð Þ*t* given *bi* , *θ* � � and *Yi t* 0 ð Þ¼ *yt*0, *t* 0 <*t*, has a strictly positive density with regard to the Lebesgue measure on E:

$$\forall y \to P\_Y(y, t - t'|\mathcal{Y}\_{t'}, b^i, \theta) > 0, \ y \in E. \tag{3}$$

*Wi* ð Þ*t* are the standard Brownian motions, and they are assumed mutually independent with *bj* for all 1≤*i*, *j*≤ *M*. The functions *μ*ð Þ� : *E* � � Θ � *B* ! and Σð Þ� : *E* � Θ � *B* ! þ represent, respectively, the drift and the diffusion term of the model and are assumed to have some properties sufficiently regular to ensure a unique solution to the model [33].

In this work, we propose an approximate method to estimate the transition density for a time-inhomogeneous NLME model with SDEs in a closed approximate form. So, we suggest to derive the transition density by solving the motion equations of the process using the Fokker-Planck equation, and we deal with the use of its solution given in [34]. Then, using the expression obtained, we get a closed form approximation of the likelihood function that we maximize to obtain the approximated MLEs ^*θ* and Ψ^. The approximated transition density obtained from the

*Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential…*

� � that we substitute in Eq. (4)

jΨ � �*dbi*

, the likelihood function Eq. (4)

2 *log* j�*H* ~

*b i* j*θ*, ΨÞj � i,

*b i* jΨ � � (8)

j*θ*, Ψ � �:

*iT* (9)

*b i* jΨ (7)

2 *log* ð Þ� 2*π* 1

j ~ *b i* , *θ* � � þ *logPB* ~

*:* (6)

*Y Yi j* , Δ*i j* j*Yi j*�1, *bi* , *θ*

> *P*ð Þ *a Y yi j* , Δ*i j* j*yi j*�1, *bi* , *θ*

For a multidimensional vector of random parameters, if the exact transition density or its closed form approximation is available, we can use the Laplace approximation method [16, 20, 35] to obtain an explicit expression of the likelihood function to maximize, despite the fact that the integral in Eq. (4) has no closed

> *b i* jΨÞ þ *q*

j*θ*, Ψ � � ¼ *logPY yi*

j ~ *b i* , *θ* � � þ *logPB* ~

**2.2 Approximate transition density for multidimensional and nonlinear SDME**

Usually, a formal general solution of the stochastic differential equations as in Eq. (1) cannot be given, which makes the calculation of the transition density for this process more complicated. Moreover, this process has a lot of fluctuations that its exact position cannot be determined but can be known given a region by its probability density; with the FP equation, such a probability density can be determined. The FP equation is a differential equation for the distribution function describing a Brownian motion by which the probability density of the stochastic process can be calculated in a much simpler way by solving this equation. This motion equation is usually used for variables describing a macroscopic but small

*∂*~ *b i ∂*~ *b*

h i � �

to obtain the following approximated likelihood function for the SDME model

� � !*PB bi*

proposed method is denoted by *P*ð Þ *a*

*DOI: http://dx.doi.org/10.5772/intechopen.90751*

*2.1.1 Laplace approximation*

can be approximated as:

*M*

*i*¼1

*logL*ð Þ *θ*, Ψ ≃ X

¼ *argmaxbi f bi*

where

**model**

~ *b i*

**95**

*L*ð Þ *a* ð Þ¼ *θ*, Ψ Y

*M*

ð Y*ni*

*j*¼1

*i*¼1

solution. So, for a q-dimensional random vector *bi*

j*θ*, Ψ � � � � *and f bi*

and ∣ � ∣ denotes the determinant of the Hessian matrix *H bi*

*∂*2 *logPY yi*

� �

*logPY yi* j ~ *b i* , *θ* � � þ *logPB* ~

*H* ~ *b i* j*θ*, Ψ � � ¼

Eq. (1):

According to the model Eq. (1), the process Y is the same and follows the same form for each individual of the population, and the SDME model of type Itô expresses the dynamic of the individual i perturbed by the Brownian motion. So, the differences between the subjects are modeled on the one hand by the different realizations of the Brownian motion paths *Wi t* � � *t*≥*t i* 0 and, on the other hand, by the

incorporation of the random parameters *bi* in the model. Therefore, the introduction of parameters varying randomly between subjects allows the model Eq. (1) to explain the variability between individuals.

The goal is to estimate the vector of fixed parameters *θ* and the parameter vector Ψ characterizing the distribution of random parameters *bi* s, but the statistical inference for such models is a difficult issue, and the level of difficulty is not the same whether the transition density is explicit or not and whether the process is observed directly or with measurement noise; in this work we assume that the process was exactly observed and no observation noise was considered.

#### **2.1 Maximum likelihood estimation in SDME model**

The likelihood function of an SDME model is expressed as follows:

$$L(\theta, \Psi) = \prod\_{i=1}^{M} P\left(\underline{\boldsymbol{\nu}}^{i} | \theta, \Psi\right) = \prod\_{i=1}^{M} \int P\_{\underline{\mathbf{\varnu}}}\left(\underline{\boldsymbol{\nu}}^{i} | b^{i}, \theta\right) P\_{B}\left(b^{i} | \Psi\right) db^{i} \tag{4}$$

with

$$P\_{\underline{Y}}\left(\underline{y}^{i}|b^{i},\theta\right) = \prod\_{j=1}^{n\_{i}} P\_{Y}\left(y\_{j}^{i},\Delta\_{j}^{i}|y\_{j-1}^{i},b^{i},\theta\right),\tag{5}$$

where *ni* is the number of observations for the subject i at discrete points of time *t i* 0, *t i* 1, … , *t i ni* n o, *i* ¼ 1, … , *M* and Δ*i j* ¼ *t i j* � *t i j*�1, *j* ¼ 1, … , *ni*. The conditional density *PY yi* j� � � is equal to the product of the transition densities Eq. (5) for given random effects and *θ*, but the availability of the explicit transition density is the second constraint for the statistical issue of model Eq. (1) to obtain an exact likelihood function and exact estimators, since computing the transition density is not always obvious and requires approximation methods. However, there are some cases where the exact likelihood function is known, and the exact MLEs of *θ* are obtained (see references in the introduction). In fact, to compute the likelihood function in a closed form of an SDME model, we can encounter two types of problems that require approximate methods to overcome them: First, when the transition density *PY yi j* , Δ*i j* j*yi j*�1, *bi* , *θ* � � is known but the integral in Eq. (4) has no solution, in this case, the numerical methods of approximation of the integral are required. Or, second, when *PY yi j* , Δ*i j* j*yi j*�1, *bi* , *θ* � � cannot even be expressed explicitly and must also be approximated, see next paragraphs. Usually, in realistic examples, we have both an unknown transition density and an integral that is difficult to solve analytically. In theory, several methods for approximating transition densities and integrals have been proposed (see references cited in the introduction).

*Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential… DOI: http://dx.doi.org/10.5772/intechopen.90751*

In this work, we propose an approximate method to estimate the transition density for a time-inhomogeneous NLME model with SDEs in a closed approximate form. So, we suggest to derive the transition density by solving the motion equations of the process using the Fokker-Planck equation, and we deal with the use of its solution given in [34]. Then, using the expression obtained, we get a closed form approximation of the likelihood function that we maximize to obtain the approximated MLEs ^*θ* and Ψ^. The approximated transition density obtained from the proposed method is denoted by *P*ð Þ *a Y Yi j* , Δ*i j* j*Yi j*�1, *bi* , *θ* � � that we substitute in Eq. (4) to obtain the following approximated likelihood function for the SDME model Eq. (1):

$$L^{(a)}(\theta,\Psi) = \prod\_{i=1}^{M} \left[ \left( \prod\_{j=1}^{n\_i} P\_Y^{(a)} \left( y\_j^i, \Delta\_j^i | y\_{j-1}^i, b^i, \theta \right) \right) P\_B \left( b^i | \Psi \right) db^i. \tag{6} \right]$$

## *2.1.1 Laplace approximation*

Σð Þ� : *E* � Θ � *B* ! þ represent, respectively, the drift and the diffusion term of the model and are assumed to have some properties sufficiently regular to ensure a

incorporation of the random parameters *bi* in the model. Therefore, the introduction of parameters varying randomly between subjects allows the model Eq. (1) to

inference for such models is a difficult issue, and the level of difficulty is not the same whether the transition density is explicit or not and whether the process is observed directly or with measurement noise; in this work we assume that the process was exactly observed and no observation noise was considered.

The likelihood function of an SDME model is expressed as follows:

¼ Y*ni j*¼1

*j* ¼ *t i j* � *t i*

random effects and *θ*, but the availability of the explicit transition density is the second constraint for the statistical issue of model Eq. (1) to obtain an exact likelihood function and exact estimators, since computing the transition density is not always obvious and requires approximation methods. However, there are some cases where the exact likelihood function is known, and the exact MLEs of *θ* are obtained (see references in the introduction). In fact, to compute the likelihood function in a closed form of an SDME model, we can encounter two types of problems that require approximate methods to overcome them: First, when the

solution, in this case, the numerical methods of approximation of the integral are

and must also be approximated, see next paragraphs. Usually, in realistic examples, we have both an unknown transition density and an integral that is difficult to solve analytically. In theory, several methods for approximating transition densities and

� �

¼ Y *M*

*i*¼1

*PY yi j* , Δ*i j* j*yi j*�1, *bi* , *θ*

where *ni* is the number of observations for the subject i at discrete points of time

is equal to the product of the transition densities Eq. (5) for given

� �

ð *PY yi* j*bi* , *θ* � �

The goal is to estimate the vector of fixed parameters *θ* and the parameter vector

form for each individual of the population, and the SDME model of type Itô expresses the dynamic of the individual i perturbed by the Brownian motion. So, the differences between the subjects are modeled on the one hand by the different

According to the model Eq. (1), the process Y is the same and follows the same

*t* � � *t*≥*t i* 0

and, on the other hand, by the

s, but the statistical

jΨ � �*dbi* (4)

, (5)

*PB bi*

*j*�1, *j* ¼ 1, … , *ni*. The conditional den-

is known but the integral in Eq. (4) has no

cannot even be expressed explicitly

unique solution to the model [33].

*Numerical Modeling and Computer Simulation*

realizations of the Brownian motion paths *Wi*

explain the variability between individuals.

*L*ð Þ¼ *θ*, Ψ Y

with

n o

j� � �

transition density *PY yi*

required. Or, second, when *PY yi*

sity *PY yi*

*t i* 0, *t i* 1, … , *t i ni*

**94**

Ψ characterizing the distribution of random parameters *bi*

**2.1 Maximum likelihood estimation in SDME model**

*M*

*P yi* j*θ*, Ψ � �

*i*¼1

*PY yi* j*bi* , *θ* � �

, *i* ¼ 1, … , *M* and Δ*i*

*j* , Δ*i j* j*yi j*�1, *bi* , *θ*

� �

*j* , Δ*i j* j*yi j*�1, *bi* , *θ*

integrals have been proposed (see references cited in the introduction).

For a multidimensional vector of random parameters, if the exact transition density or its closed form approximation is available, we can use the Laplace approximation method [16, 20, 35] to obtain an explicit expression of the likelihood function to maximize, despite the fact that the integral in Eq. (4) has no closed solution. So, for a q-dimensional random vector *bi* , the likelihood function Eq. (4) can be approximated as:

$$\log \mathcal{L}(\boldsymbol{\theta}, \boldsymbol{\Psi}) \simeq \sum\_{i=1}^{M} \left[ \log \boldsymbol{\mathcal{P}}\_{\underline{Y}} \left( \underline{\boldsymbol{y}}^{i} | \boldsymbol{\bar{b}}^{i}, \boldsymbol{\theta} \right) + \log \boldsymbol{\mathcal{P}}\_{\boldsymbol{\theta}} \left( \boldsymbol{\bar{b}}^{i} | \boldsymbol{\Psi} \right) + \frac{\boldsymbol{q}}{2} \log \left( 2\pi \right) - \frac{1}{2} \log \left| -H \left( \boldsymbol{\bar{b}}^{i} | \boldsymbol{\theta}, \boldsymbol{\Psi} \right) \right| \right], \tag{7}$$

where

$$\bar{b}^i = \arg\max\_{b^i} \left( f\left( b^i | \theta, \Psi \right) \right) \quad \text{and} \quad f\left( b^i | \theta, \Psi \right) = \log \underline{P}\_{\underline{Y}} \left( \underline{y}^i | \bar{b}^i, \theta \right) + \log P\_B \left( \bar{b}^i | \Psi \right) \tag{8}$$

and ∣ � ∣ denotes the determinant of the Hessian matrix *H bi* j*θ*, Ψ � �:

$$H\left(\tilde{b}^i|\theta,\Psi\right) = \frac{\partial^2 \left[\log P\_{\underline{Y}}\left(\underline{\boldsymbol{y}}^i|\bar{\boldsymbol{b}}^i,\theta\right) + \log P\_B\left(\bar{\boldsymbol{b}}^i|\Psi\right)\right]}{\partial \bar{\boldsymbol{b}}^i \partial \bar{\boldsymbol{b}}^{iT}}\tag{9}$$