Myopia has become one of the most significant challenges in global public health. One cross-sectional study suggested that the myopia rate among Chinese primary school children could be as high as 68.78%. Furthermore, it is projected that by 2050, up to 49.8% of the global population could have myopia, with high myopia affecting 9.8% of the population.1,2 Early engagement in near-work activities and insufficient outdoor time are major risk factors for the development of myopia.3 High myopia is associated with several irreversible blinding eye diseases.4 Consequently, preventing myopia and slowing its progression in children and adolescents is a major research focus.

Orthokeratology (OK) lenses have been proven to be safe, effective, and reversible interventions for myopia.5–7 These are corneal contact lenses featuring a reverse-geometry design. Overnight wear reshapes the cornea to create a multifocal surface, flattening the central curvature and steepening the mid-peripheral curvature. The treatment zone is the relatively flattened central area of the anterior corneal surface. It provides the best visual quality and minimal aberrations as light passes through. Previous studies have indicated that its area (Treatment Zone Area, TZA) is closely related to myopia control efficacy.8–10 Clinically, reducing the back optic (BC) zone diameter of OK lenses is a common strategy to shrink the TZA. However, in practice, the TZA alone cannot be precisely controlled. This suggests that other factors influence the TZA, which are not yet fully understood.

Currently, most studies rely on engineering software such as MATLAB for TZA quantification.11,12 While MATLAB is highly robust for matrix computations and topographic data extraction, transitioning from image processing in MATLAB to subsequent clinical statistical analysis often requires researchers to switch between different software environments. R language (R) is an open-source platform that has seen widespread adoption in biomedical and health sciences research.13 Although other open-source alternatives, such as Python, offer comparable programming capabilities, R was specifically selected for this study due to its extensive ecosystem of specialized biostatistical packages and its higher familiarity among clinical researchers and biostatisticians.13,14 Establishing a unified workflow in R that encompasses both topography data processing and downstream clinical modeling would substantially streamline the research process. However, the feasibility of using R to process raw corneal topography data and calculate the TZA has not yet been validated against established platforms like MATLAB.

Additionally, the mechanism by which OK lenses alter corneal shape remains unclear. The prevailing view is that corneal epithelial remodeling is the core process involved.15–17 Quantifying this epithelial remodeling at the cellular level is challenging. On the basis of the current literature, most studies calculate central epithelial thinning or mid-peripheral thickening to reflect the extent of remodeling. However, parameters quantifying the intrinsic relationship between these changes are rarely reported.18–20

Therefore, this study developed a TZA quantification method using R. This method aims to validate its agreement with traditional MATLAB analysis. This study also investigated the size and influencing factors of the TZA generated by different OK lens designs and different BC zones after one year of wear.

Orthokeratology is a reversible, nonsurgical myopia correction method. It can suppress axial elongation by 30% to 63%.23 Its efficacy and safety have been widely validated. OK lenses use their unique design to alter the corneal epithelial thickness profile. This reshapes the corneal curvature, creating a central treatment zone and a mid-peripheral defocus zone. This mechanism corrects myopia and slows its progression. Previous studies have demonstrated that a smaller TZA after OK lens wear enhances myopia control efficacy.9,10,24 Therefore, establishing an accurate and convenient TZA quantification method is essential for a deeper investigation of its regulatory mechanism.

Previously, measuring treatment zone parameters often involved importing differential corneal topography data into MATLAB for fitting and calculation. MATLAB remains highly advantageous for complex matrix computations and robust image processing. However, a notable drawback in the clinical research context is the fragmented workflow; researchers must often extract topographic data in MATLAB and then transition to separate statistical platforms for subsequent clinical modeling. Furthermore, recent bibliometric data demonstrate that R has achieved significantly wider adoption and familiarity than MATLAB within the health sciences and biomedical communities.13,14 While Python presents another viable open-source alternative for image processing, R offers a distinct advantage for clinical researchers due to its extensive ecosystem of specialized biostatistical packages. Our study confirms that for both TZD and TZA, the Cronbach's alpha and ICC values between the R and MATLAB calculations were greater than 0.950. Thus, rather than replacing MATLAB due to unreliability, R provides an integrated, highly accessible, and accurate analytical alternative. It enables researchers to perform both precise boundary detection and downstream statistical analysis within a single unified platform.

The morphological evaluation of the treatment zone in this study revealed that the treatment zone exhibited prominent elliptical geometric characteristics rather than an ideal circular shape, with a mean e of 0.46. Regarding the spatial distribution of lens decentration, our data were highly consistent with classic clinical observations, demonstrating that displacement was predominantly clustered within the temporal-inferior quadrant 53.97% and the temporal-superior quadrant 24.27%. This phenomenon may be attributed to Bell phenomenon during sleep and the intrinsic asymmetry of the cornea.25,26 Furthermore, chi square test results indicated that eyes with alpha between 0 and 90° primarily experienced temporal-superior decentration, whereas those with tilt angles between 90 and 180° were more susceptible to temporal-inferior decentration. Multiple linear regression analysis showed that the semimajor axis a, eccentricity e, and tilt angle alpha were independent factors influencing TZD. Notably, the vector of decentration aligned closer with the direction of the semiminor axis b of the ellipse, and a larger ellipse eccentricity e value was associated with a higher probability of decentration. Although the retrospective cross sectional nature of this study prevents the direct establishment of a causal relationship between elliptical morphology and decentration kinetics, we speculate that this localized mechanical compression may truncate the treatment zone boundary along the vector of lens displacement. The observed associations between ellipse morphology and TZD are directly supported by the present data, whereas their biomechanical interpretation remains hypothetical. One possible explanation is that asymmetric lens positioning and nonuniform tear forces may modify the treatment zone boundary along the direction of lens displacement. When an orthokeratology lens undergoes sustained directional translation during sleep, non coaxial hydrodynamic forces exert intense asymmetric shear stress on the corneal epithelium. This localized mechanical compression forcefully truncates the treatment zone boundary along the exact vector of lens displacement, physically shortening the geometric diameter in that specific meridian. Second, preexisting baseline corneal anatomical asymmetry induces an asymmetric distribution of hydrodynamic stress right from the initiation of lens wear. This uneven mechanical environment generates a lateral vector force, thereby guiding the lens to slide along the path of least resistance. However, the precise hydrodynamic and cellular level mechanisms remain hypothetical and warrant further experimental investigation.

Our results revealed significant differences in the TZA among the four groups. Further subgroup analysis revealed no significant differences in the TZA between the CRT-S and CRT-L groups, the CRT-S and VST-S groups, or the CRT-L and VST-S groups (adjusted P > 0.05). Although we observed that the overall TZA for the VST designs was greater than that for the CRT designs, the VST-S group had a BC zone of 5.5 mm, whereas the CRT-S group had a BC zone of 5.0 mm. This suggests that the differences in the TZA between different lens designs with small BC zones require further verification. However, for the VST-L and CRT-L groups, both with 6.0 mm BC zones, the TZA difference was significant (10.52 ± 2.31 vs. 8.16 ± 2.17, P < 0.001). This disparity may originate from the fundamental differences in lens geometry, where the VST features multiple curve segments joined together, whereas the CRT utilizes a mathematically continuous sigmoidal curve profile. These distinct geometric architectures may alter the volume of tear pooling within the reverse curve zone. According to Meng’s research, CRT lenses exhibit greater fluid dynamic tension within the return zone but exert less mechanical pressure on the central cornea beneath the base curve.27 Consequently, the post wear epithelial thickness within the mid periphery in the CRT group is identical to or greater than that in the VST group. Rather than expanding the treatment area, this highly concentrated fluid force and robust peripheral epithelial remodeling under the CRT design may act as a strict hydraulic barrier, effectively confining the treatment zone boundary and preventing outward expansion.27–29 This finding further suggests that, beyond the core design parameters of lens type and BC zone diameter, TZA formation and regulation were influenced by multiple interacting factors.

Previous studies on the factors influencing the TZA have often focused on static corneal morphological data. They largely assumed that preoperative corneal shape was the key predetermined factor for the TZA: Sun et al. reported that the TZA was significantly correlated with baseline SE and axial length;30 Li et al. reported that preoperative SE and flat keratometry were significantly correlated with the TZA;31 Gruhl et al. also reported a correlation between SE and treatment zone size;32 Ding et al. confirmed that the TZA was only correlated with corneal eccentricity and corneal height differences.33 Our results revealed that lens design, BC zone diameter, TZD, preoperative p, Δkm, ΔCCT, and Δp were significantly correlated with the TZA. Multiple linear regression analysis further confirmed that Δkm, TZD, Δp, preoperative p, and BC zone diameter were independent factors influencing the TZA. This model explained 73.0% of the variance in TZA size (R² = 0.730). However, because axial elongation and refractive progression were not modeled as outcomes, this R² value should not be interpreted as evidence of predictive performance for myopia control efficacy. Notably, Δkm had the most substantial effect on TZA (Standardized Coefficients = 0.963, P < 0.001). These findings confirm that, in addition to lens parameters and preoperative corneal morphology, the TZD and dynamic postoperative corneal remodeling process are highly correlated with TZA formation. Regarding the underlying mechanisms, we propose that Δkm represents the primary longitudinal expansion force, wherein a larger required myopia correction drives a greater volume of central epithelial cells to displace outwards, thereby expanding the treatment zone. However, the spatial redistribution profile of this displaced tissue is actively moderated by both preoperative p and Δp. Specifically, preoperative p serves as an anatomical predetermining factor, where a higher baseline value leads to a larger TZA. Crucially, when holding Δkm constant within the multiple regression framework, Δp emerges as a significant negative predictor of TZA size. We speculate that a larger Δp reflects a more abrupt, high density epithelial accumulation localized tightly within the paracentral return zone, which functions to confine the reshaping profile and thereby constricts the extent of the treatment zone.19,34

This finding reminds clinicians that for patients with a high preoperative p value, where significant changes in corrective curvature was expected, simply reducing the BC zone size might not achieve the ideal small TZA. Instead, further adjustments to lens fitting parameters, such as lens tightness or sagittal height, may be needed. These changes can alter the lens's corneal remodeling capability, controlling the magnitude of change in km and p, ultimately leading to a stable, smaller TZA. Our study did not find a correlation between the TZA and preoperative SE or km. This might be due to the multifactorial nature of myopia. Patients with the same SE can have different preoperative corneal curvatures and asphericities. Consequently, the reshaping ability of OK lenses also differs.

Simultaneously, the relationship between TZA and TZD should be interpreted within a broader optical and clinical context. In the present multivariable model, greater TZD was independently associated with a larger TZA, indicating a close relationship between treatment zone decentration and treatment zone morphology. Recent evidence suggests that orthokeratology generally shifts peripheral refraction toward myopic defocus, particularly at approximately 30 degrees of retinal eccentricity, and this optical change is considered one of the principal mechanisms potentially contributing to the inhibition of axial elongation.35 Previous studies have further shown that treatment zone decentration may alter the spatial distribution of peripheral refraction and higher order aberrations. In particular, temporal decentration may shift portions of the mid peripheral steepened region closer to the pupillary axis, thereby increasing myopic defocus within selected retinal regions.11,22,26,36 However, allowing unchecked decentration solely for the sake of axial inhibition is clinically unacceptable, as excessive lens displacement escalates the risks of corneal epithelial staining, localized mechanical injury, and degraded visual quality.37–39 Therefore, treatment zone decentration should be regarded as a potential clinical trade off between peripheral optical effects and the preservation of central visual quality, lens stability, patient comfort, and corneal integrity, rather than as a fitting objective that should be deliberately maximized. Although a certain degree of decentration may be clinically acceptable when satisfactory visual acuity, stable daytime vision, patient comfort, and corneal health are maintained, an optimal and safe range of decentration has not yet been established. Moreover, because peripheral refraction, higher order aberrations, visual quality, axial elongation, and myopia progression were not directly evaluated in the present study, the observed association between TZD and treatment zone morphology should be interpreted as a potential structural correlate of peripheral optical changes rather than direct evidence of improved myopia control efficacy. Future prospective studies integrating peripheral refraction, retinal image quality, corneal health, and axial length outcomes are required to determine whether specific patterns and magnitudes of decentration provide a favorable balance between potential myopia control effects and clinical safety.

Our study has several limitations. First, this was a retrospective study with a relatively limited sample size. Second, our investigation of epithelial remodeling was not based on direct epithelial measurements. This finding was inferred from changes in total corneal thickness. Third, we analyzed data from only the one-year follow-up. This cannot reflect the long-term dynamic changes in the TZA. Fourth, owing to clinical sample limitations, we did not include patients whose VST lenses had a 5.0 mm BC zone. Thus, we could not clarify TZA differences between lens designs at this specific BC diameter. The existing research findings on such differences were inconsistent.31,40 Finally, a notable limitation inherent in our retrospective design is that the allocation between VST and CRT lens designs was not randomly assigned. In routine clinical practice, practitioners typically select a specific lens design based on the baseline corneal characteristics of the patient. However, it is worth noting that all critical baseline demographic characteristics and ocular parameters were well balanced among the four groups with no statistically significant differences (Table 2). Furthermore, lens design was also excluded from the multivariate linear regression model. Future research should employ large sample, prospective cohort studies. The incorporation of technologies such as anterior segment OCT and ultrasound pachymetry for multiple time point measurements will help further explore the correlations between these parameters and the TZA.

This study confirms that R can be used to measure OK lens treatment zone parameters effectively, and the results were consistent with those of MATLAB. Treatment zone decentration tended to align with the semiminor axis of the fitted treatment zone ellipse, and greater ellipse eccentricity was associated with greater TZD. When the BC zone is 6.0 mm, VST-designed OK lenses produce a larger TZA than CRT-designed lenses do. Δkm, TZD, Δp, preoperative p, and BC zone diameter were independent factors influencing the TZA, and Δkm had the most substantial effect. Clinicians should carefully consider these factors during OK lens fitting. These findings improve the characterization of treatment zone morphology and may help clinicians anticipate morphological outcomes during individualized orthokeratology lens fitting

XL, YX, ZZ and LZ: design the study, data acquisition and analysis, manuscript writing; JZ: data acquisition and analysis; LZ: final approval of the manuscript.

This study was conducted in strict adherence to the ethical principles outlined in the World Medical Association Declaration of Helsinki (https://www.wma.net/policies-post/wma-declaration-of-helsinki/). Given the retrospective nature of the study, the Institutional Review Board (IRB) of The Third People's Hospital of Dalian University of Technology waived the requirement for written informed consent from participants, in compliance with the relevant national ethical regulations for retrospective medical research in China. All procedures involving human participants were performed in a manner that protected their privacy, confidentiality, and rights.

This study was supported by National Natural Science Foundation of China (82171032); Science and Technology Innovation Fund Project of Dalian (2023JJ12SN034); Major Dalian Peak Project (2022ky-01, 2022ky-02); Liaoning Province Science and Technology Plan Joint Program (2024JH2/102600030); the Fundamental Research Funds for the Central Universities, China (DUT25YG270).

Not Applicable. This study does not contain any identifying images, personal information, or clinical details that could compromise the anonymity of participants.