Quantitative magnetic resonance imaging and tumor prognosis

Abstract

An approach to data acquisition, analysis, and computational prediction using quantitative MRI data to predict cancer response to therapy is disclosed. An exemplary protocol details how to acquire the necessary images, followed by registration, segmentation, quantitative perfusion and diffusion analysis, model calibration, and prediction. An individual cancer patient's response to therapy is predicted by applying a biophysical reaction-diffusion model to these data. Applying the protocol, MRI data from at least two scanning visits are co-registered and the size, cellularity, and vascular characteristics of individual tumors are quantified. This allows for a spatially resolved prediction of how a particular patient's tumor will respond to therapy. Based on the predicted response, modified therapy can be determined.

Description

[Technical field]

[0001] CROSS-REFERENCE TO RELATED APPLICATIONS This application claims priority to and the benefit of U.S. Provisional Patent Application No. 63 / 257740, filed October 20, 2021, and U.S. Provisional Patent Application No. 63 / 247233, filed September 22, 2021, the entire contents of each of which are incorporated herein by reference.

[0002] STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH This invention was made with Government support under Grant Nos. U01CA142565 and U01CA174706 awarded by the National Institutes of Health. The Government has certain rights in this invention.

[0003] The present disclosure relates generally to predicting tumor response to therapy via quantitative magnetic resonance imaging (MRI) and biophysical reaction-diffusion modeling. [Background technology]

[0004] It is well known that neoadjuvant therapy (NAT) in standard care settings is not optimized for each cancer patient. Currently, therapy regimens are based on receptor status, tumor grade, body surface area, and genetic markers (each with known drawbacks) rather than spatially resolved physiological characteristics describing tumor features unique to an individual. Treatment plans may be changed in the event of lack of response, significant side effects, or consideration of the patient's quality of life, but this is done in a tentative manner. There is no mathematical theory in place to guide such decisions, forcing trial and error. Furthermore, certain existing clinical trial systems do not allow experimental evaluation of all possible combinations, timing, sequence, and administration strategies for each unique subset within a cancer population, let alone each individual patient. Thus, there is a need for an effective model that can also be used to determine individually optimized therapy regimens based on each patient's unique tumor characteristics, designed to predict tumor response in individual patients.

[0005] "Mechanism-based modeling" of cancer refers to incorporating biological mechanisms into models designed to predict the spatial and / or temporal dynamics of tumor characteristics. There is growing evidence that biophysical-mathematical models based on imaging information can accurately predict the development of kidney, prostate, brain, lung, pancreas, and breast cancers. These studies often aim to evaluate tumor growth or response to therapy on a patient-specific basis without first "training" the model on large-scale population data. That is, the model is calibrated with individual patient data, and then model-generated predictions are made of tumor response for each individual patient. Imaging data is a fundamental enabler of this process, as measurements can be collected in three dimensions (3D) at the time of diagnosis and at multiple time points during treatment, allowing patient-specific calibration and prediction. Not only does such an approach have the potential to predict individual patient responses, but this strategy may allow for establishing in silico twins for testing therapy regimens and optimizing treatment for each patient. Because the majority of oncology patients (85%) receive their care outside of research-oriented medical centers, parameterizing models with data accessible within an individual's specific local environment can have a dramatic and positive impact on patient care. Summary of the Invention

[0006] In various embodiments, disclosed herein are protocols for a complete data acquisition, analysis, and computational prediction pipeline for using quantitative magnetic resonance imaging (MRI) data to predict cancer response to therapy. The disclosed approach is applicable to heterogeneous patient populations. The protocols detail embodiments on how to acquire suitable images followed by registration, segmentation, quantitative perfusion and diffusion analysis, model calibration, and prediction. In certain embodiments, depending on the size of the tumor, the data collection portion of the protocol may require approximately 25 minutes of scanning, post-processing may require about 2-3 hours, and the model calibration and prediction components may require approximately 10 hours per patient. In various embodiments, an individual cancer patient's response to neoadjuvant therapy is predicted by applying a biophysical reaction-diffusion mathematical model to this data. In various embodiments, upon successful application of the protocol, MRI data from at least two scanning visits are co-registered and the size, cellularity, and vascular attributes of individual tumors are quantified. This allows for a spatially resolved prediction of how a particular patient's tumor will respond to therapy. In various embodiments, the disclosed protocols may employ expertise in image acquisition and analysis, as well as numerical solutions of partial differential equations, for example.

[0007] Various embodiments of the protocol involve obtaining quantitative MRI data of cancers, possibly at a community-based radiology center, analyzing this data to return spatially resolved maps of tumor physiology, and using the derived parameters to calibrate predictive, mechanism-based models of tumor growth and treatment response on a patient-specific basis. The goal of using medical imaging data to effectively inform patient-specific mathematical models of tumor growth and treatment response has been hindered by a lack of consensus on how the relevant data should be collected, processed, and modeled. The disclosed protocol embodiments also provide a (practical) learning tool that allows for the development of future strategies for using MRI data in mathematical tumor treatment.

[0008] The ability to accurately predict response and then rigorously optimize therapy regimens on a patient-specific basis would be transformative in tumor treatment. To this end, various embodiments provide an experimental mathematical framework that integrates quantitative MRI data into a biophysical model to predict patient-specific treatment response of locally advanced cancer to neoadjuvant therapy. In various embodiments, diffusion-weighted and dynamic contrast-enhanced MRI data are collected before therapy, after one cycle of therapy, and at the completion of the first therapy regimen. In one study, the first two patient-specific MRI data sets are used to initialize and calibrate the model to predict response at the third treatment and compare this to patient outcomes (N=18). Model predictions for total cellularity, total volume, and longest axis at the completion of the regimen are significant (p<0.05) within the expected measurement accuracy and strongly correlate with the measured response (p<0.01). Furthermore, the model is used to in silico explore various (practical) alternative treatment plans to achieve the greatest possible tumor control for each individual in a subgroup of patients (N=13). The model identifies alternative dosing strategies predicted to achieve better tumor control compared to standard care for 12 of 13 patients (p<0.01). In summary, the mechanism-based predictive model demonstrated the ability to identify alternative treatment regimens predicted to be clinically superior to a patient's therapy regimen, with important implications for clinical trial design and opportunities to change oncology care in the future.

[0009] Various embodiments relate to a method including acquiring, by a computing system, magnetic resonance imaging (MRI) data corresponding to a plurality of MRI scan images of an anatomical region including a tumor of a patient, the plurality of MRI scan images including a first image set obtained through a first scan performed prior to administration of a therapy to the patient and a second image set obtained through a second scan performed after administration of the therapy to the patient, the therapy including administration of a plurality of drugs; determining, by the computing system, tissue characteristics of tissue surrounding the tumor from the MRI data; registering, by the computing system, image-related data generated from the second image set to the image-related data generated from the first image set; determining, by the computing system, diffusion characteristics of the tumor based on the tissue characteristics; determining, by the computing system, growth characteristics of the tumor based on the tissue characteristics; determining, by the computing system, for each drug of the plurality of drugs, an effect of the drug on tumor cells; and generating, by the computing system, a score indicative of a predicted response of the tumor to the therapy based on the diffusion characteristics of the tumor, the growth characteristics of the tumor, and the determined effect of each drug on cells included in each voxel.

[0010] In various embodiments, the method includes performing tumor segmentation to identify a tumor region of interest (ROI) based on the MRI data prior to determining the tissue feature. In various embodiments, the tissue feature is related to vasculature within the tissue. In various embodiments, the MRI data includes dynamic contrast-enhanced MRI (DCE-MRI) data. In various embodiments, the tissue feature is quantified based on the DCE-MRI data. In various embodiments, the tissue feature is quantified based on a pharmacokinetic model and / or based on a fluid dynamic model. In various embodiments, the tissue feature is quantified based on a Kety-Tofts model and / or a variant of the Kety-Tofts model. In various embodiments, the MRI data includes diffusion weighted MRI (DW-MRI). In various embodiments, the method includes generating a map of apparent diffusion coefficient (ADC) of water. In various embodiments, registering the image-related data from the second image set to the image-related data generated from the first image set aligns the image to the map within a common domain. In various embodiments, the method includes determining drug distribution within each voxel of the tissue. In various embodiments, determining the drug distribution within each voxel of the tissue includes generating a normalized map of blood volume to define an initial drug distribution across the domain at the time of each administration of therapy. In various embodiments, the diffusion properties correspond to the diffusion of tumor cells mechanically associated with the material properties of the tissue. The diffusion of tumor cells may be mechanically associated with the material properties of the tissue via a physical stressor such as von Mises stress. The diffusion properties may represent changes in the tumor that may cause tissue deformation. In various embodiments, the growth properties are based on a carrying capacity related to the maximum number of tumor cells that can physically fit within a voxel. In various embodiments, the growth properties are based on a growth rate per voxel. In various embodiments, the growth rate is calibrated for each voxel within the patient's tumor ROI. In various embodiments, the effect of each drug on the tumor cells is based on the concentration of the drug in the tissue. In various embodiments, the effect of each drug on the tumor cells corresponds to the spatiotemporal distribution of each drug in the tissue.In various embodiments, the effect of each drug on the tumor cells is based on one or more of the drug's efficacy parameter α, the drug's washout parameter β over time after each administration, and / or the drug's initial concentration. In various embodiments, the efficacy parameter and / or washout parameter are calibrated to the patient and / or the drug. In various embodiments, the calibration of the patient's washout parameter is constrained using a boundary defined from the range of the drug's terminal elimination half-life. In various embodiments, the MRI data includes dynamic contrast-enhanced MRI (DCE-MRI) data. In various embodiments, the initial concentration is approximated using the DCE-MRI data. In various embodiments, the method includes performing rigid and / or non-rigid intra-scan image registration of the MRI data in each of the first image set and the second image set. In various embodiments, the method includes determining a modified therapy based on a score indicative of the tumor's expected response to the therapy. In various embodiments, the method includes administering the modified therapy based on a score indicative of the tumor's expected response to the therapy. In various embodiments, the therapy is a neoadjuvant therapy (NAT). Alternatively or additionally, in various embodiments, the correction therapy is NAT. In various embodiments, the therapy is a first NAT (such as chemotherapy, radiation therapy, or hormone therapy) and the correction therapy is a second NAT (a different type of therapy than the type of therapy of the first NAT, or a different therapy of the same type as the first NAT). In various embodiments, both the therapy and the correction therapy can be before the surgical intervention, or both the therapy and the correction therapy can be after the surgical intervention. In various embodiments, the therapy can be before the surgical intervention and the correction therapy can be after the surgical intervention.

[0011] Various embodiments relate to a computing system comprising one or more processors and a computer readable memory having instructions configured to cause the one or more processors to: acquire MRI data corresponding to a plurality of MRI scan images of an anatomical region including a tumor, the plurality of MRI scan images including a first image set obtained through a first scan performed prior to administration of a therapy to the patient and a second image set obtained through a second scan performed after administration of the therapy to the patient, the therapy including administration of a plurality of drugs; determine from the MRI data one or more tissue attributes of tissue surrounding the tumor; register image-related data generated from the second image set to the image-related data generated from the first image set; determine diffusion characteristics of the tumor based on the tissue attributes; determine growth characteristics of the tumor based on the tissue attributes; determine, for each drug of the plurality of drugs, an effect of the drug on tumor cells; and generate a score indicative of a predicted response of the tumor to the therapy based on the diffusion characteristics of the tumor, the growth characteristics of the tumor, and the determined effect of each drug on cells included in each voxel.

[0012] In various embodiments, the instructions are configured to cause the one or more processors to perform tumor segmentation to identify a tumor ROI based on the MRI data prior to determining the tissue feature. In various embodiments, the tissue feature is related to vasculature within the tissue. In various embodiments, the MRI data includes dynamic contrast-enhanced MRI (DCE-MRI) data. In various embodiments, the tissue feature is quantified based on the DCE-MRI data. In various embodiments, the tissue feature is quantified based on a pharmacokinetic model and / or a fluid dynamic model. In various embodiments, the tissue feature is quantified based on a Kety-Tofts model and / or a variant of the Kety-Tofts model. In various embodiments, the MRI data includes diffusion weighted MRI (DW-MRI). In various embodiments, the instructions are configured to cause the one or more processors to generate a map of apparent diffusion coefficient (ADC) of water. In various embodiments, registering the image-related data from the second set of images to the image-related data generated from the first set of images aligns the images to the map within a common domain. In various embodiments, the instructions are configured to cause the one or more processors to estimate the drug distribution within each voxel of the tissue. In various embodiments, determining the drug distribution within each voxel of the tissue includes generating a normalized map of blood volume to define an initial drug distribution across the domain at the time of each administration of the therapy. In various embodiments, the diffusion characteristics correspond to the diffusion of tumor cells mechanically associated with the material properties of the tissue via a physical stressor (e.g., force per area), thereby representing changes in the tumor that may cause deformation in the tissue. In various embodiments, the growth characteristics are based on a carrying capacity related to the maximum number of tumor cells that can physically fit within a voxel. In various embodiments, the growth characteristics are based on a growth rate per voxel. In various embodiments, the growth rate is calibrated for each voxel within the patient's tumor ROI. In various embodiments, the effect of each drug on the tumor cells is based on the concentration of the drug in the tissue. In various embodiments, the effect of each drug on the tumor cells corresponds to the spatiotemporal distribution of each drug in the tissue.In various embodiments, the effect of each drug on the tumor cells is based on one, more than one, or all of the following: an efficacy parameter α of the drug, a washout parameter β of the drug over time after each administration, and / or an initial concentration of the drug. In various embodiments, the efficacy parameter and / or the washout parameter are calibrated to the patient and / or the drug. In various embodiments, the calibration of the patient's washout parameter is constrained using a boundary defined from a range of terminal elimination half-lives of the drug. In various embodiments, the MRI data includes DCE-MRI data. In various embodiments, the initial concentration is approximated using the DCE-MRI data. In various embodiments, the instructions are configured to cause the one or more processors to perform rigid and / or non-rigid intra-scan image registration of the MRI data in the first image set and / or the second image set. In various embodiments, the instructions are configured to determine a modified therapy based on a score indicative of a predicted response of the tumor to the therapy. In various embodiments, one or both of the therapy and / or the modified therapy is / are a neoadjuvant therapy (NAT).

[0013] The foregoing summary is illustrative and is not intended to be in any way limiting. In addition to the exemplary aspects, embodiments, and features described above, further aspects, embodiments, and features will become apparent by reference to the following drawings and detailed description. [Brief description of the drawings]

[0014] [Figure 1A] 1 is an exemplary system for implementing the disclosed tumor prediction approach, according to various potential embodiments. [Figure 1B] 1 is an exemplary process for tumor prediction, according to various potential embodiments. [Figure 1C]Overview of an exemplary protocol, according to various potential embodiments. Each panel includes summary keywords for each section. The procedure has five main sections: patient population definition (step 1, not shown), image acquisition (steps 2-9), data analysis (steps 10-25), mapping of imaging data to a mathematical model (steps 26-36), and tumor prediction (steps 37-40). Each section of the procedure focuses on a specific area of ​​the protocol, and each section can be adapted to alternative investigations or used independently given specific circumstances of other studies. For example, the image acquisition section can be adapted to MRI studies of other organs. Also, the mapping and prediction sections can be applied given imaging data that has already been acquired and analyzed. DW-MRI: Diffusion-weighted magnetic resonance imaging, DCE-MRI: Dynamic contrast-enhanced MRI. [Diagram 2] Timeline of MRI acquisition with an exemplary standard of care NAT regimen for triple-negative breast cancer consisting of two therapy regimens, according to various potential embodiments. Panel (a) depicts the exemplary NAT regimen alone, and panel (b) depicts the exemplary regimen with protocol calibration and prediction strategies. In both panels, the thick arrows indicate the first administration and the start of each cycle (i.e., administration of a single drug or combination of drugs over a specified period of time, usually 2-4 weeks), and the thin arrows indicate any additional administrations within each cycle. In the first regimen, the dotted arrow represents a combination of doxorubicin and cyclophosphamide (usually consisting of 4 cycles, with the drugs administered once at 2-week intervals). In the second regimen, the solid arrow represents paclitaxel (usually consisting of 4 cycles, with the therapy administered weekly, with each cycle lasting 3 weeks. The thin solid arrow represents an additional administration). Some patients are treated with a combination of carboplatin and paclitaxel, with carboplatin administered only during the first week of each paclitaxel cycle (solid arrows). After NAT is completed, patients undergo surgery as part of the standard of care to determine pathological response. The protocol has MRI data collected before and after the first cycle of each therapy regimen. [Diagram 3] Fuzzy c-means (FCM) clustering to generate tumor regions of interest (ROIs) according to various potential embodiments. Sagittal breast sections of average DCE-MRI data from one patient (all panels) are depicted. The center panel shows a manually drawn ROI identifying the bounding polygon of the unobtrusive tumor. The right panel depicts the resulting ROI generated from the FCM algorithm within the manually drawn bounding polygon. [Figure 4] Comparison of inter-visit registration results with and without incorporating tumor ROI penalty into the registration scheme according to various potential embodiments. Panel (a) is the target image (scan image 2 defined in FIG. 2) and panel (b) is the "moving" image that is deformed / shifted to match the target image (scan image 3). In panel (b), the tumor ROI in the moving image is shown with a black outline. Panel (c) represents the result of registering the moving image to the target image using the approach described in steps 27-31 of the main text (i.e., rigid + non-rigid B-spline with tumor ROI penalty). To show the resulting deformation due to the registration process, panel (d) depicts a representative grid of the original moving image, while panel (e) represents the resulting deformed grid after registration (corresponding to the registered image in panel (c)). Panels (f) and (g) depict the deformation of a representative grid after rigid registration only (translation and rotation only) and non-rigid B-spline registration without tumor ROI penalty, respectively. White circles have been added across panels (b-g) to aid in comparison of the fields of area surrounding the tumor. Note that the inclusion of tumor ROI penalty results in less deformation of the tumor ROI; that is, there is less deformation within the white circle in panel (e) compared to panel (g). This procedure was applied to the MRI datasets of scans 1, 3, and 4. [Diagram 5]FIG. 1 is a flow chart of data analysis steps of an exemplary protocol, according to various potential embodiments. Both step numbers and step names are provided. Early data analysis steps (10-19) are dependent on the completion of previous steps, while some of the later steps (20-24) can be performed in parallel. [Figure 6] Exemplary image acquisition results according to various potential embodiments. A central slice of an exemplary patient depicting 200s / mm2b values ​​from DW-MRI (a), flip angle (FA) ratio from B1 map (b), a 10° T1 weighted acquisition (c), and the mean signal intensity of the DCE-MRI data across all kinetics (d). Tumor volume is indicated by a box. [Figure 7] Exemplary results from data analysis according to various potential embodiments. Each of the MRI data is aligned using a rigid body registration algorithm, interpolated to the same resolution, and aligned across visits. The DW-MRI data is used to calculate the ADC map (panels a and e). The DCE-MRI data is used to identify the tumor ROI (panels b and f). Multi-flip angle (MFA) T1 scan images and B1 map correction are used to calculate the T1 map (panels c, d, and h, respectively). The DCE-MRI data (panel b) and the T1 map (panel h) are used to calculate the Kety-Tofts model parameters (panel g) within the tumor ROI (panel f). [Figure 8] Conversion of imaging data into physical quantities for mathematical modeling, according to various potential embodiments. Before deriving modeling quantities, inter-visit registration is required to align images across all visits (panel a, details provided in FIG. 4). Once aligned, the ADC map (panel b) is used to calculate tumor cellularity (panel c). The DCE-MRI data (panel d) is used to identify fibroglandular and adipose tissue (panel e) using a fuzzy k-means algorithm. The patient-specific Kety-Tofts model parameters are used together with each patient's individual therapy regimen (panel f) to derive the approximate drug distribution within the tumor tissue (panel g). [Figure 9]Comparison of 3D model predicted results (across three central slices, left column) with observed results at the third scan time (right column) for one exemplary patient, according to various potential embodiments. The number of tumor cells is shown overlaid in color on each anatomical image. Although the model captures the correct shape of the tumor (Dice coefficient = 0.79), and areas of higher and lower cellularity may not match directly, there is a similar scale of cell density between the model prediction and the patient's actual tumor, resulting in Pearson and concordance correlation coefficients of 0.80 and 0.78, respectively. Comparison of summary measures of predicted and measured tumors in scan image 3 yields differences of 13, 17, and 4% in total cellularity, volume, and longest axis, respectively. [Figure 10] Graphical depictions showing the integration of breast MRI data with a mathematical model to predict tumor response and simulate alternative treatment regimens, according to various potential embodiments. Data from MRI scans obtained before and after one cycle of the initial NAT regimen were used to generate spatially resolved maps of tumor cell count and drug delivery (panels a and b, respectively). This initial NAT data was then used to calibrate the model parameters (panel c), and the model was run up to the time of the patient's third scan (panel d). Model predictions for total cellularity, total volume, and longest axis measures were then directly compared to each patient's actual tumor outcome as determined by scan 3 data. Model predictions were also compared to RECIST specifications to determine accuracy in terms of clinical measures. Using the patient-specific model parameters (panel c), alternative treatment regimens that adjust the frequency and dosage of treatment ("Tx1", panel e) were evaluated using model predictions from scan 2 to scan 3 for each of these regimens to determine the optimal treatment schedule for each patient (panel f). [Figure 11]Exemplary expected response results to NAT regimens for one patient (Patient 4) who received a standard of care regimen consisting of a combination of doxorubicin and cyclophosphamide every 3 weeks for 12 weeks (4 cycles total) according to various potential embodiments. The figure depicts total tumor cellularity overlaid in color on an anatomical image of a central slice of the breast. Although there is not a perfect match between the patient's actual scan 3 data (panel a) and the expected values ​​for the standard of care regimen (i.e., the standard of care regimen the patient actually received, panel b), the percentage difference between the expected and measured tumor response to the standard of care regimen is 1%, 16%, and 1% for total cellularity, volume, and longest axis. Also depicted are expected total cellularity for two alternative regimens: 1 / 3 dose administered weekly (panel c) and 1 / 21 dose administered daily (panel d). (Note: for each alternative regimen, the total drug amount during the treatment period was the same as the standard of care regimen). Across all regimens, the model predicted the greatest tumor cytopenia when a daily dosing regimen was implemented. Compared to the standard of care regimen, the model predicted that the daily dosing regimen would reduce total tumor cellularity by an additional 45%. [Figure 12]Distributions of sample means between expected and measured outcomes (i.e., total cellularity, volume, and longest axis) generated by Monte Carlo resampling (N=500 for each distribution) across patient cohorts according to various potential embodiments. (For details on the construction of these plots, see Supplementary Material below.) Each panel depicts the distribution of randomly sampled differences between each expected and measured outcome. The red vertical lines indicate the mean difference for the cohort. (Note that the red lines are not in the middle of the distributions because they do not represent the mean of the sampled distributions.) For example, panel a represents the sample mean difference between expected and measured total cellularity assuming an absolute difference between expected and measured total cellularity of 10%. Panels b and c display similar data for volume and longest dimension, respectively. Panels d-f correspond to panels a-c, but assume an absolute difference of 15% between expected and measured total cellularity, volume, and longest axis, respectively. Panels g-i correspond to panels a-c, but assume an absolute difference of 20% between the expected and measured total cellularity, volume, and longest axis, respectively. [Figure 13]Scatter plots comparing actual measurements from scan 3 for each patient (N=18) with model predicted values ​​according to various potential embodiments, and corresponding correlation coefficients (CCC, PCC, and KCC) by group. The dashed lines indicate the 45° line of agreement. Note that there is a greater correlation in the chemotherapy subgroup compared to chemotherapy+ (i.e., patients who received targeted therapy or immunotherapy in addition to chemotherapy) for all three tumor response measures. Panel (a) depicts the predicted total cell count compared to the measured total cell count estimated from DW-MRI data. Panel (b) depicts the predicted total tumor volume compared to the measured tumor volume, as determined by the total number of voxels within the tumor ROI. Panel (c) depicts the predicted longest axis of the tumor compared to the actual measured longest axis, as determined from the tumor ROI of scan 3. Note that using the logarithmic scale in panels (a) and (b), the patient with zero measured tumor at scan 3 (chemotherapy subgroup, patient 4) is not shown, although the corresponding CC value includes this data. [Figure 14] Comparison of percent change in total tumor cellularity from scan 2 to scan 3 achieved by standard regimen, 2x frequency regimen at 1 / 2 dose, 3x frequency regimen at 1 / 3 dose, 4x frequency regimen at 1 / 4 dose, and equivalent daily dose when administered daily to each patient according to various potential embodiments (panel a). Across patients, the difference in percent change between the standard therapy regimen and the most effective regimen identified in the model was a maximum difference of 46%, a minimum difference of 0%, and a median difference of 17% (panel b). Note: Positive differences indicate a potential additional percent reduction in total tumor cellularity achieved by the alternative regimen compared to the standard dose received by the patient. The standard regimen group was found to be statistically inferior to the group of the optimal regimen selected for each patient in terms of tumor control / reduction (p<0.001 for percent change from scan 2 to 3 and expected cell count at scan 3). Most effective regimen by patient population: standard N=1, halved N=4, third N=2, quarter N=2, daily N=4. [Figure 15] The MRI data undergoes four distinct steps before being incorporated into the mathematical model in various potential embodiments: intra-scan registration, data processing, inter-scan registration, and calculation of modeling quantities. Details are provided in the text and figures below, but briefly, inter-scan registration involves aligning different MRI data types (panels a and b) within each scanning session to compensate for any possible motion that occurs during the scanning visit. The data processing step involves deriving values ​​from each of the MRI data sets to identify the region of interest (panel c), quantify the vascular structure (panel d), and define values ​​that determine the tumor cell density (panel e). Inter-scan registration involves aligning all MRI scans of each patient over time (panel f). The final step involves taking the aligned data and defining specific quantities to be used in the mathematical model, including breast area for modeling (panel g), tissue map (panel h), tumor cell count (panel i), and drug concentration (panel j). [Figure 16] FIG. 1 is a simplified block diagram of representative server and client computer systems that can be used to implement certain embodiments of the present disclosure. [Figure 17]Two frameworks for predicting patient-specific response to NAST, according to various potential embodiments. Panel A shows the timeline of treatment administration and data acquisition for each patient. Panel B shows the processing modeling pipeline for generating patient-specific digital twins. Two frameworks are established to evaluate the predictive capabilities of digital twins. Framework 1 (Panel C) uses digital twins to predict the outcome of doxorubicin and cyclophosphamide (A / C) regimens. Patient-specific images from Visit 1 (V1) and Visit 2 (V2) and the A / C schedule provide the inputs on which the digital twin is calibrated. Once calibrated, the digital twin outputs a prediction of the tumor spatiotemporal evolution according to A / C. The prediction is then directly compared with the V3 images. Framework 2 (Panel D) uses digital twins to predict the overall NAST outcome. Images from V1, V2, and V3 and both A / C and paclitaxel schedules are given as inputs, and the digital twin outputs a prediction of whether the patient will achieve pCR. The prognosis is then directly compared with the pathologic response after surgery. [Figure 18] FIG. 1 shows a flow chart of MRI data processing, according to various potential embodiments. Panel A shows an example of DW-MRI and DCE-MRI data acquired from a single patient visit. Panels B-D show three steps of the processing pipeline, respectively. In panel B, the multi-parametric images are cropped and aligned to the same field of view (FOV). In panel C, the V1 and V3 images are aligned to the V2 image. In panel D, tissue segmentation and tumor cellularity calculation (from the DW-MRI data) are performed. These steps prepare the data for calibration with a biology-based mathematical model and establishing a digital twin for each patient. [Figure 19]Temporal accuracy of patient-specific prediction of TNBC response to A / C according to various potential embodiments. Panels A and B show the time course of calibrated therapy efficacy for A / C regimen in two representative patients, respectively. Panels C and D show the temporal dynamics predicted by the digital twin of the same two patients, with subpanels (i) and (ii) showing the changes in tumor cellularity (TTC) and tumor volume (TTV), respectively. In each panel, the red circle represents the measured TTC or TTV at a particular time point, and the blue curve and shading represent the median and range of predicted dynamics, respectively. Very small differences are observed between the measured and predicted changes in TTC and TTV over time in the exemplary patient. Panels E and F plot the correlation (CCC=0.95 and 0.94) between the measured and predicted changes in TTC and TTV in the cohort, respectively. These results show high accuracy and precision for predicting the temporal dynamics of patient-specific TNBC in response to A / C. [Figure 20] Spatial accuracy of patient-specific prediction of TNBC response to A / C according to various potential embodiments. Panels A and B show measured and predicted tumor cell distributions on the median tumor slices of two patients. Panels C and D show 3D renderings of the measured and predicted changes in these two tumor shapes. Very small differences are observed between the measured and predicted tumor cell distributions or tumor shapes of the patients. Panel E represents the difference between the measured and predicted changes in tumor cell distributions in the cohort. The median (red circles) and interquartile range (blue bars) of the difference within the tumor area of ​​each patient are represented. The mean (95% Cl) of the difference across all patients was 0.20% (-20.35% to 20.75%). These results demonstrate the high accuracy of the digital twin to predict the spatial dynamics of patient-specific TNBC in response to A / C. [Figure 21]Accuracy of patient-specific prediction of final pathological response according to various potential embodiments. Panels A and B show the time course of calibrated therapy effect during NAST in two exemplary patients, respectively. Panels C and D show the temporal dynamics predicted by the digital twin of the same two exemplary patients, with subpanels (i) and (ii) representing the changes in tumor cellularity (TTC) and tumor volume (TTV), respectively. Panel E shows the ROC analysis distinguishing pCR from non-pCR based on predicted TTC (blue) and measured TTC (red). Similarly, Panel F shows the ROC analysis based on predicted and measured TTV. In both Panels E and F, the larger AUC of the blue curve compared to the red curve indicates a better accuracy in predicting final pathological response. [Figure 22] 19E shows the predicted TTC time course (median and range) for patients showing the maximum range in FIG. 19E, according to various potential embodiments. [Diagram 23] Predicted TTV at end of treatment (EoT), predicted TTC at EoT, measured TTV at V3, and measured TTC at V3 for each RCB class according to various potential embodiments. [Figure 24] 20E shows measured and predicted tumor cell distribution on a central tumor slice of patient 26 according to various potential embodiments. [Diagram 25] 20E shows measured and predicted tumor cell distribution on a central tumor slice of patient 49 according to various potential embodiments. [Figure 26] 13 plots NAST overall predictions using two-scan and three-scan calibration models according to various potential embodiments.

[0015] The foregoing and other features of the present disclosure will become apparent from the following description and appended claims, taken in conjunction with the accompanying drawings, in which: The present disclosure is described with further specificity and detail using the accompanying drawings, with the understanding that these drawings depict only some embodiments according to the present disclosure and are therefore not to be considered limiting of its scope. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0016] In the following detailed description, reference is made to the accompanying drawings, which form a part of this specification. In the drawings, similar symbols generally identify similar components unless the context dictates otherwise. The exemplary embodiments described in the detailed description, drawings, and claims are not meant to be limiting. Other embodiments may be utilized, and other changes may be made, without departing from the spirit or scope of the subject matter presented herein. It is readily understood that the aspects of the present disclosure, as generally described and illustrated in the figures herein, can be arranged, substituted, combined, and designed in a wide variety of different configurations, all of which are expressly contemplated and made a part of this disclosure.

[0017] Traditionally, high spatial resolution imaging data enabling anatomical and morphological assessments have been acquired in standard care settings. Such data typically do not provide insight into the underlying physiological, cellular, or molecular characteristics of cancer, limiting their use for mechanism-based mathematical modeling. Developing more specific and quantitative measurements related to tumor biology, such as (for example) vascular status, perfusion, cellularity, hypoxia, metabolism, and proliferation, is a major endeavor in MRI research.

[0018] In various embodiments, the disclosed protocol uses two quantitative MRI modalities: dynamic contrast-enhanced MRI (DCE-MRI) and diffusion-weighted MRI (DW-MRI). DCE-MRI acquires images continuously before, during, and after the injection of contrast agent. If the data is acquired with a sufficiently high temporal resolution (e.g., in various embodiments, preferably 15 seconds or less per frame), the data can be analyzed with an appropriate pharmacokinetic model to estimate distinct tissue vascular features on a voxel-specific basis within the imaging volume. DCE-MRI is repeatable and reproducible, and the output of the DCE-MRI analysis has a statistical relationship with the response of tumors (e.g., breast tumors) to neoadjuvant therapy (NAT). In parallel, DW-MRI acquires data on water mobility within tissues associated with the number and quality of cellular barriers present, thereby providing a non-invasive readout on the cellularity of the tissue. DW-MRI is also repeatable and reproducible, and can predict the response of tumors to NAT. These two methods can therefore be used for mechanism-based mathematical modeling.

[0019] To predict individual prognosis for cancer patients (e.g., breast cancer patients), various embodiments use models that use patient-specific MRI data to initialize and constrain model parameters and expectations. That is, model parameters are calibrated to the unique characteristics of each patient. In some embodiments, a simpler logistic model that utilizes DW-MRI data may be used to estimate tumor cellularity. In certain embodiments, the logistic model may be defined in 2D, both in time and space, incorporating baseline measurements of the tumor, tumor cell migration, and mechanical properties of the tissue (e.g., breast tissue) to constrain the model's expectations for tumor growth and shape according to each individual patient's anatomy. The expectations of such models are superior to standard measures (such as Response Evaluation Criteria in Solid Tumors, RECIST) as prognostic indicators of response to therapy, but do not explicitly account for each individual patient's therapy. Thus, various embodiments extend the model to include estimates of drug delivery to each voxel via DCE-MRI, allowing for more accurate assessment of local tumor cell death due to therapy on a patient-specific basis. In various embodiments, the model can be used to identify theoretical treatment regimens that are hypothesized to be superior to the standard of care regimen that the patient actually received.

[0020] Model Details In various embodiments, the model used is based on a 3D mathematical model that includes mechanistic coupling of tissue properties to tumor growth and delivery of systemic therapy. The model is designed to be initialized with patient-specific imaging data to predict the response of cancer patients to NAT. Spatiotemporal evolution of tumor cells

number

number

number

[0021] function

number

number

number

number

number

number

[0022] In various embodiments, the second term on the right hand side of equation (1) is a response term that describes tumor growth and therapy response. Due to the nature of the data, logistic growth is defined on a voxel-by-voxel basis. Specifically, for MRI data, measurements are defined on a voxel-by-voxel basis, with a known volume for each voxel. Thus, the maximum number of tumor cells can be estimated using approximate cell size and packing density (see "Approximating Tumor Cellularity" below). With logistic growth, the carrying capacity θ is defined for each voxel as one value for all voxels, and the growth rate k is calculated in a spatially resolved manner:

number

[0023] In various embodiments, the response portion of the model also includes a term for tumor cells due to therapy. The parameter α is a global parameter that represents the effectiveness of the therapy,

number

number

[0024] Overview of the System and Method Referring to FIG. 1A, in various embodiments, a system 100 can be used to implement an exemplary protocol 160 (see FIG. 1C and functions 162, 164, 166, and 168) and the overall approach disclosed herein. The system 100 can include a computing device 102 (or multiple computing devices at the same location or remote from each other), an imaging system 132, and a motion sensor 134. The imaging system 132 can include one or more MRI scanners and / or other imaging devices and sensors capable of capturing various data, for example, to provide DW-MRI and / or DCE-MRI data. In various embodiments, the imaging system 132 and the motion sensor 134 can be integrated into one condition detection system 130. In certain embodiments, the computing device 102 (or components thereof) can be integrated with one or more of the condition detection system 130, the imaging system 132, and / or the motion sensor 134. In various potential setups, referring to FIG. 16 , one or more computing devices 102 may correspond to a server system 1600 that receives MRI data from a client computing system 1614 (which may be an imaging system or may include an imaging system) and / or a client computing system 1614 that transmits MRI data and analysis to the server system 1600.

[0025] The condition detection system 130, the imaging system 132, and / or the motion sensor 134 may be directed to a platform 136 on which a patient or other object may be positioned (to image the object, apply a treatment or therapy to the object, and / or detect motion by the object). In various embodiments, the platform 136 may be movable (e.g., using any combination of motors, magnets, etc.) to allow for positioning and repositioning of the object (such as for minor adjustments due to motion of the object).

[0026] The computing device 102 (or multiple computing devices) may be used to directly control and / or receive signals acquired via the imaging system 132 and / or the motion sensor 134. In certain embodiments, the computing system 102 may be used to control and / or receive signals acquired via the condition detection system 130. The computing device 102 may include one or more processors and one or more volatile and non-volatile memories for storing computing code and data that is captured, acquired, recorded, and / or generated. The computing device 102 may include a controller 104 configured to exchange control signals with the condition detection system 130, the imaging system 132, the motion sensor 134, and / or the platform 136, and the computing device 102 may be used to control the capture of images and / or signals via its sensors and to position or reposition an object as needed. The computing device 102 may also include, for example, an image acquisition unit 106 configured to perform image acquisition functionality 162 (steps 2-9 discussed below), a data analyzer configured to perform data analysis functionality 164 (steps 10-25 below), a model generator 110 configured to map the imaging data to a model by performing modeling functionality 166 (steps 26-36), and a tumor predictor 112 configured to perform tumor prediction functionality 168 (steps 37-40). As discussed with respect to FIG. 1C, the exemplary protocol 160 may include five components: patient population definition (e.g., step 1, not shown), image acquisition (e.g., steps 2-9), data analysis (e.g., steps 10-25), mapping of imaging data to a mathematical model (e.g., steps 26-36), and tumor prediction (e.g., steps 37-40).

[0027] The transceiver 114 allows the computing device 102 to exchange readings, control commands, and / or other data with the condition detection system 130, the imaging system 132, the motion sensor 134, and / or the platform 136 wirelessly or via wires. One or more user interfaces 116 allow the computing system to receive user input (e.g., via a keyboard, touch screen, microphone, camera, etc.) and provide output (e.g., via a display screen, audio speakers, etc.). The computing device 102 may additionally include one or more databases 118 for storing, for example, signals acquired via one or more sensors, raw and processed MRI data, and analysis results. In some embodiments, the database 118 (or a portion thereof) may alternatively or additionally be part of another computing device, either co-located or remote, in communication with the computing device 102, the condition detection system 130, the imaging system 132, the motion sensor 134, and / or the platform 136.

[0028] Referring to FIG. 1B, an exemplary tumor prediction process 150 is shown, according to various potential embodiments. The process 150 may be performed by or via one or more computing devices 102. At 152, the computing device 102 may acquire imaging data corresponding to a scan of an anatomical region of a patient having a tumor. The imaging data may be MRI data corresponding to multiple MRI scan images of the anatomical region including. The multiple scan images may include a first image set obtained by a first scan performed before administration of a therapy to the patient and a second image set obtained by a second scan performed after administration of the therapy to the patient. The therapy may include one or more therapies, such as chemotherapy, radiation therapy, and / or hormone therapy. The therapy may include administration of multiple drugs. In various embodiments, image-related data generated from the second image set may be registered to image-related data generated from the first image set.

[0029] At 154, the process 150 may include determining one or more characteristics of the tissue surrounding the tumor. The characteristics may be determined based on the MRI data acquired at 152. The tissue characteristics may include mechanical or other material characteristics of the tissue (e.g., breast tissue for breast cancer). Exemplary tissue characteristics may relate to tissue stiffness and / or tissue vasculature. In certain embodiments, the tissue characteristics may correspond to shear modulus, Young's modulus, and the like. In certain embodiments, the tissue characteristics may be known for a particular tissue type, and in other embodiments, patient-specific tissue characteristics are determined from imaging data from a patient scan. In certain embodiments, tumor segmentation may be performed to identify a region of interest (ROI) of the tumor. In various embodiments, the imaging data may include DCE-MRI data, and the tissue characteristics may be quantified based on the DCE-MRI data. The tissue characteristics may be quantified based on a pharmacokinetic model and / or a fluid dynamics model. In certain embodiments, tissue characteristics may be quantified based on the Kety-Tofts model or variations thereof.

[0030] At 156, the process 150 may include determining one or more characteristics of the tumor. The characteristics may include tumor diffusion characteristics and / or tumor growth characteristics. In various embodiments, the tumor characteristics may correspond to tumor cellularity and / or tumor vasculature. In various embodiments, the characteristics may be determined based on tissue characteristics.

[0031] At 158, the process 150 may include predicting the tumor's response to the therapy. This may include generating a score indicative of the predicted response. The score may be on any scale deemed suitable. In some embodiments, the score may be the likelihood (e.g., in percentage) that the therapy will have the intended or desired effect. Predicting the tumor's response to the therapy may include determining the effect of the therapy on the tumor cells. The effect of the therapy may be determined for the tumor cells in each voxel. If the therapy involved multiple drugs, the effect of each drug administered as part of the therapy may be determined. The predicted response (e.g., score) may be based on characteristics of the tumor, such as diffusion characteristics, growth characteristics, etc., and / or the determined effect of the therapy (e.g., each drug in the therapy) on the tumor.

[0032] Illustrative Protocol Overview Various embodiments of the protocol / procedure 160 may be divided into five main components (see FIG. 1C): identification of patients who would benefit (step 1); functions related to image acquisition 162 (steps 2-9); functions related to data pre-processing and analysis 164 (steps 10-25); functions related to mapping imaging data to a mathematical model 166 (steps 26-36); and functions related to tumor prediction 168 (steps 37-40). Each component is divided into multiple steps for clarity of presentation. The discussion of procedure 160 below provides a detailed description of each component, as well as guidance to avoid potential pitfalls and suggestions for troubleshooting.

[0033] Patient Selection: In one study, an embodiment of a protocol was developed for patients recruited from a community-based care center and eligible for NAT as a component of their care. Such patients are heterogeneous in tumor size, receptor status, age, obesity, and ethnicity. NAT (i.e., any treatment given prior to surgical intervention) is usually indicated for patients with locally advanced breast cancer and consists of one or more regimens given over the course of 4-6 months. For example, in the case of triple-negative breast cancer, the standard of care may include doxorubicin and cyclophosphamide (first regimen), followed by paclitaxel (second regimen). However, there are many variations in these protocols, as determined by the treating physician. (This is the main motivation for developing a mathematical prediction system so that treatment can be optimized on a patient-specific basis). The designation of clinical response to NAT of pathological complete response (pCR) or residual disease (non-pCR) is determined by the surgical pathology diagnosis. Specifically, pCR was defined and reported as no residual invasive disease in the breast or axillary lymph nodes following NAT.

[0034] Image Acquisition: In the study, embodiments of the protocol require the acquisition of quantitative MRI data of breast cancer patients before and during NAT to calibrate a predictive mechanism-based mathematical model designed to predict individual responses. The timing of the imaging time points before and during NAT is particularly important as they are used to calibrate, simulate, and evaluate the predictive values ​​of tumor response by the mathematical modeling system. In certain embodiments, MRI data can be acquired at four time points: 1) before NAT initiation, 2) after one cycle of NAT, 3) after two to four cycles of NAT, and 4) scan image 3 to one cycle after NAT. (Note: a "cycle" refers to the administration of a single drug or combination of drugs over a specified period of time, e.g., 2 to 4 weeks). These four time points provide data corresponding to the first cycle of each therapy regimen for patients receiving two consecutive regimens (see FIG. 2). Although three or more imaging time points are recommended, only two imaging time points are required to calibrate the model system and then directly test the modeling predictive values ​​by comparing the predictive values ​​to standard clinical measures (e.g., pathological data from biopsy or surgery).

[0035] In various embodiments, all image acquisition and patient care (imaging, oncology treatment, etc.) can be performed in a community care environment (i.e., not an academic, research-oriented medical center). However, to work within the scope of imaging in a standard care environment, certain factors must be considered. For example, in the figures and examples included in this disclosure, two scanners were used: the first at an outpatient imaging facility and the second at a community hospital that provides both inpatient and outpatient services. Although both imaging facilities performed breast MRIs as part of routine clinical practice in the study (a complete diagnostic scan is approximately 20 minutes), both facilities are located in different facilities, have different service contracts, and have different quality control guidelines. The MRI technicians at such facilities are typically responsible for patient positioning and execution of the study protocol, but may not have prior experience or expertise. Therefore, it is important to establish the repeatability and reproducibility of the required MRI measurements in each new environment, and the implementation of the acquisition protocol requires providing clear (step-by-step) instructions to the MRI technicians performing the scans.

[0036] Although working with local physicians allows reaching a broader demographic of the population, scheduling a study MRI scan requires close coordination and frequent communication regarding availability of the patient, treating oncologist, referring physician, nurses, imaging center staff, and study staff. It is not uncommon for time points or data to be missing due to equipment failure, patient health status, scheduling issues, and / or parties unwilling to donate their time. Furthermore, research-oriented nurses are not always employed in community settings. Thus, lines of communication need to be clearly defined at the beginning of the study. Despite these challenges, working in a community-based radiology environment may be easier for patients, with better access to different facilities and more convenient travel to participate in the study.

[0037] In a specific embodiment, the protocol involves the acquisition of five MRI data types in each scanning session: 1) DW-MRI, 2) B1 field maps to correct for high frequency inhomogeneity, 3) variable flip angle T1-weighted data to generate pre-contrast T1 maps, 4) dynamic high temporal resolution T1-weighted data before, during, and after injection of a gadolinium-based contrast agent (DCE-MRI), and 5) high resolution pre- and post-contrast T1-weighted anatomical scans. These MRI data types were selected to provide reliable quantitative values ​​for individual breast cancer tumors because they are well established in the literature. This imaging protocol utilizes standard sequences available on all clinical MRI scanners, eliminating the need for work-in-progress sequences (WIPs) or new sequences that are not universally available.

[0038] DW-MRI provides information about tissue microstructure by quantifying the movement of water molecules, which at 37°C number approximately 3 × 10 -3 mm 2 Diffusion rates are free to diffuse at 0 s / mm / s, but as various tissue barriers, including densely packed cells, are encountered, this diffusion rate, known as the apparent diffusion coefficient (ADC), is reduced. Estimation of the ADC requires at least two b-values ​​(coefficients that reflect the strength, duration, and timing of the diffusion encoding gradient within a scan) (in this protocol, we use 0 s / mm, which is commonly utilized in breast tissue). 2 , 200s / mm 2 , and 800s / mm 2 Three b values ​​were used: 2 DW-MRI acquired at 100 s / mm or greater may result in low signal-to-noise ratios (SNRs) that adversely affect ADC estimates, while low b-values ​​(100 s / mm or greater) may result in low signal-to-noise ratios (SNRs) that adversely affect ADC estimates. 2 Smaller than 100 nm) may be affected by tissue perfusion, where blood flow in the smallest vessels mimics diffusion, thereby altering image interpretation.

[0039] Since different tissues exhibit different T1 relaxation times, T1 mapping provides a means to distinguish between tissue types (e.g., fat, muscle, parenchyma, and / or tumor) and provides native T1 values ​​necessary for downstream pharmacokinetic analysis of the DCE-MRI data. Various embodiments may use a standard approach for clinical breast T1 mapping, which involves the acquisition of multiple T1-weighted images at variable flip angles. Various embodiments acquire images at 10 flip angles ranging from 2° to 20° (in 2° increments) to estimate typical breast tissue T1 values ​​(more flip angles provide more data points for better curve fitting, and this range allows accurate estimation of various tissues within the breast from fat to tumor). However, this approach is susceptible to inhomogeneity in the high-frequency B1 field used to "tilt" the magnetization at various flip angles, which can lead to inaccurate estimation of native T1. To address this issue, during acquisition of the T1-weighted images used to map the T1 parameters, a B1 map is acquired to quantify and correct any spatial deviations of the nominal flip angle. In various embodiments, other T1 mapping approaches include inversion sequences (which are less affected by B1 inhomogeneity) and saturation recovery sequences. These methods are the "gold standard" for T1 calculation, but the time required to acquire these sequences in multislice or 3D may be clinically prohibitive and therefore may not be incorporated into protocols.

[0040] In DCE-MRI, a paramagnetic contrast agent is injected into the bloodstream through a peripheral vein. It travels throughout the circulatory system and may extravasate into the tumor, resulting in a decrease in T1 relaxation time and a corresponding increase in signal intensity in the T1-weighted images. DCE-MRI data is acquired by collecting T1-weighted images before, during, and after delivery of the contrast agent. The DCE-MRI data can then be analyzed to segment different tissues with different contrasts and to extract measures that characterize the pharmacokinetics of the contrast agent (details are provided below under "DCE-MRI Analysis"). In various embodiments, the acquisition parameters for the DCE-MRI measurements were selected to provide a time resolution <10 s (7.27 s) for accurate estimation of pharmacokinetic parameters. In various embodiments, the protocol can be adjusted for other tissue types, keeping in mind that an appropriate flip angle that minimizes saturation effects of the contrast agent can be selected for optimal DCE-MRI results, which may vary depending on the tissue imaged (i.e., breast vs. brain).

[0041] Data Pre-Processing and Analysis: Image processing begins with quality control, image correction, and image registration, then proceeds to the extraction of tumor-specific characteristics and quantitative descriptions of each tumor's cellular density and vasculature. Various embodiments involve methods that include analysis of quantitative MRI data to return maps of quantities reporting blood flow and water diffusion, as well as segmentation via clustering techniques.

[0042] Tumor Segmentation: In various embodiments, a tumor region of interest (ROI) is obtained for each patient and scan session to analyze and process the data. Some embodiments may use an expertly drawn ROI for the tumor volume. However, if a conservatively drawn "bounding box" (i.e., a hand-drawn polygon that surrounds the tumor but not its specific outline) is provided, augmentation-based thresholding may be used to determine the tumor boundary from the DCE-MRI data. The threshold is a value selected such that any voxel with a signal intensity above that threshold in the post-contrast image is considered to be part of the tumor. Because thresholding techniques may require manual adjustment of each patient scan and additional information to define the patient-specific threshold (which varies with the type and amount of contrast), it is desirable to use an automated approach. Various embodiments use a fuzzy c-means (FCM) clustering algorithm. The FCM algorithm outputs the probability that a voxel is tumor or non-tumor based on the DCE-MRI contrast pattern. Because FCM clustering does not divide voxels into clusters, it is more adjustable compared to other "hard" clustering methods (such as k-means clustering). See Figure 3 for a representative image of generating ROIs using FCM.

[0043] Registration techniques: Various embodiments employ an approach to imaging-based modeling that requires all image sets for each patient to be registered to one common spatial coordinate system. That is, the images are registered to each other. Note that all registration processes do not fully preserve voxel information even when rigid registration is used due to multiple interpolations and resampling. To achieve image alignment, various embodiments may perform two types of registration: intra-visit registration (aligning all data collected within one scanning session) and inter-visit registration (aligning each of the data sets across all scanning sessions for each patient). In various embodiments, intra-visit registration is performed to correct for patient motion during the scanning session and is achieved through rigid registration before calculating quantitative parameters from the data.

[0044] In the study, for each patient, all image data sets are registered to a common space across time (between visits) via a constrained non-rigid registration algorithm that preserves tumor volume at each time point. If MRI data is obtained at four time points, various embodiments may choose to register scans 1, 3, and 4 to scan 2 (the target) because scan 2 is not an end scan and is acquired early enough that the patient (likely) has tumor volume present to help guide the alignment. In some embodiments, the registration algorithm may include or consist of a rigid registration of the tumor ROI, followed by a deformable b-spline registration with a stiffness penalty on the tumor region. This stiffness penalty is imposed to maintain the tumor volume, size, and shape across all scan times. With a fully deformable registration, the tumor ROI in scans 1, 3, and 4 can be deformed to match the tumor ROI in scan 2. See FIG. 4 for a comparison example of various registration results.

[0045] DCE-MRI Analysis: In various embodiments, a pharmacokinetic model of the contrast agent is used to analyze the DCE-MRI data and derive quantitative parameters of vascular perfusion and tissue volume fraction. For example, an extended Kety-Tofts model (or a variation thereof) is used to perform quantitative analysis of the DCE-MRI data. In various embodiments, the approximately 7.27 s time resolution of the DCE-MRI data provides sufficient SNR and time sampling to apply the extended Kety-Tofts model. However, various embodiments include evaluation of voxel time course fits to determine which model adequately captures the time course behavior of a particular dataset. Pharmacokinetic modeling requires characterization of the time rate of change of contrast agent concentration in the feeding artery, i.e., the arterial input function (AIF). Various embodiments estimate the AIF from population mean signal intensity time courses extracted from the axillary artery. Additionally, various embodiments calculate the bolus arrival time (BAT) and shift the population AIF on a voxel-by-voxel basis to align the enhancement time of the AIF with that of individual voxels, allowing for improved fitting and more accurate parameter estimation.

[0046] In various embodiments, the reference region model is an alternative approach to quantitative DCE-MRI analysis, eliminating the need for direct measurement of AIF (which requires additional tissue segmentation work), where the tumor enhancement curve is compared to that of a reference region, such as pectoralis major muscle tissue. In various embodiments, simpler methods for analyzing DCE-MRI data, such as calculation of signal enhancement ratios and area under the signal intensity time course curve, can provide semi-quantitative measures of vascular properties that have proven valuable in distinguishing benign from malignant lesions and predicting disease recurrence.

[0047] Approximation of tumor cellularity: In various embodiments, ADC is calculated from DW-MRI data, which represents the rate at which water molecules diffuse into tissue. It has been shown to approximate the cell density of tissue. ADC values ​​may be calculated for each voxel via standard methods, and there is an inverse correlation between the measured ADC and tumor cellularity. However, the causes of ADC changes may vary, since many other factors in addition to cellularity (e.g., cell membrane permeability, cell size, tissue tortuosity) may also affect the measured ADC. Therefore, the approach of using ADC to evaluate cellularity may be an approximation, and there may be other methods to estimate ADC with less ambiguity to improve the outcome of the application of the protocol.

[0048] Mapping the imaging data to the mathematical model: In various embodiments, after generating quantitative maps from each MRI data type, the maps are registered across all scanning sessions for each patient so that all imaging data is aligned to a common image space (i.e., inter-visit registration). Once aligned, a tissue domain (such as a breast domain) is defined within which each individual patient's data is used to calculate the physical properties of each patient's tumor that are utilized by the mathematical model. Steps for generating tumor properties include calculating tumor cellularity, defining a mask that outlines fibroglandular and adipose tissue, approximating drug distribution, and estimating summary measures for each tumor across all scans.

[0049] In various embodiments, the mathematical model uses patient-specific characterization of breast (or other tissue) anatomy, cellularity, vascular features, and approximate drug distribution. If tumor response could be reliably predicted using a mathematical model initialized and constrained by individual patient non-invasive imaging data collected early in the course of therapy, oncologists could potentially intervene and modify therapy on a patient-specific basis. Additionally, such models could be used to optimize a patient's therapy regimen using more robust mathematical methods such as optimal control theory.

[0050] Tumor prediction: In various embodiments, the quantitative maps are then used to initialize and calibrate tumor cell proliferation, drug efficacy, and cell migration in a mathematical model. That is, each patient's unique imaging series is used to parameterize or identify growth and response parameters specific to that individual patient. Once parameterized, the model can be advanced for patient-specific prediction of the spatiotemporal evolution of tumor cellularity, allowing prediction of treatment response that can be directly compared to each patient's observed outcome.

[0051] Protocol implementation Various embodiments of this protocol may be performed by a multidisciplinary team with experience and expertise across multiple disciplines, including both advanced non-standard care MRI data acquisition and analysis, image processing (including segmentation and registration), and numerical solution of partial differential equations (PDEs) in a clinical setting. When performed in a community setting, its implementation may involve coordination with local healthcare providers.

[0052] In various embodiments, the procedures described below address the challenges presented for incorporating cancer (e.g., breast cancer) MRI data into predictive mathematical models. Thus, these methods are applicable to data collected in a community environment, multiple imaging facilities, within the constraints of a specific standard of care, and to a heterogeneous patient population. Repeatability and reproducibility of this quantitative MRI protocol in a community-based imaging center has been established. Furthermore, embodiments of the protocol utilize several strategies to reduce bias and dependency on user interaction, including detailed data acquisition strategies and automated or semi-automated computational algorithms. Certain tasks benefit from expert assessment (e.g., radiologists outlining tumor regions of interest), while others rely on operator input (such as positioning the field of view for MRI acquisition). Embodiments of the disclosed approach include automated quality checks and data quantification to reduce overall user / operator influence.

[0053] Protocol Application In various embodiments, the protocol is well suited to predict the response of cancer patients undergoing NAT as a component of clinical care. Although the details of the protocol presented here are specific to breast cancer, the method is generally applicable to any solid tumor for which the necessary data are accessible. There are also applications of the disclosed approach in preclinical settings.

[0054] In various embodiments, in addition to mechanism-based mathematical modeling for which the protocol is well suited, researchers who have previously collected suitable imaging data may use portions of the protocol for data processing and analysis to obtain (for example) longitudinally aligned quantitative maps of tumor features. The data returned upon completion of step 24 can be used for more conventional statistical studies to assess long-term changes in tumor characteristics and differentiate between responders and non-responders. Such data can also provide input data for analyses in the fields of radiomics, "habitat imaging," as well as applications to artificial intelligence-based models.

[0055] considerations The time sampling requirements of the DCE-MRI protocol constrain the achievable spatial resolution of the data that can be collected while limiting noise. All images are acquired in the sagittal plane. While the transverse plane is the standard care choice with a bilateral field of view (FOV), embodiments of the disclosed imaging protocol use sagittal slices to obtain better resolution across the diseased tissue in a minimal amount of time. For example, a standard care T1-weighted contrast acquisition requires 80-90 seconds, whereas the same scan in an embodiment of the disclosed approach takes less than 10 seconds. At this coarse spatial resolution, important anatomical features such as small nutrient vessels are potentially missed. This in turn may impact modeling strategies aimed at incorporating estimates of nutrient / oxygen delivery and distribution of systemic therapies. However, it would be envisioned that faster acquisitions at the spatial resolution and FOV typically acquired in a standard care setting would improve results.

[0056] Since some degree of error is expected in quantitative measurements, various embodiments employ best practices to encounter and address physiologically unlikely results, such as when interpreting diffusion weighted data and pharmacokinetic analyses. This is particularly important because mathematical models project / propagate errors derived from the data employed for calibration and prediction (which is addressed in steps 37-40 of the Tumor Prediction section). MRI data and corresponding quantitative analyses are ultimately approximations of breast cancer characteristics at the tissue and cellular scale, and alternative methods and procedures may be used to reduce error and ambiguity in these quantities.

[0057] It is emphasized that the protocols and approaches presented herein are examples, and alternative embodiments may employ one of many other possible modeling formulations applicable to the datasets generated by this pipeline, and that a variety of other cancer modeling modalities exist. Also, while the disclosed protocol embodiments are constructed to accommodate a heterogeneous patient population, there are certain exclusions that apply to make the mathematical modeling as practical as possible. For example, one study excluded stage I and stage IV tumors due to the fact that stage I tumors may be too small (diameter <2 cm) to be reliably measured by MRI techniques, and stage IV tumors may require alternative modeling techniques to account for tumor invasion and metastasis. Also, stage I and stage IV patients do not usually undergo NAT. Stage I tumors are treated with surgery and radiation, and stage IV tumors are treated with palliative care. Furthermore, it is noted that many programming languages ​​and numerical schemes (such as the finite element method) are available to determine the solution of the PDE, and for those less interested in deriving numerical codes, there are specific software programs (such as FEniCS and MATLAB's PDE solver) that assist in the implementation of PDEs.

[0058] research materials Subjects: Females over 18 years of age presenting with intermediate to high grade invasive primary breast cancer and considering NAT as a component of clinical care. If the patient's disease is stage II and stage III cancer, patients with disease of all subtypes and treatment regimens (including immunotherapy, targeted therapy, and cytotoxic therapy) may be included. Patients with a history of renal disease, abnormal creatinine or estimated glomerular filtration rate, pregnant or lactating patients must be excluded. Also excluded are any patients with any MR incompatible ferromagnetic material, acute illness, and / or technically incapable of MRI (e.g., due to breast volume or obesity). In the exemplary data presented here, health information related to each participant's disease and MRI scan images was obtained. This information includes clinical test results, medical imaging reports, and diagnostic and treatment codes.

[0059] Reagent: Gadolinium-based contrast agent (e.g., Multihance (Bracco, Monroe Township, NJ) or Gadovist (Bayer, Leverkusen, Germany)). Used for contrast scans.

[0060] Equipment: Powered injector for contrast administration (e.g., Medrad, Warrendale, PA), MRI scanner, and personal computer or server. For data processing and analysis, as well as model simulation, a personal computer may be able to perform model calibration via the software described below, but in various embodiments, a server as specified here with parallelized scripts for computational efficiency may be preferred. 40 nodes. CPU per node: 2 / 8 Xeon E5-2680 2.7GHz (Turbo, 3.5) 1 / 61 Xeon Phi SE10P 1.1GHz. Memory: 32GB per node.

[0061] MRI Scanner Setup: Breast MRI may be acquired using a 3T scanner equipped with an 8-channel or 16-channel dual breast receive coil. The study considered here used a Siemens scanner (Siemens Healthcare, Erlangen, Germany) equipped with a Sentinelle coil (Invivo, Gainesville, FL). Other 3T scanners (Philips, GE) and breast coils may require slight adjustments to reproduce the pulse sequences described in section 3.2, but since these are standard sequences, the scanners have the capability to collect comparable imaging protocols. After scanning, Digital Imaging and Communications in Medicine (DICOM) files are stripped of protected health information (PHI), labeled with an assigned study identifier, and uploaded to a firewall-protected server.

[0062] Software setup: REDCap (https: / / www.project-redcap.org / ) or similar. Clinical information regarding patient demographics, diagnoses, treatments, and surgical outcomes may be communicated to the study team and stored in REDCap, a secure, HIPAA-compliant web application for building and managing online databases. REDCap is a widely available and used service for managing single and multi-center clinical studies.

[0063] MATLAB (MathWorks, Natick, Massachusetts) or similar. In various embodiments, MATLAB software may be used for much of the process due to its portability and ease of use across a diverse group of biomedical researchers, including experimentalists, computational scientists, engineers, physicists, mathematicians, and medical professionals. However, alternative software programs are available for data processing, and with sufficient programming background the functionality described throughout the disclosed protocols can be replicated in those environments.

[0064] Elastix (https: / / elastix.lumc.nl / ). Elastix is ​​a toolbox for rigid and non-rigid registration of images in 3D. Elastix is ​​an open source toolbox and its functions can be run from the command line via MATLAB, which allows it to be integrated into existing data processing scripts.

[0065] Study Procedure Patient selection, timing ~15 min.

[0066] Step (1) Ensure that patients are eligible for the study. Patients will be females over 18 years of age presenting with intermediate- to high-grade invasive primary breast cancer and considering NAT as a component of their clinical care. Each patient will confirm with their treating oncologist that they have no history of renal disease and have normal creatinine and estimated glomerular filtration rate within 30 days of the imaging study. Exclude pregnant and lactating women. Also exclude any patients with any MR-incompatible ferromagnetic material, patients with acute illnesses, and / or patients for whom MRI is technically not possible (e.g., due to breast volume or obesity).

[0067] Step (2) Obtain informed consent from patients. Patients must be asked for consent to participate in the study. Quantitative MRI, PHI, and therapy regimens will be obtained for patients diagnosed with intermediate to high grade invasive breast cancer who are eligible for NAT as a component of clinical care at a community care oncology facility. The study population will consist of women over 18 years of age presenting with primary breast cancer and considering NAT as a component of their clinical care. Importantly, for stage II and stage III cancer, patients with disease of all subtypes and treatment regimens (including immunotherapy, targeted therapy, and cytotoxic therapy) may be included.

[0068] Image acquisition, timing is approximately 1.5-2 hours total per session (from patient arrival to departure), 40 minutes (patient positioning and removal from the scanner), and approximately 25 minutes (MRI scan). The 40 minutes allows for sufficient time for the patient to be positioned in the breast coil (prone position), scanned, and exit the scanner. Preparation for the MRI examination is typical of a standard of care examination, including placement and removal of an intravenous line to administer contrast.

[0069] Step (3) Insert an intravenous line. A short peripheral intravenous catheter (20-22 gauge) is inserted into the antecubital or forearm area. Correct positioning of the catheter tip should be checked for venous reflux by withdrawing blood and flushing with saline. Please refer to Table 2 for a complete summary of the imaging parameters described in the next steps.

[0070] Step (4) Establish the FOV. Obtain several localizer scan images at the beginning of the MRI session. First, obtain a localizer scan image by acquiring a three-plane (cross-sectional, sagittal, and coronal) scout image series covering both breasts. Adjust the FOV to ensure coverage of the tumor within the affected breast without changing scan parameters that would affect resolution or scan time. Since all patients have been diagnosed with breast cancer and the location of the tumor is known, identifying the approximate location of the lesion (e.g., 7 o'clock, right breast) from the referring physician's office prior to the MRI will aid in locating the FOV. Often, metal-induced sensitivity artifacts around the biopsy clip / marker can be used to identify the approximate tumor center location on the localizer scan. While finding the center of the FOV, adjust the anterior-posterior coverage to include the axillary artery as much as possible (to allow for characterization of the AIF). Acquire a second localizer scan image with a sagittal-only sequence with 30 slices, each 5 mm thick (15 cm coverage from right to left of the patient). The imaging volume is placed at the center of the tumor and includes all or as much of the tumor as possible.

[0071] Step (5) Acquire DW-MRI data. Acquire DW-MRI over 10 slices with a thickness of 5 mm and no slice gap. The following parameters may be used: repetition time / echo time TR / TE=3000 / 52 ms, flip angle 90°, matrix 128×128 (256×256 mm). 2 100 s / mm 2 ), with automated calibrated partial parallel acquisition (GRAPPA) acceleration. The study included spectrally selective adiabatic inversion recovery (SPAIR) fat suppression. To allow for approximately equal SNR ratios at all b-values, a monopolar, single-shot spin-echo, echo-planar imaging sequence in a 3D diagonal diffusion-encoding orientation was performed at five values ​​between 0 and 200 s / mm 2 . 2 Average of 6 times, 5 values ​​800s / mm 2 Average 18 acquisitions at 100x the mean. If using other scanners, use sequences equivalent to GRAPPA (e.g., ARC (Automatic Calibration Reconstruction for Cartesian Imaging)) for GE and image-domain SENSE (Sensitivity Encoding) for Phillips.

[0072] Step (6) Map the B1 field. Use the Siemens TurboFLASH sequence with preconditioning RF pulses to map the B1 field with the following acquisition parameters: TRITE=8680 / 2 ms, flip angle=8°, matrix 96×96, slice thickness 5 mm. As the B1 mapping protocol includes a slice gap, perform two acquisitions with a slice gap of 5 mm to cover the same FOV as the measurements above without losing coverage in the slice direction. This pulse sequence is only available as a standard product sequence on Siemens scanners at the time of printing. If a Philips or GE Healthcare scanner is used, an approximate B1 map can be calculated by assuming a uniform T1 over the adipose tissue and interpolating the remaining images, as previously published.

[0073] Step (7) Acquire high-resolution T1-weighted images. Use a VIBE (Volumetric Interpolated Breath-hold Examination, no breath-hold required) sequence with the following acquisition parameters: TRITE=5.3 / 2.3 ms, flip angle 10°, acquisition matrix 256×256, slice thickness 1 mm (96 slices), GRAPPA acceleration 2 in phase encoding direction, and SPAIR (Spectral Selective Adiabatic Inversion Recovery) fat suppression. Next, acquire images of pre-contrast T1 maps without fat suppression using a 3D gradient echo, FLASH (Fast Low-Angle Shot) sequence T1-weighted image, also known as SPGRE (Spoiled Gradient Recovered Echo) sequence, at 10 flip angles with the following parameters (2, 4, 6, … 20°): TRITE=7.9 / 2.71 ms and GFtAPPA acceleration factor 3 in phase encoding direction. Use an acquisition matrix of 192×192×10 on a sagittal square FOV (256 mm2) with a slice thickness of 5 mm. This pre-contrast T1-weighted MRI scan with fat suppression is necessary not only for anatomical visualization purposes but also for the radiologist to identify the tumor ROI. The pre-contrast T1 map is necessary for pharmacokinetic modeling of the DCE-MRI data, which will be described in the next step.

[0074] Step (8) Acquire a dynamic set of T1-weighted VIBE (no breath-hold) images (this is the DCE-MRI acquisition). Use the following acquisition parameters: TRITE=7.02 / 4.6 ms flip angle=15°, matrix 192×192, 10 slices each 5 mm thick, and GRAPPA acceleration factor 2 in the phase-encoding direction, resulting in a time resolution of 7.27 seconds. (Note that with further development of fast acquisition, slices thinner than 5 mm may be accessible without sacrificing too much SNR.) The equivalent sequences on the GE and Philips scanners are FAME (Fast Acquisition with Multilayer EFGRE3D) and THRIVE (T1W High Resolution Isotropic Volumetric Examination), respectively. Start this DCE-MRI data acquisition. After 1 minute, while continuing to acquire DCE-MRI data, start administration of gadolinium-based contrast agent via a power injector at the dosage recommended in the product insert and at a flow rate of 2 mL / sec through the IV catheter placed in step 3. After contrast administration, a saline flush (20 mL) is (again) administered through the IV catheter placed in step 3 at a flow rate of 2 mL / sec. Note that different contrast agents or injection rates may result in different pharmacokinetics. It is important to use a consistent contrast agent and injection protocol if pharmacokinetic parameters are to be calculated, used, or compared across patient populations. It is also important to record the amount of contrast administered to be used as a reference at subsequent visits.

[0075] Step (9) Acquire a post-contrast high-resolution T1-weighted image. Use a VIBE sequence with the following acquisition parameters: TRITE=5.3 / 2.3 ms, flip angle 10°, acquisition matrix 256×256, slice thickness 1 mm (96 slices), GRAPPA acceleration 2, and SPAIR (Spectral Selective Adiabatic Inversion Recovery) fat suppression. Post-contrast T1-weighted MRI scan images with fat suppression are necessary not only for anatomical visualization purposes but also for the radiologist to identify the tumor ROI.

[0076] Image analysis and timing takes approximately 2-3 hours.

[0077] Step (10) Upload and store data. Store as non-anonymized DICOM files on a Picture Archiving and Communication System (PACS) server associated with the radiology clinic where the scans are obtained. Store a separate copy of the DICOM files with PHI replaced with unique subject identifiers and upload these anonymized files to a firewall-protected server. For each patient, images are preferably processed and analyzed as a set across all visits. See Figure 5 for a flow chart of all steps in this data analysis pipeline.

[0078] Step (11) Organize the data. Copy the DICOM data into individual patient-visit directories to ensure that original copies of the data acquired from the scanner are preserved. Import the DICOM data into the MATLAB workspace using the built-in MATLAB functions dicomread and dicominfo. Place the DICOM slices in ascending order of slice location in the scanner (i.e., spatial order, not acquisition order; slices acquired in an interleaved manner need to be sorted by spatial location). Save the final images for each type of MRI scan as a matrix of MATLAB structures, along with the header information for each DICOM (pulled in using dicominfo). Save the structures for each scan as .mat files with floating-point accuracy to integrate with the image processing and analysis steps of the pipeline. Save the DW-MRI, variable flip angle T1-weighted, and DCE-MRI data as 4D matrices, where the fourth dimension represents each b-value, flip angle, and repetition, respectively. Save the final MATLAB structures in each patient / visit directory.

[0079] Step (12) Verify slice position. DICOM header information is used to ensure that each scan aligns with that of the DCE-MRI data to ensure the same FOV is analyzed across scans for each patient and visit. Several types of errors can occur in the scanner, including a shifted FOV (compared to the previous scan) or slices with slice thickness offset by a value different than 5 mm.

[0080] Step (13) Upsample the DW-MRI and B1 map data to match the spatial resolution of the DCE-MRI. Use nearest neighbor interpolation on the 2D grid data (interp2, MATLAB). Save the resulting interpolated data in a .mat file in the MATLAB structure for that patient and visit. In this protocol, the B1 map and DW-MRI data are acquired at a lower spatial resolution than the DCE-MRI data (due to time constraints and SNR considerations). These scan images need to match the resolution of the target images for intra-visit registration.

[0081] Step (14) Align the DW-MRI data with the DCE-MRI data. Use the b=0 image (b=200 or 800 s / mm 2 Compared with the data, b=0s / mm 2 (Because the SNR of the data is high). We used the MATLAB rigid body registration algorithm with the imregtform function, b=0s / mm 2 Align each slice of the diffusion-weighted scan image to the corresponding slice of the first iteration of the DCE-MRI data. Use the MATLAB function imwarp to apply the transformation obtained from the function's output to the MRI images with b = 200 and 800 s / mm. 2 Apply this to the DW-MRI data in . Save the resulting 4D matrix of the DW-MRI data aligned to the DCE-MRI data as a separate .mat file.

[0082] Step (15) Align variable flip angle T1-weighted MRI to DCE-MRI data. Repeat step 14 but use variable flip angle T1-weighted MRI data instead of DW-MRI data. Align each slice of each T1-weighted image to the first iteration of DCE-MRI data (imregtform, imwarp, MATLAB). Save the variable flip angle data as a 4D matrix and export as a .mat file.

[0083] Step (16) Align B1 Map. The B1 mapping sequence utilized in this study outputs a proton density weighted image and a calculated map of the estimated flip angle experienced by each voxel in the image. To align the B1 map to the DCE-MRI data, due to the higher SNR of the proton density weighted image compared to the flip angle map, align each slice of the proton density weighted image to the first iteration of the DCE-MRI data (imregtform, imwarp, MATLAB - see step 14). Apply the resulting geometric transformation to the flip angle map. Save the resulting data and export it as a .mat file.

[0084] Step (17) Correct for motion in the DCE-MRI scan images. Align each slice of each repetition of the DCE-MRI sequence to the corresponding slice of the first repetition (imregtform, imwarp, MATLAB - see step 14). Save the resulting data as a 4D matrix and export as a .mat file. Breast motion due to breathing and patient movement can occur during the DCE-MRI acquisition.

[0085] Step (18) Generate the tumor ROI. Manually draw the bulk ROI (this can be done by a board-certified radiologist) to unobtrusively outline the lesion under investigation. Then, apply FCM clustering to the voxels within the drawn ROI, using MATLAB's fcm function, with class number set to 2 (one each for diseased and non-diseased tissue types). After the binary mask has been identified, post-process it to fill holes (i.e., zeros in the mask surrounded by ones) via the MATLAB function imfill. Additionally, a 8.45 mm 3 Remove regions smaller than 1×1×1 voxels (i.e., 1×1×1 voxels). Save the resulting ROIs in a .mat file within the MATLAB structure for that patient and visit.

[0086] Step (19) Verify successful registration and identify any artifacts (due to silicon implants, cardiac motion, etc.). Manually visualize each set of scan images before (steps 11-13) and after registration (steps 14-17). Remove patient datasets containing substantial artifacts (image deformation) from the study. To evaluate the tumor segmentation results, visualize the ROIs overlaid on background-subtracted DCE-MRI data and high-resolution post-contrast T1-weighted images.

[0087] Step (20): Calculate the ADC map by fitting the DW-MRI data to Equation (4). S(b) = S0 exp(-ACD b) (4) where S(b) is the signal intensity in the presence of a diffusion gradient of strength b, S0 is the signal intensity in the absence of a diffusion gradient, and b is the strength of the diffusion gradient. 2The ADC map is fitted voxel-wise to the signal intensity from the DW-MRI data of 1000 s / mm. In particular, for data sets with only two b values, the ADC values ​​can be calculated directly (i.e., no fitting). Alternatively, the ADC map is determined via the MATLAB function regression using ln(S(b)) and a linear regression with b, where b=0 s / mm. 2 The data is acquired but not used for the ADC calculation in this protocol since the low b-values ​​are affected by tissue perfusion, but is used to register the diffusion weighted data (details in step 14). Save the resulting ADC map in a .mat file within the MATLAB structure for that patient and visit.

[0088] Step (21) Calculate the B1-corrected T1 value for each voxel. The measured multi-flip angle and signal intensity data are fitted to the gradient-recalled echo signal S (Isqcurvefit, MATLAB), and the formula is as follows:

number

[0089] Step (22) Calculate the population AIF. For each patient dataset, generate a difference image from the contrast data, where the average pre-contrast kinetics is subtracted from the average post-contrast kinetics. For each scan image, identify the axillary artery in the difference image. Manually select one voxel in the axillary artery as the "seed" for a 3 × 3 kernel centered on the seed. Save the location of the seed in a vector. For adjacent slices, select a 5 × 5 ROI centered on the voxel that corresponds to the location of the seed voxel. For each voxel in the 5 × 5 ROI, define a 3 × 3 window centered on each voxel. Compare the time course of the average signal intensity of each window (25 windows) to the time course of the average signal intensity of the kernel of the seed voxel using a correlation coefficient (MATLAB function corr). For the voxel window with the highest correlation with the seed, save the location of the voxel in the vector and set that voxel as the new seed. Repeat this process for all slices. For the voxels stored in the vector, remove any voxels that meet the following criteria (i)-(iii): (i) Maximum signal intensity does not occur within the first (approximately) 45 seconds of post-contrast dynamics. (ii) The maximum signal intensity is < (approximately) 20 times the standard deviation of the first three precontrast dynamics. (iii) The average signal strength over the last (approximately) 120 seconds is > (approximately) 40% of the maximum.

[0090] For voxels remaining in the vector after the three bulleted steps, average them at each time point to obtain one average time course. Save the resulting time course as an individual AIF in a .mat file in the MATLAB structure for that patient and visit. To generate a population AIF, average the resulting individual AIFs across all available patient sets. This population AIF calculation utilizes a semi-automated algorithm developed by Li et al. The population AIF calculation can be updated whenever new patient data is acquired.

[0091] Step (23) Determine the BAT by fitting the time course of signal intensity for each tumor voxel to a half-logistic function (using Isqcurvefit in MATLAB).

number

number

[0092] Step (24) converts the signal intensity measured in the DCE-MRI experiment to a concentration of contrast agent to enable pharmacokinetic modeling as described below. The signal intensity measured in the DCE-MRI experiment is expressed by equation (5), where T1 ≡ T1 (t) is how the measured T1 value at time t varies with the concentration of contrast agent according to the following equation: 1 / T1(t)=r1C t (t)+1 / T 10 (8) where r1 is the relaxation constant specific to the contrast agent, and T 10 is the native T1 (from the pre-contrast T1 map obtained in step 21). Finally, C t (t) is the BAT shift AIF (i.e., C from Eq. (7) p (t)), and fitting the concentration time course of each voxel to the extended Kety-Tofts model gives:

number

number

number

number

number

number

number

number

number

number

[0093] Step (25) Generate an enhanced anatomical image. The contrast is enhanced by applying a local statistics-based transfer function to each voxel of the difference anatomical image (average pre-contrast image minus average post-contrast image). Specifically, a transfer function is used:

number

[0094] Mapping of imaging data to a mathematical model, timing is approximately 1 day (steps 28-30 are rate-limiting steps due to the need to use all patient data in the study)

[0095] The following steps (i.e., steps 26-36) describe how to transform the processed MRI data into a single modeling domain and derive the physical quantities relevant for model simulation starting from inter-visit registration.

[0096] Step (26) Manually compare slices across different visits and evaluate each patient's anatomy to determine rough slice alignment. (We note that this is the first time we are comparing data from different visits, i.e., we cannot complete this step or proceed to the next step until we have data from at least two visits.) Using the augmented images, we exploit the patient's anatomy (visible structures within the tissue and blood vessels) as well as unique characteristics in the slices with the tumor to determine which of the 10 slices for each visit correspond to slices across all visits. Performing this initial alignment improves the ability of the registration algorithm to converge to a solution. Save the corresponding slices in a .mat file in the MATLAB structure for that patient and visit. This is performed before applying the registration algorithm.

[0097] Step (27) Convert the files to .mhd files. This can be achieved using the MATLAB Medical Imaging Toolbox function write_mhd. To use this function, save the tumor ROI and the enhanced anatomical images of the slices determined in step 26 for registration. Additionally, define a parameter file for each patient and scan image to be registered (detailed instructions on the format of these files can be found at http: / / elastix.bigr.nl / wiki / index.php / Parameter_file_database). Once these files are generated, the registration is performed by calling the elastix function, which first rigidly aligns the tumor ROI with the corresponding anatomical image, and then performs a b-spline registration with stiffness penalty on the tumor ROI. Once the deformation fields are generated, apply the deformation fields to all remaining maps and images of each patient set by calling the transformix function. This step is performed before running Elastix.

[0098] Step (28) Consider Range of Stiffness Penalty Weights: To select an appropriate weight for the stiffness penalty, apply the above registration function (step 27) to a range of values ​​for a subset of patients, where Scan Image 2 is the "target" image and Scan Image 1 is the image to be registered.

[0099] Step (29) Evaluate Stiffness Penalty Weights: The following metrics are used to evaluate the results of the various penalty weights (from step 28): Metric 1: Similarity or normalized mutual information between the reference image and the alignment target image Metric 2: Sum of the Jacobian determinants of the final deformation field within the tumor region Metric 3: Consistency between distributions of quantitative parameters before and after registration

[0100] We choose weights that minimize metric 2 while maximizing metrics 1 and 3. For metric 3, we divide the histograms (probability density functions) of the original target and alignment map into the same 100 segments (i.e., convert the histograms into vectors), then normalize and calculate their dot product as the similarity.

[0101] Step (30) Select stiffness penalty weights. The resulting registered images are compared (step 28) to identify an "acceptable range" for the studied configuration parameter (i.e., the stiffness penalty term weight). Within this tolerance range, the performance of the registration is reasonable and similar. Outside the range, the registration fails or the values ​​within the tumor ROI change significantly (as per the metrics described in step 29). Note that the acceptable range not only varies from patient to patient, but also from individual patient visit to individual patient visit. Therefore, choose values ​​that fall within the acceptable range for all patients used in this parameter study.

[0102] Step (31) Align Images. Using the details of step 27 and the penalty weights derived in steps 28-30, align the augmented images of scans 1, 3, and 4 to scan 2 (the target). Apply the resulting deformation fields to all corresponding patient parameter maps. Save the resulting alignment data in a new .mat file (separate from the previous data analysis files) within the MATLAB structure for that patient and visit.

[0103] The following steps describe the definition and calculation of modeling quantities.

[0104] Step (32) Define the modeling domain. Manually draw a breast ROI to define the domain for mathematical modeling. Using the augmented image, carefully outline the breast area outside the chest wall and within the skin of the organ for each slice using the MATLAB function roipoly. Save the outlined shape as a binary mask with zeros outside the breast ROI. If the tumor is not near the chest wall, the segmentation can simply remove the breast area. Save the resulting breast mask in a .mat file within the MATLAB structure for that patient and visit.

[0105] Step (33) Calculate tumor cells per voxel. Using the MATLAB function imfilter, smooth the aligned ADC map for each slice using a Gaussian filter of size 3 × 3 voxels. Via established methods, the ADC value of each voxel in the tumor (segmented using the method above) is calculated for each location.

number

number

number

number

number

number

[0106] Step (34) Segment fibroglandular and adipose tissues. Generate initial masks of fibroglandular and adipose tissues using the coordinated imaging data and two-class K-means clustering (MATLAB function kmeans). To suppress noise and remove voxels containing both tissues, apply a second K-means clustering to the adipose regions segmented by the first clustering to erode edges. Finally, remove small "islands" (<10 connected voxels) from the fibroglandular tissue mask using the MATLAB function bwareaopen. Save the resulting segmentation masks in a .mat file in the MATLAB structure for that patient and visit.

[0107] Step (35) Approximate Drug Delivery. To approximate the concentration of drug delivered throughout the tumor tissue, we utilize physiological parameters derived from perfusion / diffusion analysis and each patient's individual therapy regimen. We use equation (9) (Kety-Tofts model), which approximates the concentration of drug delivered throughout the tumor tissue using the parameters (K trans , v e , and v p ) is available for each voxel (from the analysis described in step 24), and C p The (t) term is replaced with the measured plasma drug concentration from the population curve for each drug the patient received, administered repeatedly according to each patient's specific therapy regimen. The resulting drug distribution is saved as a 4D matrix (time is the fourth dimension) in a .mat file within the MATLAB structure for that patient and visit.

[0108] Step (36) Calculate summary measures of the tumor. Calculate the total tumor cellularity by summing the tumor cell counts across all voxels within the tumor ROI. The total number of voxels within the segmented tumor ROI and the voxel volume (8.45 mm using the DCE-MRI spatial resolution described in section 7 above) are then calculated. 3 ) to approximate tumor volume. To calculate the longest axis of each tumor, the 3D tumor ROI is evaluated using the function regionprops3 in MATLAB, which approximates the longest axis possible within the 3D object. These measures of cellularity, volume, and longest axis are applied to all model predictions to allow direct comparison with clinically measured data. Note that by performing this process of image segmentation and calculation of the longest diameter, the evaluation is as rigorous as possible.

[0109] The following steps provide a method for implementing a mathematical model to utilize the patient-specific MRI derived quantities described above to generate an individual patient response prediction.

[0110] Tumor prediction,timing is approximately 10 hours per patient and may be up to,a few days depending on the size of the tumor.

[0111] The step (37) model (i.e., equations (1-3)) is implemented in 3D. A fully explicit finite difference scheme with central differencing in space and forward differencing in time is used. Voxel dimensions are defined by the size of the DCE-MRI voxel grid, with a maximum diffusion coefficient of 0.001 mm. 2 / day, with a time step of Δt = 0.25 to ensure numerical stability. We set the size of the computational domain in square form, whose dimensions are determined by the size of the breast domain for each patient. For mechanical coupling, we assign the tissue stiffnesses in Table 1 to the tumor, fibroglandular and adipose tissue ROIs. We apply a no-flux boundary condition to the tumor cells and set the tissue displacements on the boundaries of the breast domain to zero (i.e., x on the boundaries).

number

[0112] Step (38) Calibrate the model parameters for each patient. Two MRI data sets for each patient are used to calibrate the mathematical model and simulate the model at the time of subsequent scans or surgery to predict tumor response. For example, the data sets from visits 1 and 2 are used to calibrate a first therapy regimen, which allows simulating the model and predicting the tumor response measured at visit 3. Similarly, the data sets from visits 3 and 4 are used to calibrate a second therapy regimen and simulate the model and predict the tumor response (determined by pathology) at the time of surgery. See FIG. 2 for an illustration of this calibration and prediction strategy.

[0113] The cellularity maps (derived from ADC of DW-MRI data - see steps 20 and 33) of two scans of each patient are used to calibrate the model parameters. Here, the tumor from the previous imaging visit initializes the calibration (together with tissue maps corresponding to mechanical binding and drug distribution) and the later tumor ROI is the target of the calibration. The D0 and α parameters are globally calibrated and the k parameter map is spatially calibrated, where the remaining parameters are assigned to the literature values ​​in Table 1 and θ is calculated directly (step 33). The calibration uses Levenberg-Marquardt least squares, a nonlinear optimization. Here, the sum of squared error between the tumor cell count simulated from the model and the tumor cell count calculated from the imaging data is calculated using the following parameters: number of iterations = 200, initial lambda = 10 -20 , lambda increment factors = 9 and 11, and the Jacobian calculation is performed every 25 iterations. Additionally, the global parameters D0 and α are constrained to be greater than zero, with D0 < 0.001 mm 2 Stopping criteria are a sum of squared errors <0.001, a match correlation coefficient between the model simulation and the target distribution of tumor cells of 1, and / or after a maximum number of iterations has been achieved.

[0114] Step (39) Evaluate Uncertainty. The following predictions and evaluations may take into account uncertainties associated with the calibrated parameters and data measurements. Choose three representative data sets from the cohort. Simulate the model from scan image 1 to scan image 2 (the calibration target scan image) for each tumor using the calibrated parameters. Add an appropriate range of noise (e.g., as determined by a repeatability study of the DW-MRI data) to voxels in each tumor cell map of the model-generated scan image 2 results using a normal distribution (MATLAB function randn). Calibrate the model to the tumor cell map of the noisy scan image 2 and save the resulting parameter values. Repeat for a total of 100 sets (total N=300) for each tumor in the three representative data sets. Calculate the percent difference for each parameter between the corresponding original parameter value and each result from fitting the noisy data. Calculate the 95% confidence interval for each parameter using all samples. Determine whether a uniform or normal distribution represents the resulting parameter difference distribution.

[0115] For each patient in the cohort, generate 50 random parameter sets by sampling the distribution of parameters defined by the calculated 95% confidence intervals, centered around each patient's calibrated parameter set (using MATLAB's rand for uniform distribution and randn for normally distributed values). For each of the 50 randomized parameter sets for each patient, simulate the model from scan image 2 to scan image 3. For each patient, calculate 95% confidence intervals across all 50 resulting tumor predictions for total cellularity, volume, and longest axis. Thus, the simulation results will include confidence intervals based on the uncertainty of the parameter estimates.

[0116] Step (40) Predict tumor response. Using the corresponding calibrated parameters and maps (tissue, tumor cellularity, drug distribution) from scan image 2, simulate a model from scan image 2 to scan image 3. Using summary measures (total tumor cellularity, total tumor volume, and longest axis, step 36), evaluate the 3D prediction of the outcome using the measured tumor response. Additionally, compare the predicted tumor and the measured tumor in scan image 3 with Dice coefficient (which measures the overlap between the predicted and measured ROIs; a Dice of zero indicates no overlap and a Dice of one indicates complete overlap) and / or concordance correlation coefficient to directly compare the prediction and the measurement for each patient. The predicted percent change from baseline to scan image 3 can be compared to the actual response as defined by RECIST.

[0117] Simulate the model from scan image 4 to scan surgery using the corresponding calibrated parameters and maps (tissue, cellularity, drug distribution) from scan image 4. Compare the resulting summary measures and their corresponding expected percent change to the surgery-defined response of the cohort (i.e., pCR group vs. non-pCR group) using statistical comparison tests appropriate for the cohort size (i.e., t-tests, Wilcoxon rank sum tests, Kendall correlation coefficient, and Pearson correlation coefficient, etc.). See Table 3 for troubleshooting guidance.

[0118] Timing: Here we provide a detailed breakdown of the timing of each individual step within the protocol. All times are approximate.

[0119] Step 1 - Varies widely depending on local patient acquisition method. Step 2 - 15 minutes for patient consent. Step 3 - 10 minutes for placing IV lines. Step 4 - 2 minutes for determining FOV of MRI exam. Step 5 - 1 minute 39 seconds for acquiring DW-MRI data. Step 6 - 34 seconds for acquiring B1 map data. Step 7 - 3 minutes 13 seconds for pre-contrast high-resolution T1 weighted scan and 99 seconds for variable flip angle T1 weighted images. Step 8 - 8 minutes for acquiring DCE-MRI data. Step 9 - 3 minutes 13 seconds for pre-contrast high-resolution T1 weighted scan. Step 10 - 5-10 minutes for uploading to PACS. Step 11 - 6 minutes for reading and converting DICOM. Step 12 - 1 minute for confirming slice position using DICOM header information. Step 13 - 3-4 minutes for upsampling DW-MRI and B1 map data to the resolution of DCE-MRI data. Step 14 - 1 minute for aligning DW-MRI data to DCE-MRI data. Step 15 - 7 min to align the variable flip angle T1 weighted images with the DCE-MRI data. Step 16 - 1 min to align the B1 data with the DCE-MRI data. Step 17 - 9 min to correct for motion in the DCE-MRI data. Step 18 - 20 min to segment the tumor. Step 19 - 5 min to perform quality control of the acquired data. Step 20 - 4 min to calculate the ADC map. Step 21 - 22 min to calculate the B1 corrected T1 map. Step 22 - 5 min to construct the AIF for the patient. Step 23 - 5 min to determine the BAT for each voxel in the DCE-MRI dataset. Step 24 - 25 min to perform pharmacokinetic analysis for each voxel in the DCE-MRI data. Step 25 - 2 min to enhance the anatomical images. Step 26 - 10 min to compare slices across all visits for each patient. Step 27 - 1 min to convert the file type to .mhd. Steps 28-30 - 1 day to identify the longitudinal alignment weights for all patient data from all patient visits. Step 31 - 15 minutes to align parameter maps from all visits for each patient.Step 32 - 5-10 min to select the domain for mathematical modeling for each patient set (i.e. all scan sessions from one patient). Step 33 - 1 min to convert the ADC maps into cell count estimates. Step 34 - 1 min to segment breast tissue into fibroglandular and adipose tissue. Step 35 - 1 min to estimate spatiotemporal distribution of drug concentrations. Step 36 - 1 min to calculate estimates of tumor cells, total tumor volume and longest axis. Step 37 - 5 min to simulate the model along spatial and temporal directions. Step 38 - Several hours to 2 days depending on the size of the domain to calibrate the model to the patient data. Step 39 - 1 day to evaluate the uncertainty of the model predictions. Step 40 - 5 min to evaluate the predictive ability.

[0120] result The following sections describe and present exemplary results associated with processing one dataset throughout the protocol. The exemplary dataset used is that of a breast cancer patient with triple-negative (estrogen receptor, progesterone receptor, and human epidermal growth factor receptor 2 negative) invasive ductal carcinoma in the left breast. At the time of the first imaging session, the patient was 59 years old with a body mass index of 28.3. All four scans were acquired over a six-month period during which the patient was treated with NAT. By surgery, the patient was determined to have residual disease (i.e., non-pCR outcome) at the end of NAT.

[0121] Image Acquisition: See FIG. 6 for an example of the resulting MRI data collected on the patient's first scan from the data acquisition step. Panel a shows a diffusion weighted image with a low b value (see step 5), and panel b shows a B1 map that quantifies the difference between the prescribed flip angle and the flip angle each voxel actually experiences (step 6). Panel c represents a single 10° flip angle image from the multi-flip angle data acquired to estimate the T1 map (step 7), and panel d represents the average of all images acquired during the DCE-MRI sequence (step 8). Note the differences between the various MRI modalities, and in particular observe how the tumor and other tissue structures are visualized by the DCE-MRI data (panel d).

[0122] Data Processing and Analysis: See FIG. 7 for how the exemplary images of FIG. 6 are processed using the protocol (steps 10-25) to identify various attributes of a patient's tumor. In particular, an ADC map is derived from the DW-MRI data acquisition (step 20). A tumor ROI is defined using the DCE-MRI data and the FCM algorithm (step 18), a corrected T1 map is generated from the variable flip angle T1 images and the B1 map (step 21), and the resulting Kety-Tofts parameters characterizing the vascular structure within the tumor are derived utilizing the DCE-MRI data, the tumor ROI, and the corrected T1 map (step 23).

[0123] Mapping the imaging data to a mathematical model: see FIG. 8 for exemplary intermediate and subsequent images of the step of converting the quantitative data map into quantities utilized within the mathematical modeling system. Here, for example, images of all four scans can be visualized within the inter-visit registration panel (steps 26-31). After inter-visit registration, a modeling domain is identified across the breast (step 32), ADC values ​​within the tumor ROI are converted to tumor cellularity (step 33), a mask of fibroglandular and adipose tissue is generated using a k-means clustering algorithm (step 34), and the drug distribution within the patient's tissues is approximated using the Kety-Tofts model and the plasma concentration curve of the patient's therapy regimen (step 35). For this patient, the resulting summary measures described in step 36 for visits 1-4 are as follows: Total cellularity (cell count) is 1.53×10 9 , 1.42×10 9 , 1.13×10 9 , and 1.20×10 9 Volume (mm 3 ) are 14,462, 11,831, 5,813, and 4,551. The longest axis (mm) is 36, 34, 32, and 29, respectively.

[0124] Tumor Prediction: See FIG. 9 for an exemplary comparison of predicted values ​​from the mathematical model with experimentally measured data for the three central slices at the time point of scan image 3 for the same patient data shown in FIGS. 6-8. Specifically, the mathematical model was calibrated using the patient's data from the first two scan images, and then the model was simulated in time from the time point of scan image 2 to scan image 3 with the resulting patient-specific parameters (as described in steps 37-40). Note that the model is able to predict the patient's tumor response with less than 17% error for three summary measures of the tumor (total cellularity, volume, longest axis) and has strong statistical correlation with tumor shape and cellularity (see figure caption for details). [Table 1] [Table 2] [Table 3-1] [Table 3-2] [Table 3-3] [Table 3-4]

[0125] A model constrained by MRI data to evaluate patient-specific neoadjuvant regimens One approach to individualizing the mathematical model is to leverage physiological information from non-invasive imaging data (obtained in 3D and at multiple time points) to initialize and constrain model parameters, thereby enabling patient-specific predictions. Biologically derived mathematical models based on image information can provide accurate predictions of tumor development in kidneys, brain, lungs, pancreas, etc. Various embodiments use DW-MRI data to estimate tumor cellularity, mechanistic coupling between breast tissue characteristics and tumor growth based on each individual patient's anatomy, and estimates of drug delivery to each voxel via DCE-MRI. Various embodiments incorporate mathematical descriptions of chemotherapy attenuation and efficacy (calibrated from each patient's data set). Various embodiments allow for the identification of anticipated alternative personalized dosing strategies that are superior to the therapy regimen each patient would receive as standard of care. In various embodiments, the necessary imaging data can be acquired at a community-based radiology center (i.e., not a research-oriented academic medical center) using widely available hardware. Since the various embodiments can be used in such facilities where the majority (85%) of oncology patients receive care, the disclosed approaches dramatically increase the population that these technologies may serve in the future.

[0126] Study data and methods Patient population: Quantitative MRI data were acquired in a cohort of 21 patients diagnosed with intermediate- to high-grade invasive breast cancer who were eligible for NAT as a component of their clinical care. Participants were cared for and imaged in community care settings (i.e., not academic research-oriented medical centers). Table 4 summarizes the main clinical characteristics of the patient population. [Table 4]

[0127] For the majority of patients, NAT consisted of two regimens. For example, patient 1 received doxorubicin and cyclophosphamide (regimen 1) followed by paclitaxel (regimen 2). MRI data were acquired four times throughout NAT: 1) pre-therapy, 2) after one cycle of the first therapy regimen, 3) at the completion of the first therapy regimen, and 4) after one cycle of the second therapy regimen. In this study, we utilize the first three data sets that summarize the first regimen. The NAT regimen included cycles of doxorubicin and cyclophosphamide administered approximately every 2 weeks for 4 cycles, paclitaxel (with or without carboplatin and / or targeted therapy) administered every 3 weeks for 4 cycles (total of 12 doses), and docetaxel (with carboplatin and targeted therapy) administered every 3 weeks for 6 cycles. As is common in standard care settings, there were some variations in the regimens prescribed by the treating physicians (variations are reported in Table 4). Note that 14 patients received only cytotoxic therapy in their initial NAT regimen, which we refer to as the "chemotherapy" group, whereas 7 patients received additional therapy (targeted therapy or immunotherapy), which we refer to as the "chemotherapy+" group.

[0128] MRI Data Acquisition: MRI was performed at two regional imaging facilities. MRI technologists at each facility were directly involved and responsible for patient positioning and implementation of the study imaging protocol. Thus, the image acquisition protocol was designed to be practical for routine imaging using widely available hardware and expertise.

[0129] Five MRI data types were acquired in each scanning session: 1) pre-contrast T1 map, 2) pre-contrast B1 field map to correct for radio frequency (RF) inhomogeneity, 3) DW-MRI data, 4) high-temporal resolution T1-weighted DCE-MRI data before, during, and after injection of gadolinium-based contrast agents (Gadovist, Bayer, Ontario, Canada, or Multihance, Bracco Diagnostics, New Jersey, USA), and 5) high-resolution T1-weighted anatomical scan images (post-contrast). The Supplementary Material provides a more detailed summary of the acquisition parameters for each of these measurements.

[0130] Data Analysis: Here, we briefly summarize various embodiments of the data processing method (additional details are provided in the Supplementary Materials section below). The first step involves intra-scan registration of MRI data within each scan session to correct for motion via rigid body registration. Next, tumor regions of interest (ROIs) are identified based on post-contrast scans, and estimates of tissue attributes related to vasculature perfusion permeability are quantified by analyzing the DCE-MRI data using a standard Kety-Tofts model. The DW-MRI data is analyzed to return a map of apparent water diffusion coefficients (ADCs). The third step is inter-scan registration to align images and calculated maps across all of a patient's imaging sessions into a common domain. The final step involves the calculation of certain quantities used within the mathematical model. These quantities include approximating the number of tumor cells from voxel-based ADC values, and segmenting fibroglandular and adipose tissues (based on the augmentation of the DCE-MRI data). To approximate the drug distribution within each voxel of tissue, a normalized map of blood volume is calculated and then scaled by the peak concentration of the drug (as estimated from the Kety-Tofts model; see Supplementary Materials below) to define the initial drug distribution across the entire domain at each administration time point of therapy.

[0131] Tumor volumes and longest axis from each patient's scans are automatically calculated on the inter-scan registration images, and Response Evaluation Criteria in Solid Tumors (RECIST) is used to assign each patient's response. The inventors note that this process of image segmentation and longest diameter calculation is performed to make the RECIST assessment as rigorous as possible and not systematically penalize in comparison with mathematical modeling. The inventors note that when calculating the longest axis of the tumor, only the central bulk tumor is considered, ignoring the smaller cut sections. Using the RECIST assessment, the tumor volumes of scans 1 and 3 were compared, and responders were patients with complete or partial responses (CR and PR, respectively), and non-responders were patients with stable disease or progressive disease (SD and PD, respectively).

[0132] Model: In various embodiments, a 3D mathematical model that includes the mechanistic coupling of tissue properties to tumor growth and therapy delivery may be used. The model may be initialized with patient-specific quantitative MRI data to predict therapy response. This approach may be extended to consider multiple chemotherapy regimens. The governing equation for the spatiotemporal evolution of tumor cells (NTC) is:

number

number

number

[0133] The therapy term describes the spatiotemporal distribution of each drug within the tissue and its effect on the cells in each voxel. Here, we extend the model to recognize different potencies and decay rates using the following equation:

number

number

number

[0134] Model Parameterization and Predictive Ability Assessment: An embodiment of the disclosed approach uses the first two MRI data sets of each patient to calibrate the mathematical model, and the third MRI data set to assess the mathematical model's ability to predict tumor response to therapy. Specifically, various embodiments use the cellularity maps of each tumor before and after the first cycle of NAT (scans 1 and 2, respectively) to calibrate the model parameters (D0, α, and β parameters are global, and k is spatially determined). Using these calibrated patient-specific parameters, the model is reinitialized with the tumor cellularity, tissue, and drug distribution maps of scan 2 and run up to the time of scan 3 to predict tumor response. See FIG. 1 for a graphical depiction of this strategy and the supplementary material below for additional details of the numerical implementation.

[0135] The predictive ability of the model is evaluated by comparing three measures that quantify tumor response: total tumor cellularity, tumor volume, and longest axis of the tumor. Three evaluations are performed to evaluate the predictive ability of the model using these measures. First, the error between the model's predicted values ​​and the patient's actual tumor values ​​(determined from the corresponding third scan image) is calculated for the three measures, and the significance of the model's predictive accuracy is also calculated. Second, total cellularity, volume, and longest axis are evaluated in the entire cohort to determine the level of agreement between the model's predicted values ​​and the measurements from the data as a group. Third, the predicted percent change in the three tumor measures from scan image 1 to scan image 3 is compared between the two RECIST-defined response groups.

[0136] Simulation of modified therapy regimens for individual patients: Alternative regimens are proposed for the group of patients (chemotherapy group) with the highest correlation between the model predictions and the observed values. The proposed alternative regimens use the same total doses that each patient received during the NAT regimen from scan 2 to scan 3, but the alternative regimens differ in individual doses and frequency between their second and third scans. More specifically, the embodiment simulates the effect of alternative regimens consisting of doses that are 1 / 2, 1 / 3, or 1 / 4 of the doses that the patient received, but administered twice, three times, or four times as frequently, respectively, as the patient received. Daily dose ratios were also investigated. For example, for a patient who received chemotherapy every 2 weeks, the alternative regimens would be 1 / 2 dose every week, 1 / 3 dose every 4-5 days, 1 / 4 dose every 3-4 days, and 1 / 14 dose every day.

[0137] An embodiment runs the model from scan 2 to scan 3 time points using each patient's pre-calibrated parameters and alternative regimens administered to assess differences in tumor response across all regimens (see FIG. 1). The model simulation serves as an in silico twin to determine each patient's therapeutic response to alternative regimens. This allows for the identification of alternative individualized therapeutic regimens that are hypothesized to lead to better tumor control compared to the current standard "one size fits all" approach.

[0138] Analysis: To summarize the predictive performance of the model, the absolute percent difference between the tumor response measured from the third scan data and the corresponding predicted value from the model is calculated for each patient (median and interquartile range reported). The accuracy of patient-specific predictions for total cellularity, volume, and longest axis is tested by generating Monte Carlo estimated p-values ​​(see supplementary materials) to determine the significance of 10, 15, or 20% absolute differences between the model predicted value and the measured outcome. (The range of 10-20% corresponds to the accuracy of the measurement.) Due to the modest sample size, it is difficult to accurately assess the normality of the data. Therefore, we calculate the Kendall correlation coefficient (KCC) to determine the agreement between the model predicted tumor response and the observed value at the scan image 3 time point across the entire cohort. However, to allow comparisons between previous and future efforts with larger datasets, we also report the concordance correlation coefficient and Pearson correlation coefficient (CCC and PCC, respectively). A two-tailed Wilcoxon rank sum test was used to determine significant differences in median percent change from scan 1 to scan 3 between responding and non-responding groups, with p<0.05 considered significant.

[0139] To determine which therapy regimen provides the greatest tumor control for each patient, embodiments calculate the percent change from the start of the alternative regimen (scan 2) to the expected tumor at the time of scan 3. These percent changes are compared across each patient's regimens, and the inventors determine the "optimal" regimen based on the greatest tumor shrinkage / control. Then, embodiments calculate the difference in percent change between the standard regimen and the optimal regimen to determine the additional percent tumor shrinkage potentially achieved if the alternative regimen is administered. The median and interquartile range of these values ​​are reported. Using a paired two-sided Wilcoxon signed rank test, the inventors determine whether there is a significant difference between the tumor control achieved by the standard regimen group compared to the alternative regimen group by comparing the percent reduction from scan 2 to scan 3 and the expected total tumor burden at scan 3 between the two groups.

[0140] Results: Three patients are excluded from the analysis. For patient 1, the tumor invaded the chest wall, violating the flux-free boundary condition assumption adopted in the numerical implementation (see Supplementary Material). For patient 7, a silicone breast implant caused significant artifacts in the DW-MRI data. For patient 16, a third scan was not collected due to a scheduling conflict. Exclusion of these three datasets reduces the total cohort to N=18 and the chemotherapy subgroup to N=13.

[0141] FIG. 11 shows a comparison of predicted and experimentally measured cellularity maps for a representative patient when the calibration model was run up to scan 3 using information from scans 1 and 2. The number of tumor cells predicted by the model is overlaid (in color) on a grayscale anatomical image of the breast. While areas of higher and lower cellularity may not directly match between predicted and observed values, the model is able to capture the general shape of individual tumors, and there is a small error between predicted and measured values ​​for the patient scans (listed in the figure caption). Table 5 summarizes the results of the absolute percent error for all three tumor measures (i.e., total cellularity, total volume, and longest axis) compared to measurements from the third scan for each patient. FIG. 12 depicts the distribution of sample means generated by the Monte Carlo method to determine the significance of the model's predictive accuracy for 10, 15, and 20% absolute differences between the model's predictions and the measured outcomes. For the cohort (N=18), model predictions for the longest axis are significant (p<0.05) for all three thresholds. Measures of total cellularity and volume tend to be significant (i.e., p<0.1) at the 10% and 15% thresholds, achieving p-values ​​<0.05 at the 20% threshold. [Table 5]

[0142] Figure 13 is a scatter plot comparing predicted and actual tumor response. For all three measures, the model's predictions are strongly correlated with the actual tumor response. In particular, we found that across the entire cohort, cellularity CCC / PCC=0.91 / 0.92, volume CCC / PCC=0.90 / 0.90, and longest axis CCC / PCC=0.86 / 0.88 (p<0.01). Considering the KCC measures, for these tumor measures, the predictions are correlated with the actual tumor response. For total cellularity, total volume, and longest axis, KCC=0.59, 0.65, 0.76 (p<0.01), respectively. Note that for all three tumor measures, the correlation is higher in the chemotherapy subgroup compared to the chemotherapy+ group. In particular, we found that the chemotherapy subgroup had greater CC values ​​than the entire cohort, with CCC / PCC=0.92 / 0.94 for cellularity, CCC / PCC=0.90 / 0.91 for volume, and CCC / PCC=0.92 / 0.96 for longest axis (p<0.01). For KCC measures, the predicted values ​​are strongly correlated with actual tumor response for these tumor measures: KCC=0.72, 0.77, 0.85 for total cellularity, total volume, and longest axis, respectively (p<0.01).

[0143] The model predictions are also compared to the tumor response status determined by RECIST at the end of the NAT regimen. Application of RECIST to each tumor from scan 1 to scan 3 resulted in eight patients being labeled as responders and ten patients as non-responders (see Table 1). See Supplementary Table S.4 for the percent change in the longest axis (as well as the percent change in total cellularity and volume) from scan 1 to scan 2 and 3 for each patient. When evaluating the observed percent change from scan 1 to scan 3, the percent change in each of the three tumor measures resulted in significantly different medians between the responder and non-responder groups for total cellularity (p<0.05), total volume (p<0.04), and longest axis (p<0.001). In contrast, the model predicted a significant median percent change between responders and non-responders for the longest axis (p<0.002), while for total cellularity and total volume, p=0.08. See Table 3 for median and interquartile range of measured and expected percent change in all three measures. Note that when comparing responder and non-responder values ​​in the chemotherapy subgroup (N=6 and N=7, respectively), the model predicts significantly different median changes for all three tumor measures (p<0.04 for total cellularity and volume, p<0.01 for longest axis). The model-predicted percent change for this subgroup is also listed in Table 6. [Table 6]

[0144] As the model predictions show the greatest correlation with chemotherapy subgroups, an in silico study of alternative regimens is performed for this subset of patients (see Table 1, N=13). This choice is discussed further below. Using each patient's parameters calibrated from scans 1 and 2, the model is simulated from scan 2 to scan 3 to determine whether greater or less tumor burden control / reduction would be achieved with the proposed hypothetical alternative treatment schedule (i.e., 1 / 2, 1 / 3, 1 / 4, and daily dose fractions given daily at 2x, 3x, and 4x the frequency of the standard dose, respectively). Comparing the percent change from scan 2 to scan 3 for two measures of tumor size (volume and longest axis) between standard and alternative regimens, the median total volume differed by 8% and the median longest axis by less than 1%. However, the difference between standard and alternative regimens that achieved the greatest reduction in total cellularity per individual patient ranged from 0% to 46%, with a median difference of 17% (interquartile range [6%, 35%]). Positive differences represent additional reduction / percent control that could have been achieved if the patient had received the alternative regimen. The 0% example is one patient for whom the standard regimen was optimal among all patients tested. However, all other patients had further reductions in total cellularity when comparing the standard regimen to the best alternative regimen. Figure 14 summarizes the predicted percent change from the start of the alternative regimen (scan 2) to the completion of therapy (scan 3) and a heat map showing the difference between the alternative regimen compared to the standard care regimen received by each patient. Compared to the standard regimen, the optimal regimen predicted by the model resulted in a significant reduction in total tumor cell count from scan 2 to scan 3 (p<0.001) as well as a significant reduction in the predicted total tumor cell count at scan 3 (p<0.001). Note that no single regimen was most effective across all patients: the daily dosing regimen was best for 4 patients, while the quarter, third, half, and standard dosing regimens were best for 2, 2, 4, and 1 patient, respectively.Figure 2 shows an example of the changes in total cellularity resulting from two alternative regimens for one patient.

[0145] considerations In an exemplary embodiment, a mathematical model describing patient-specific growth and spread of tumor cells, mechanical properties of breast tissue, and treatment regimen was individually calibrated using quantitative MRI data acquired in a community-based care environment. After patient-specific calibration, the model was run forward in time to predict tumor response upon completion of the first therapy regimen. The significance of the model's predictive accuracy was established by analyzing the various absolute differences between the model predictions and the measured outcomes (Figure 12). The model predictions were also found to be significantly correlated with the actual tumor outcomes of the cohort for three different measures of tumor response (total cellularity, volume, and longest axis; Figure 13). Notably, the CCC values ​​for the entire cohort were > 0.86 when all three tumor measurements were compared. Furthermore, the model predictions were significantly different for the change in the longest axis when compared between RECIST-defined responding and non-responding groups (Table 6). It is important to note that RECIST is only an assessment of response and is not intended to be used to predict response. In fact, the RECIST designations (i.e., CR, PR, SD, PD) for scans 1 and 2 changed in 8 of 18 patients when compared to the RECIST designations from scans 1 and 3. However, recall that the mathematical model only required data from scans 1 and 2 (after only one cycle of therapy) to predict the response observed at the end of the initial NAT regimen.

[0146] The results show that various embodiments of the model are superior in predicting tumor response to chemotherapy-only regimens. For the chemotherapy-only subgroup, the predicted values ​​of the mathematical model were more highly correlated with the measured tumor response of each patient compared to the chemotherapy+ subgroup, and the model had significantly different predicted values ​​between responder and non-responder patients for all three tumor measures. In some embodiments, the model does not explicitly consider the effects of targeted therapy and / or immunotherapy, but only implicitly through the calibration of parameters (such as growth maps). Thus, in these embodiments, the chemotherapy subgroup was selected for in silico alternative regimen investigation. In various embodiments, modeling efforts can be made to consider the effects of targeted therapy to increase the generalizability of the methodology to all breast cancer subtypes.

[0147] Using parameters derived from each individual patient in the chemotherapy subgroup, several alternative treatment regimens were evaluated that differed in frequency and dosage of the drug from the standard of care regimen received by each patient (but the total amount of drug remained constant). Thus, the mathematical model served as an in silico twin for each patient, and an in silico clinical trial was performed in four branches, each branch testing alternative schedules of total dosage and time course of the standard of care. The overall tumor control achieved by the group of standard regimens was found to be statistically inferior to the group of the optimal regimen selected for each patient. Since no single regimen is optimal for all patients, the success of any particular treatment approach depends on the tumor characteristics of the individual patient and the importance of identifying a patient-specific treatment protocol to maximize or otherwise enhance the response. Various embodiments of the disclosed patient-specific mathematical modeling systematically exploit this realization to individually optimize the response. Furthermore, it is noted that other possible alternative regimens can be used to identify a patient-specific therapy regimen that maximizes efficacy.

[0148] The aim of this study was to demonstrate the use of biophysical and mathematical methods to optimize therapy on a patient-specific basis. The above analysis presented showed that the predictions of the mathematical model were significantly correlated with and predicted the actual response. It is also important to note that the results of this study would be very difficult to achieve using artificial intelligence methods, which typically require huge training databases to identify patterns that emerge at a population level (rather than at an individual level).

[0149] In various embodiments, the number of data points (i.e., scan time) used to calibrate the model can be increased to enhance the predictive capabilities of the approach. In particular, if the change in the tumor between calibration scans is minimal compared to the overall change observed at the end of therapy, the model's ability to capture the global dynamics may be compromised. Acquiring additional data before and during therapy will allow the model calibration scheme to more accurately determine patient-specific parameter values.

[0150] Various embodiments use biologically-based mathematical models in predicting tumor response using data obtained from individual patients, starting from the earliest time points during neoadjuvant regimens. In silico results show how therapeutic regimens can be tailored and even optimized for each patient using mathematical models and simulation studies. These results individualize patient regimens through quantitative imaging and mathematical modeling.

[0151] Supplementary Material MRI data acquisition: The two imaging facilities are an outpatient imaging facility and a regional hospital that routinely perform breast MRI, but with different service agreements and quality control guidelines. It is again noted that the repeatability and reproducibility of quantitative MRI at these centers has been previously established. For breast MRI, a Siemens 3T scanner (Erlangen, Germany) equipped with an 8- or 16-channel receive dual breast coil (Sentinelle, Invivo, Gainesville, FL) was used. All images were acquired in the sagittal plane. Diffusion-weighted MRI (DW-MRI) was acquired using a monopolar, single-shot spin-echo, echo-planar imaging sequence with a diagonal diffusion encoding orientation. The six acquisitions were performed at 0 and 200 s / mm 2 The 18 acquisitions were averaged at a b-value of 800 s / mm 2 The images were averaged over b-values ​​of 100x1000, ...

[0152] The MRI protocol also included high-resolution T1-weighted, 3D gradient echo, FLASH (fast low-angle shot) acquisition with the following parameters: TRITE=5.3 / 2.3 ms, flip angle=10°, acquisition matrix 256 × 256, slice thickness=1 mm, GRAPPA acceleration 2, and SPAIR (spectrally selective adiabatic inversion recovery) fat suppression. The acquisition time for anatomical images was 3 min 11 s. The DCE-MRI protocol consisted of a T1-weighted, VIBE (volume-interpolated breath-hold examination; however, breath-holding was not used in these studies) acquisition with TRITE=7.02 / 4.6 ms, flip angle 6°, matrix=192 × 192, 10 slices of 5 mm thickness, and GRAPPA acceleration factor 3, resulting in a temporal resolution of 7.27 s, 1 min before and 6 min after administration of gadolinium-based contrast agent (Multihance (Bracco, Monroe Township, NJ) or Gadovist (Bayer, Leverkusen, Germany)) via a subsequent saline flush using a power injector. The B1 field was mapped using a Siemens TurboFLASH sequence and corrected for transmission inhomogeneity with the following acquisition parameters: TRITE=8680 / 2 ms, flip angle=8°, matrix=96 × 96, and slice thickness=5 mm. Because the B1 mapping protocol included a slice gap, two acquisitions were performed to cover the same field of view as the measurements above, resulting in a total acquisition time of 34 s. For a summary of all imaging parameters, see Table S.1. [Table 7]

[0153] Data Analysis Techniques: Figure 15 provides a summary of the data analysis and processing steps. For each data set, two types of image registration were performed: within scans (registration within a single visit) and between scans (registration across multiple visits). For each patient, the MR images acquired in each session (within scans) were aligned using a rigid registration algorithm where the B1, T1, and diffusion weighted images were registered to the DCE images. Images acquired at different resolutions compared to the DCE images were upsampled using a nearest neighbor approach via the function interp3 in MATLAB (MathWorks, Natick, MA). The rigid algorithm used for within scan registration was implemented using the function imregister in MATLAB. For each patient, all image data sets were registered to a common space across time (between scans) via a constrained non-rigid registration algorithm that preserved the tumor volume at each time point. Between scans registration uses an adaptive basis algorithm implemented using the software Elastix.

[0154] DCE-MRI data were used for both tissue segmentation and characterization of the vasculature of each patient. DCE-MRI data were used to segment tumor regions of interest (ROIs) at each time point using a fuzzy C-means based clustering algorithm. The clustering algorithm divides voxels into classes based on the probability weighting of voxels that are likely to belong to the tumor ROI. To generate masks of fibroglandular and adipose tissue, the MATLAB function adapthisteq was first applied to augment the DCE images after inter-scan registration using a contrast-limited adaptive histogram equalization algorithm. Fibroglandular and adipose tissues were segmented using a k-means clustering algorithm (these masks are used to assign tissue stiffness properties in the mathematical model detailfidbfilow).

[0155] DCE-MRI data were analyzed using the standard Kety-Tofts model.

number

number

number

number

number

[0156] To approximate the drug distribution within each voxel of tissue, a normalized map of blood volume is calculated by calculating the dynamic area under the curve (AUC) of the time course subtracted from baseline for each voxel, then normalizing by the maximum AUC value from the entire tumor ROI. This normalized blood volume map is then scaled by the peak concentration of the drug (estimated from the Kety-Tofts model) to define the initial drug distribution across the domain upon administration of each therapy. Specifically, the Kety-Tofts model (i.e., Eq. S1)) and physiological parameters derived from DCE-MRI per voxel are used to replace the concentration of contrast agent in plasma with the concentration of drug in plasma from the population curve measured for each therapy. Thus, the drug concentration within the tumor tissue is spatially heterogeneous and varies in time based on the individual patient's NAT schedule and vascular characteristics.

[0157] Apparent diffusion coefficient (ADC) values ​​for individual voxels were calculated using DW-MRI data via standard methods. The ADC values ​​for each voxel within the tumor (segmented using the methods described above) were then calculated for each 3D location.

number

number

number

[0158] Tumor volume was approximated as the product of the total number of voxels within the segmented tumor ROI and the DCE-MRI voxel volume. To calculate the longest axis of each tumor, the 3D tumor ROI was evaluated by MATLAB's regionprops3 function, which approximates the longest possible axis within a 3D object. These measures of tumor volume and longest axis were also applied to all model predictions (the models use the same domain as the MR images) for direct comparison with experimentally measured values.

[0159] Model: A 3D mathematical model including mechanical coupling of tissue properties to tumor growth and therapy delivery was designed, initialized with patient-specific quantitative MRI data of breast cancer to predict therapy response. This previous approach was extended to consider multiple chemotherapy terms. See Eq. (12) for the governing equation of the tumor spatiotemporal evolution. The first term on the right hand side describes the migration of tumor cells, the second term describes the logistic growth of cells, and the third term describes the effect of chemotherapy. (See Tables S.1 and S.2 for a description of the variables and parameters and how they are assigned.) For the growth terms,

number

number

number

[0160] The therapy term in equation (12) describes the spatiotemporal distribution of each drug in the tissue and its effect on the cells in each voxel. Here, we extend the model using equation (13) to recognize different efficacy and decay rates. (13). All simulation codes and numerical calculations were written and performed in MATLAB (MathWorks, Natick, MA). The model was implemented in three dimensions (3D) using a fully explicit finite difference scheme with Δt = 0.25 days, with mesh dimensions defined by the size of the DCE-MRI voxels. The size of the computational domain was set to be square, whose dimensions were determined by the size of each patient's breast domain. To reduce the computational time of all calibrations, the voxel matrix within this specified rectangular domain was downsampled by a factor of 2. A no-flux boundary condition was prescribed at the breast boundary. To calibrate the model parameters, Levenberg-Marquardt (LM) least-squares nonlinear optimization was used to minimize the sum of squared errors between the tumor cell count simulated from the model and the tumor cell density calculated from the imaging data. [Table 8] [Table 9]

[0161] Statistical Analysis: To test the accuracy of the model predictions, we generate Monte Carlo estimated p-values ​​to test the significance of the mean absolute difference between the model predictions and the measured outcomes versus random sampling. This is a standard bootstrap random resampling approach. In this particular application, the null hypothesis that the mean difference is not equal to a delta value of 10%, 15%, or 20% error (determined using established variability in MRI measurements) is tested via the following steps:

[0162] (1) For each patient, calculate the absolute difference between the model prediction and the tumor response measurement (N=18). (2) Calculate the cohort mean absolute difference. (3) Sample the absolute difference from step 1 500 times with replacement (sample size N=18). (4) Calculate the mean of each of the 500 samples created in step 3. (5) Subtract the corresponding mean delta from each of the 500 sample means. (6) Calculate the percentage of time that each delta-adjusted sample mean is less than the cohort mean calculated in step 2. (7) Subtract the percentage calculated in step 6 from 1.0 to determine a Monte Carlo estimated p-value.

[0163] These p-values ​​are calculated for each measure of tumor response (ie, total cellularity, volume, longest axis). [Table 10]

[0164] Various operations described herein can be implemented on computer systems having a variety of design features. Figure 16 shows a simplified block diagram of a representative server system 1600 (e.g., a computing device that analyzes MRI data) and a client computer system 1614 (e.g., a computing device that receives and presents imaging devices and sensors and / or analysis results) that can be used to implement certain embodiments of the present disclosure. In various embodiments, the server system 1600 or a similar system can implement the services or servers described herein, or portions thereof. The client computer system 1614 or a similar system can implement the clients described herein.

[0165] The server system 1600 may have a modular design incorporating multiple modules 1602 (e.g., blades in a blade server embodiment). Although two modules 1602 are shown, any number may be provided. Each module 1602 may include a processing unit 1604 and local storage 1606.

[0166] The processing unit 1604 may include a single processor that may have one or more cores, or multiple processors. In some embodiments, the processing unit 1604 may include a general-purpose primary processor and one or more special-purpose co-processors, such as a graphics processor or digital signal processor. In some embodiments, some or all of the processing unit 1604 may be implemented using customized circuitry, such as an application specific integrated circuit (ASIC) or a field programmable gate array (FPGA). In some embodiments, such integrated circuits execute instructions stored in the circuitry itself. In other embodiments, the processing unit 1604 may execute instructions stored in local storage 1606. Any type of processor may be included in the processing unit 1604 in any combination.

[0167] The local storage 1606 may include a volatile storage medium (such as, for example, a conventional DRAM, SRAM, or SDRAM) and / or a non-volatile storage medium (such as, for example, a magnetic or optical disk, or a flash memory). The storage medium incorporated in the local storage 1606 may be fixed, removable, or upgradeable, as desired. The local storage 1606 may be physically or logically divided into various sub-units, such as a system memory, a read-only memory (ROM), and a permanent storage device. The system memory may be a volatile read-write memory, such as a read-write memory device or a dynamic random access memory. The system memory may store some or all of the instructions and data that the processing unit 1604 needs at run time. The ROM may store static data and instructions needed by the processing unit 1604. The permanent storage device may be a non-volatile read-write memory device that can store instructions and data even when the module 1602 is powered off. As used herein, the term "storage medium" includes any medium capable of storing data indefinitely (subject to overwriting, electrical disturbances, or power loss, etc.), and does not include carrier waves and transitory electronic signals propagated over wireless or wired connections.

[0168] In some embodiments, local storage 1606 can store one or more software programs executed by processing unit 1604, such as an operating system and / or various server or computing functions, such as any of the components of FIGS. 1 and 12, or any other computing device, computing system, and / or sensor identified in this disclosure.

[0169] "Software" generally refers to a sequence of instructions that, when executed by the processing unit 1604, causes the server system 1600 (or portions thereof) to perform various operations, thus defining one or more specific machine embodiments that execute and perform the operations of the software programs. The instructions may be stored as firmware resident in a read-only memory and / or as program code stored in a non-volatile storage medium that can be loaded into a volatile working memory for execution by the processing unit 1604. The software may be implemented as a single program or as a collection of separate programs or program modules that interact as necessary. From the local storage 1606 (or non-local storage, described below), the processing unit 1604 may retrieve program instructions to execute and data to process in order to perform the various operations described above.

[0170] In some server systems 1600, multiple modules 1602 may be interconnected via a bus or other interconnect 1608 to form a local area network that supports communication between the modules 1602 and other components of the server system 1600. The interconnect 1608 may be implemented using a variety of technologies, including server racks, hubs, routers, etc.

[0171] A wide area network (WAN) interface 1610 may provide data communication capabilities between a local area network (interconnect 1608) and a larger network such as the Internet. Conventional or other active technologies may be used, including wired technologies (e.g., Ethernet, IEEE 802.3 standard) and / or wireless technologies (e.g., Wi-Fi, IEEE 802.11 standard).

[0172] In some embodiments, the local storage 1606 is intended to provide working memory for the processing units 1604, providing fast access to programs and / or data being processed while reducing traffic on the interconnect 1608. Storage for large amounts of data can be provided on a local area network by one or more mass storage subsystems 1612, which can be connected to the interconnect 1608. The mass storage subsystems 1612 can be based on magnetic, optical, semiconductor, or other data storage media. Direct attached storage, storage area networks, network attached storage, etc. can be used. Any data store or other collection of data described herein as generated, consumed, or maintained by a service or server can be stored in the mass storage subsystem 1612. In some embodiments, additional data storage resources can be accessible via the WAN interface 1610 (potentially increasing latency).

[0173] The server system 1600 can operate in response to requests received via the WAN interface 1610. For example, one of the modules 1602 can implement a monitoring function and assign separate tasks to the other modules 1602 in response to received requests. Conventional work allocation techniques can be used. Once the request is processed, the results are returned to the requester via the WAN interface 1610. Such operations can generally be automated. Furthermore, in some embodiments, the WAN interface 1610 can connect multiple server systems 1600 together, providing a scalable system capable of managing large volumes of activity. Conventional or other techniques for managing server systems and server farms (collections of cooperating server systems), including dynamic resource allocation and reallocation, can be used.

[0174] The server system 1600 can interact with a variety of user-owned or user-operated devices over a wide area network such as the Internet. One example of a user-operated device is shown in FIG. 16 as a client computing system 1614. The client computing system 1614 can be implemented as a consumer device such as, for example, a smartphone, other mobile phone, a tablet computer, a wearable computing device (e.g., smart watch, glasses), a desktop computer, a laptop computer, etc.

[0175] For example, client computing system 1614 may communicate over WAN interface 1610. Client computing system 1614 may include traditional computer components such as a processing unit 1616, a storage device 1618, a network interface 1620, user input devices 1622, and user output devices 1624. Client computing system 1614 may be computing devices implemented in a variety of form factors, such as a desktop computer, a laptop computer, a tablet computer, a smartphone, other mobile computing device, or a wearable computing device.

[0176] The processor 1616 and storage device 1618 may be similar to the processing unit 1604 and local storage 1606 described above. Suitable devices may be selected based on the demands placed on the client computing system 1614. For example, the client computing system 1614 may be implemented as a "thin" client with limited processing power or as a high performance computing device. The client computing system 1614 may include program code executable by the processing unit 1616 to enable various interactions with the server system 1600 of the message management service, such as accessing messages, performing actions on messages, and other interactions described above. Some client computing systems 1614 may also interact with the messaging service independent of the message management service.

[0177] Network interface 1620 may provide a connection to a wide area network (e.g., the Internet) to which WAN interface 1610 of server system 1600 is also connected. In various embodiments, network interface 1620 may include a wired interface (e.g., Ethernet) and / or a wireless interface implementing various RF data communication standards, such as Wi-Fi, Bluetooth, or cellular data network standards (e.g., 3G, 4G, LTE, 5G, etc.).

[0178] User input device 1622 may include any device (or devices) that allows a user to provide a signal to client computing system 1614. Client computing system 1614 may interpret the signal as indicating a particular user request or information. In various embodiments, user input device 1622 may include any or all of a keyboard, a touch pad, a touch screen, a mouse or other pointing device, a scroll wheel, a click wheel, a dial, a button, a switch, a keypad, a microphone, and the like.

[0179] The user output device 1624 may include any device that allows the client computing system 1614 to provide information to a user. For example, the user output device 1624 may include a display-to-display image generated by or delivered to the client computing system 1614. The display may incorporate a variety of image generating technologies, such as liquid crystal display (LCD), light emitting diode (LED) including organic light emitting diode (OLED), projection system, or cathode ray tube (CRT), along with supporting electronics (e.g., digital-to-analog or analog-to-digital converter, signal processor, etc.). Some embodiments may include devices such as a touch screen that function as both an input device and an output device. In some embodiments, other user output devices 1624 may be provided in addition to or instead of a display. Examples include indicator lights, speakers, haptic "display" devices, printers, and haptic devices (e.g., haptic sensory devices may vibrate at different rates or intensities at various times), etc.

[0180] Some embodiments include electronic components such as a microprocessor, storage, and memory that store computer program instructions on a computer-readable storage medium. Many of the functions described herein can be implemented as a process specified as a set of program instructions encoded on a computer-readable storage medium. When these program instructions are executed by one or more processing units, the program instructions cause the processing units to perform various operations indicated in the program instructions. Examples of program instructions or computer code include machine code, such as produced by a compiler, and files containing high-level code that are executed by a computer, electronic component, or microprocessor using an interpreter. Through suitable programming, the processing units 1604 and 1616 can provide various functions to the server system 1600 and the client computing system 1614, including any of the functions described herein as being performed by a server or client, or other functions associated with a message management service.

[0181] It will be understood that the server system 1600 and the client computing system 1614 are exemplary and that variations and modifications are possible. Computer systems used in connection with embodiments of the present disclosure may have other features not specifically described herein. Furthermore, while the server system 1600 and the client computing system 1614 are described with reference to certain blocks, it should be understood that these blocks are defined for convenience of description and are not intended to imply a particular physical arrangement of component parts. For example, different blocks may, but need not, be located in the same facility, in the same server rack, or on the same motherboard. Furthermore, the blocks need not correspond to physically separate components. The blocks may be configured to perform various operations, for example, by programming a processor or providing appropriate control circuitry, and the various blocks may or may not be reconfigurable depending on how the initial configuration is obtained. The embodiments of the present disclosure may be realized in a variety of apparatuses, including electronic devices implemented using any combination of circuitry and software.

[0182] In various embodiments, the code may be implemented using CPU-based programming logic. A multi-core CPU may be beneficial if some steps are to be performed in parallel, but this is not required (nor is multi-node necessary). For some datasets, all images from one visit may require, for example, approximately 3 gigabytes (GB) of memory; if each patient has three visits, approximately 10 GB may be allocated to maintain similar performance in terms of computation time for the particular embodiment discussed above.

[0183] Additional Example: MRI-Based Digital Models Predict Patient-Specific Treatment Response to Neoadjuvant Chemotherapy in Triple-Negative Breast Cancer Triple-negative breast cancer (TNBC) is persistently refractory to therapy, and methods to improve therapy targeting and response assessment in this disease are needed. According to potential embodiments, applicants integrate quantitative magnetic resonance imaging (MRI) data with biology-based mathematical modeling to accurately predict TNBC response to neoadjuvant systemic therapy (NAST) on an individualized basis. Specifically, 56 TNBC patients enrolled in the ARTEMIS trial (NCT02276443) received standard of care doxorubicin / cyclophosphamide (A / C) followed by paclitaxel for NAST, and dynamic contrast-enhanced MRI and diffusion-weighted MRI were acquired before treatment and after 2 and 4 cycles of A / C therapy. Biology-based models were established to characterize tumor cell migration, proliferation, and treatment-induced cell death. Two evaluation frameworks were investigated using: 1) Images acquired before and after two cycles of A / C for calibration and prediction of tumor status after A / C, and 2) Images acquired before, after, and after two and four cycles of A / C for calibration and prediction of response after NAST. In framework 1, the concordance correlation coefficients between patient-specific post-A / C changes in predicted and measured tumor cellularity and volume were 0.95 and 0.94, respectively. In framework 2, the biology-based model achieved an area under the receiver operating characteristic curve of 0.89 (sensitivity / specificity=0.72 / 0.95) for discriminating between pathological complete response (pCR) and non-pCR, which was statistically superior (P<0.05) to the value of 0.78 (sensitivity / specificity=0.72 / 0.79) achieved by tumor volume measured after four cycles of A / C. Overall, the model successfully captures the patient-specific spatiotemporal dynamics of TNBC response to NAST and provides highly accurate prediction of NAST response.

[0184] By integrating MRI data with biology-based mathematical modeling, we were able to successfully predict breast cancer response to chemotherapy, demonstrating that digital twins have the potential to facilitate a paradigm shift from simply assessing response to predicting and optimizing therapy efficacy.

[0185] Introduction: Neoadjuvant systemic therapy (NAST) is widely considered the standard of care for the treatment of stage II-III locally advanced triple-negative breast cancer (TNBC). NAST increases the success rate of breast-conserving surgery by reducing tumor burden and provides an opportunity to treat micrometastases in a naïve setting, thereby improving patients' progression-free survival. Importantly, TNBC patients who achieve a pathological complete response (pCR) in the neoadjuvant setting have good recurrence-free survival rates. In contrast, patients with residual disease after NAST are at higher risk of early recurrence and death. Unfortunately, based on a recent pooled analysis of 52 studies from 1999 to 2016, only 32.6% of TNBC patients treated with standard taxane / anthracycline-based NAST had pCR or minimal residual disease at the time of surgical resection.

[0186] The development of new neoadjuvant treatment regimens has provided an opportunity to tailor treatment to individual patients to improve TNBC outcomes. It is therefore increasingly important to develop techniques that can accurately predict the response to NAST in individual TNBC patients. If it can be conclusively determined that a therapy regimen is unlikely to achieve a patient's pCR, a risk-adapted therapy, with rationale-based addition of treatment and removal of unnecessary components, can be employed, potentially improving outcomes and reducing side effects. The importance of being able to remove patients from ineffective therapy as early as possible cannot be overemphasized, especially considering the significant toxicities, including increased chances of hospitalization, cardiac damage, leukemia, and even death. Furthermore, accurate prediction of response to NAST may allow the identification of exceptional responders who may benefit from de-escalation of treatment, including the possibility of non-surgical management of the disease.

[0187] Much effort has been devoted to investigating approaches that can accurately distinguish pCR and non-pCR patients at the early stage of NAST. Imaging biomarkers derived from magnetic resonance imaging (MRI), positron emission tomography, and ultrasound imaging have been shown to be strongly correlated with breast tumor response to NAST. In particular, MRI measurements before and during NAST are valuable predictors of pCR, especially functional tumor volume (FTV) derived from dynamic contrast-enhanced (DCE-) MRI and apparent diffusion coefficient (ADC) derived from diffusion-weighted (DW-) MRI. More recently, artificial intelligence methods have been used to extract features from high-dimensional data to build predictive models for distinguishing pCR from non-pCR in breast cancer. However, most of these approaches to predict or assess response have the inherent limitation of being population-based. Importantly, population-based approaches that rely solely on statistical inference from the characteristics of a large population inevitably obscure the specific status of individual patients over time, especially for a heterogeneous disease such as cancer. Conversely, biology-based models using patient-specific data could shift the paradigm from population-based to individual-based approaches. Furthermore, biology-based mathematical modeling of tumor response can not only predict changes in global indices summarizing tumor burden (e.g., total tumor volume and cellularity), but also reveal biologically specific information (e.g., spatially resolved maps of proliferation, pharmacokinetics, and individual patient sensitivity to administered therapy). This promises unique opportunities to characterize tumor pathophysiology, rigorously predict long-term outcomes, and even optimize treatment plans on a patient-specific basis.

[0188] As discussed herein, in accordance with various potential embodiments, the applicants have developed a clinical computational approach to establish a patient-specific model for early prediction of spatiotemporal development and response of individual TNBC patients to standard-of-care NAST. Because the model can represent a physical object (i.e., tumor), predict the behavior of the object given an influence (i.e., treatment), and enable decision-making to optimize the object's future behavior (i.e., improve treatment outcomes), the inventors propose that the methodology is a practical manifestation of digital twins in tumor treatment. In this approach, population-based training of models is not required. Instead, the approach integrates individual patient multi-parametric MRI data (Figure 17A) obtained at multiple time points during treatment with biology-based mathematical modeling. Specifically, two frameworks were constructed to determine the predictive utility of each patient's digital twin (Figure 17B). In framework 1, the inventors seek to use digital twins to predict the outcome of a single NAST regimen (i.e., A / C). Specifically, we evaluate the accuracy of the digital twin to predict global indices related to changes in tumor burden, and spatiotemporally resolved tumor dynamics at the end of A / C (Figure 17C). In Framework 2, we seek to use the digital twin to predict the outcome of the entire NAST (i.e., both A / C and paclitaxel). Specifically, we evaluate the accuracy of the digital twin to predict individual pCR or non-pCR status at the completion of NAST (Figure 17D).

[0189] Patients and MRI Data: Treatment-naive patients with biopsy-confirmed TNBC were enrolled in the Institutional Review Board-approved prospective clinical trial "A Robust TNBC Evaluation FraMework to Improve Survival" (ARTEMIS, NCT02276442). Patients who signed informed consent between June 2018 and January 2020, had clinical stage I-III disease, had completed serial MRI scans, completed NAST, and had known postoperative pathology were included in the study (n=56).

[0190] Each patient underwent multiparametric MRI scans before (baseline / V1), after 2 cycles (V2), and after 4 cycles (V3) of standard-of-care (neoadjuvant) doxorubicin and cyclophosphamide (A / C). (Each "cycle" of A / C is 2 weeks; see Figure 17A.) Patients with progressive disease or <70% reduction in tumor volume at the end of A / C were offered the opportunity to enroll in a biomarker-guided clinical trial using targeted biotherapy / chemotherapy to complete therapy (n=9). Patients not meeting criteria for suboptimal response to A / C were recommended to continue standard-of-care paclitaxel weekly for 12 cycles (n=43) or twice every 3 weeks for 4 cycles (n=4). (The exact paclitaxel regimen for each patient was determined by the physician.) All patients underwent surgery after NAST. Postoperative pathology was used to classify patients as pCR or non-pCR. pCR was defined as the absence of residual invasive carcinoma and carcinoma in situ on hematoxylin and eosin evaluation of the completely resected breast specimen and all sampled regional lymph nodes after completion of NAST.

[0191] MRI was performed on a GE Discovery MR750 or MR750w whole-body scanner (GE Healthcare) equipped with an 8-channel bilateral breast coil. In particular, DCE-MRI data were acquired using a 3D DISCO sequence with a bipolar readout (FIG. 18A) and the following scan parameters: field of view = 30 × 30 cm 2, matrix size = 320 × 320, slice thickness / spacing = 3.2 / -1.6 mm, number of slices = 140, flip angle = 12°, repetition / echo time = 8 / 2 ms. After one pre-contrast phase was obtained, a single bolus injection (0.1 mL / kg at approximately 2 s, followed by a saline flush) of contrast agent (Gadovist, Bayer Healthcare) was performed at the beginning of the post-contrast acquisition. The time resolution of the DISCO series ranged from 8 to 15.5 s (median = 11 s) depending on the slice coverage, resulting in a number of post-contrast phases that varied from 32 to 64. DW-MRI was acquired using a 2D spin-echo sequence of the affected breast with the following scan parameters: field of view = 16 × 16 cm 2 , matrix size = 80 × 80, slice thickness / spacing = 4 / 0 mm, number of slices = 16, flip angle = 90°, repetition / echo time = 4000 / 70 ms. The b values ​​used were 100 and 800 s / mm 2 Apparent diffusion coefficient (ADC) maps were calculated using a GE AW server (v3.2, GE Healthcare, Milwaukee, WI). Tumors were manually segmented by two board-certified breast radiologists with 4–12 years of experience (authors MB, RMM). Tumor segmentations were reviewed by two breast specialist-trained radiologists with 19 and 20 years of experience (authors GMR, BEA). Whole tumor volumes were segmented at two initial stages of DCE-MRI using a home-based software package developed in MATLAB (R2021b, Mathworks, Natick, MA). All segmentations were further refined using the package's thresholding tools to exclude non-tumor voxels as determined by the radiologists. Necrotic regions and artifacts from biopsy clips were manually segmented and excluded by two radiologists (authors MB, RMM).

[0192] Image Processing: According to various potential embodiments, all MRI data from each patient was processed through a pipeline consisting of three components: 1) pre-processing, 2) inter-visit registration, and 3) post-processing (Figure 18B-D). This highly automated pipeline allows for efficient processing of multi-visit, multi-parametric MRIs with minimal user input.

[0193] First, the multiparametric images were colocalized to the same imaging grid, and slice and voxel locations were aligned. Specifically, the DCE-MRI acquired from both sides was cropped to the DW-MRI field of view covering the diseased breast, and the slice and ADC maps of the DW-MRI were linearly interpolated to match the slice locations of the DCE-MRI. Rigid registration was applied to the DCE-MRI to align all phases within one scan image to the pre-contrast phase (MATLAB function, imregtform). Rigid registration was applied between the colocalized DCE-MRI and DW-MRI data to remove small discrepancies between the image volumes.

[0194] Next, image registration between visits was performed to account for changes in breast tissue shape and patient position between MRI visits. Specifically, registration was performed to align the V1 and V3 images to the V2 image. The algorithm consisted of a rigid registration of the tumor ROI for initial alignment, followed by a non-rigid registration of a deformable b-spline of the entire breast with a stiffness penalty on the tumor region. This stiffness penalty was imposed to maintain the tumor volume and shape across all visits. The registration was developed based on MATLAB and the open-source command line software Elastix.

[0195] Third, post-processing was performed in preparation for subsequent predictive modeling. Specifically, semi-automated segmentation of the breast contour was performed in pre-contrast frames of DCE-MRI based on a manually chosen intensity threshold, followed by smoothing of the segmented mask edges (MATLAB function, imgaussfilt). Two-class k-means clustering (MATLAB function, kmeans) was used to segment fibroglandular and adipose tissues in each pre-contrast DCE-MRI. Enhancement for each DCE-MRI was calculated by subtracting the pre-contrast phase from the average of the post-contrast phase. Tumor cellularity maps N(x,t) were estimated based on the measured ADC maps of each MRI visit.

number

[0196] In the formula, ADC W is the ADC of free water (3 × 10 -3 mm 2 / s), ADC(x,t) is the ADC value of a voxel at position x and time t, and ADC min is the minimum (positive) ADC value in a patient's tumor across all visits. The capacity, θ, describes the maximum number of tumor cells that can physically fit within a voxel, which was determined by assuming a spherical packing density of 0.7405 and a nominal tumor cell radius of 10 μm.

[0197] Image-guided biology-based modeling: Applicants developed a biology-based mathematical model that represents the spatiotemporal resolved dynamics of tumor growth and response to NAST. Specifically, a reaction-diffusion type partial differential equation is used to describe the evolution of tumor cells N(x,t) in response to therapy, as shown in Equation (T2).

number

[0198] where a detailed list of variables, parameters, their definitions and assignments are given in Table T3 below. The first term on the right hand side of equation (T2) describes the migration of tumor cells and the compression of the surrounding tissue by a diffusion process coupled with the mechanical properties of the tissue via equation (T3). D(x,t)=D0e -rσ(x,t) (T3)

[0199] where Do is the diffusion coefficient of tumor cells in the absence of external forces. The exponential term reduces the mobility of tumor cells via the von Mises stress σ(x,t) and the stiffness of the surrounding tissues due to the empirical coupling constant γ. The von Mises stress is calculated for fibroglandular and adipose tissues in the breast, with fibroglandular tissue being assigned a greater stiffness than adipose tissue. The technical details regarding mechanically coupled diffusion are found in “Tumor Cell Growth” (second term on the right hand side of equation (T2)), which is described by logistic growth with a spatially varying growth rate k(x) and an overall carrying capacity θ. The effect of the administered therapy (third term on the right hand side of equation (T3)) is modeled as the treatment-induced mortality rate λi(x,t) of tumor cells. Mortality is determined by the exponential decay of the concentration of the administered drug and its effectiveness.

number

[0200] In the formula, a i i 番目 is the efficacy of the drug administered at time t, where i=1, 2, and 3 refer to doxorubicin, cyclophosphamide, and paclitaxel, respectively. Since each drug is administered multiple times during a therapeutic regimen, the total effect of the drug at time t is the cumulative effect of all administration cycles and their decays. That is, i 番目 Drugs j 番目 The administration of i,j In total, J i The decay in efficacy of the drug from each dose was β i The decay rate of this therapy is expressed as 番目 caused by administration of i 番目 The spatial distribution of the drug C(x,τ i,j) is determined by the enhancement of DCE-MRI. Specifically, τ i,j From the last DCE-MRI data collected prior to the injection at , the area under the DCE-MRI time course is calculated at each voxel and normalized by the maximum value within the tumor. The voxel-wise normalized area under the curve represents a map of the drug concentration induced by this injection (see Supplementary Section 1.1 for details).

[0201] Equations (T2-T4) are personalized by the imaging and clinical data of each patient. Specifically, the outline of the computational domain is determined by the segmented breast contour, tumor, and fibroglandular / adipose tissue. The tumor cellularity map at each imaging time point is determined from the ADC map obtained at the corresponding time point via equation (T1). Successive cellularity maps are calculated by the model parameters (i.e., k(x), D0, a i , and β i ) for patient-specific calibration. Additionally, DCE-MRI acquired during NAST was used to update the spatial distribution of the administered drug and the mechanical properties of the breast tissue. The model constrained by the patient-specific data was implemented in MATLAB and solved with finite differences. Details of the specific numerical implementation can be found in Jarrett AM, Kazerouni AS, Wu C, Virostko J, Sorace AG, DiCarlo JC, et al. Quantitative magnetic resonance imaging and tumor prognosis in breast cancer patients in the community setting. Nat Protoc 2021;16(11):5309-5338 and Jarrett AM, Hormuth II DA, Wu C, Kazerouni AS, Ekrut DA, Virostko J, et al. Evaluation of patient-specific neoadjuvant regimens for breast cancer via mathematical models constrained by quantitative magnetic resonance imaging data. Neoplasia 2020;22(12):820-830.

[0202] Digital Twin Framework: As shown in Figure 17C, framework 1 focuses on using digital twins to predict outcomes of a single NAST regimen (i.e., A / C). Specifically, for each patient, processed images from V1 and V2 are imported along with the treatment regimen and calibrated to the formula (T2-T4). The calibrated model is a digital twin that represents patient-specific pathophysiological attributes of tumor growth and response, including pre-treatment tumor shape and cellularity, tumor cell proliferation rate and mobility, and efficacy and decay of administered drugs (A / C). Because efficacy and decay rates are strongly coupled and difficult to calibrate simultaneously at only two time points, the decay rates of A / C were randomly sampled five times from the range described in the published literature.

number

[0203] The prediction accuracy of Framework 1 is evaluated temporally and spatially. Temporal accuracy is evaluated by the agreement between predicted and measured global metrics (i.e., TTC and TTV). In particular, the concordance correlation coefficient (CCC, see Supplementary Section 1.2) is calculated between the predicted and measured changes in TTC at the end of A / C. Similarly, the CCC is calculated between the predicted and measured changes in TTV. Spatial accuracy is evaluated by the difference between the predicted and measured spatially resolved tumor cell distributions. In particular, for each patient, the percent change in tumor cell count from baseline (V1) to the end of A / C (V3) can be calculated at each location x. We use the predicted and measured tumor cell distributions (i.e., ΔTCD p (x) and ΔTCD m The change from (x) was calculated, and the spatially resolved difference was ΔTCD p (x)-ΔTCD m The mean and 95% confidence interval (Cl) of the median difference from each patient is given by (x). For each patient, this difference is reported as the median and interquartile range within the original tumor area. The difference across the cohort is assessed by the mean and 95% confidence interval (Cl) of the median difference from each patient. (Considering that five predicted values ​​are given for each patient, these estimates are based on the median of the five predicted values.)

[0204] As shown in FIG. 17D, Framework 2 focuses on using the digital twin to predict the outcome of the entire NAST (i.e., both A / C and paclitaxel). Specifically, for each patient, processed images from V1, V2, and V3 are imported into the mechanism-based model for initialization and calibration. The calibrated model is the digital twin. (Note that the sampling scheme of Framework 1 is not necessary for Framework 2, since the three time points of data allow the efficacy and decay rate to be simultaneously calibrated within the same range as assumed in Framework 1.) Because imaging data were not available during the paclitaxel regimen, the efficacy of paclitaxel was calculated using the literature value, α3=0.3 days. -1and assumed the decay rate of paclitaxel to be the average of the calibrated A / C decay rates.) β3 = (β1 + β2) / 2. (For details, see Supplementary Section 2.1.) By applying paclitaxel to the digital twin, we predict patient-specific spatiotemporal resolved tumor cell distribution, TTC, and TTV at the end of NAST.

[0205] The output of Framework 2 is evaluated by the digital twin's ability to distinguish between pCR and non-pCR. Specifically, a receiver operating characteristic (ROC) analysis is performed on the predicted TTC and TTV. We report the area under the ROC curve (AUC), sensitivity (i.e., ability to correctly identify residual tumor in the final surgical pathology), and specificity (i.e., ability to correctly identify pCR in the final surgical pathology) based on optimal cutoffs. Additionally, the AUC / sensitivity / specificity from the predicted TTC and TTV was compared to that obtained by the measured TTC and TTV. The 95% CI of the AUC was calculated and compared via the DeLong method, and P<0.05 was considered statistically significant.

[0206] Results,Framework 1: Patient-specific prediction of spatiotemporal response to A / C:A cohort of 50 patients was used in framework 1. Six patients were excluded from the overall patient cohort (n=56) due to image acquisition errors or artifacts (n=1), inter-visit registration failure due to large changes in breast shape between visits (n=1), and complete tumor response in V2 that did not provide data for model calibration (n=4).

[0207] Each patient's digital twin provided a range of estimated treatment efficacy of A / C (e.g., Figures 19A-B) and simulated a range of treatment outcomes (e.g., Figures 19C-D). Specifically, Figure 19C shows a patient who had a suboptimal response to A / C (i.e., V3 imaging showed less than 70% reduction in tumor volume). The digital twin showed that in V3, the mean (range) TTC and TTV were 3.17 x 10, respectively. 8 (3.05×10 8 ~3.29×108 ) cells and 3.19 × 10 3 (3.15×10 3 ~3.23×10 3 )mm 3 This corresponds to an expected percent reduction in TTC and TTV in V3 of 49.55% (47.60%-51.50%) and 36.14% (35.44%-36.84%), respectively. In comparison, the measured percent reduction in TTC and TTV in V3 was 38.99% and 32.58%, respectively. In contrast, Figure 19D shows a patient who responded well to A / C (i.e., V3 imaging showed less than 70% reduction in tumor volume). The digital twin predicted that TTC and TTV in V3 were 3.82 × 10 6 (1.69×10 6 ~5.94×10 6 ) cells and 12.63 (0.00-25.27) mm 3 19E). The CCC between the predicted and measured changes in TTC in V3 was 0.95 (Figure 19E). The CCC between the predicted and measured changes in TTV was 0.94 (Figure 19F). These results indicate high prediction accuracy and precision (i.e., uncertainty in the model predictions; see Supplementary Section 2.2 for interpretation) of the time dynamics of TNBC response to neoadjuvant A / C (Table T1). [Table 11]

[0208] Importantly, the personalized digital twin provides not only the overall metrics summarized in the previous paragraph, but also the spatiotemporal evolution of each patient's tumor. Figures 20A and 20B show the measured and predicted tumor cell distributions, respectively, from the central slices of the same two exemplary patients. Figures 20C and 20D show 3D renderings of the measured and predicted tumor volumes, respectively. In both cases, the digital twin successfully captures the lack of response (Figure 20A) or the presence of response (Figure 20B). Quantitatively, for the first patient, the median (interquartile range) of the difference between the predicted and measured changes in tumor cell distribution in V3 is -3.30% (-22.07% to 0.00%). For the second patient, the median (interquartile range) of the difference is 0.00% (0.00% to 0.00%). The median difference across all patients had a mean (95% CI) of 0.20% (-20.35% to 20.75%) in V3 (Figure 20E). These results demonstrate the high predictive precision and accuracy of spatially resolved prediction of tumor cell distribution in TNBC patients who responded to neoadjuvant A / C (Table T1).

[0209] Results,Framework 2: Patient-specific prediction of final pathological response:,For framework 2, a cohort of 37 patients (18 pCR, 19 non-pCR) was used.,After excluding 6 patients from the entire cohort (n=56) as done for,framework 1, 13 more patients were excluded due to presumed,nonresponse to further standard of care chemotherapy and enrollment in,clinical trials (n=9, part of the ARTEMIS schema), missing,schedule for paclitaxel (n=2), and image acquisition errors (n=2,,ADC maps not covering the whole tumor in V3).

[0210] For each patient, the personalized digital twin estimated treatment efficacy based on V1-V3 images (Figure 21A) and depicted tumor dynamics in response to A / C. These results were then used to predict response to paclitaxel and thus final treatment outcomes after all NASTs (Figure 21C). For example, Figure 21C shows that the digital twin can be calibrated during the A / C regimen and used to predict the absence of regrowth during paclitaxel in patients who actually achieved pCR after NAST. In contrast, Figure 21D shows that the digital twin captured the initial response and subsequent regrowth during A / C and predicted robust regrowth pre- and during paclitaxel, resulting in predicted TTC and TTV values ​​after NAST of 9.03 × 10, respectively. 8 , TTV is 5.21 × 10 3 mm 3 This shows that:

[0211] As shown in Table T2, in the entire cohort, both predicted TTC and TTV from the digital twin (based on models calibrated with V1-V3 data) achieved an AUC of 0.89 (0.78-0.99) for discriminating between pCR and non-pCR (based on postoperative pathology). In contrast, measured TTV or TTC (based on V3 data) achieved an AUC of 0.78 (0.62-0.94) for discriminating between pCR and non-pCR. Furthermore, using predicted TTC and TTV, the specificity was 0.95 and 0.89, respectively. In contrast, using measured TTV or TTC, the specificity was only 0.79. These results indicate that the digital twin improved the prediction of final response compared to directly measured data. Specifically, the AUC improved by 14.28% (P = 0.04, significant) for TTC and 13.83% (P = 0.07) for TTV. Specificity improved by 20.25% and 12.66% for TTC and TTV, respectively. Sensitivity was unchanged. [Table 12]

[0212] Discussion of Various Embodiments: Potential embodiments of the present disclosure provide a digital twin approach to achieve early, patient-specific, spatio-temporal resolved prediction of TNBC patient response to neoadjuvant doxorubicin, cyclophosphamide, and paclitaxel. The approach was based on a biology-based mathematical model calibrated with multi-visit, multi-parametric MRI acquired for individual patients. Thus, this methodology represents a major step away from population-based predictions and towards individual-based predictions.

[0213] Framework 1 shows how the digital twin utilizes imaging data from individual patients acquired early in neoadjuvant A / C to predict tumor status at the end of A / C with great accuracy. The CCC between predicted and measured values ​​for total tumor burden and total tumor volume were 0.95 and 0.94, respectively. This strongly indicates that early changes during the A / C regimen contain enough information to calibrate the digital twin and confidently predict tumor response at the end of A / C. This observation is consistent with previous reports that demonstrated that indices from early treatment MRI are strong predictors of NAST response in breast cancer. This provides strong support that our approach can be used to adjust treatment regimens on a patient-specific basis. For example, our approach can be applied after the first two cycles of A / C to predict whether further A / C administration should be continued or alternative interventions should be considered.

[0214] Framework 2 shows how the digital twin utilizes imaging data from individual patients acquired during the A / C portion of NAST to accurately predict the patient's final pathological status (i.e., pCR or non-pCR) upon completion of all NAST, with an AUC of 0.89. Importantly, using tumor volume measured at the end of A / C only yielded an AUC of 0.78. Thus, the digital twin provided a significant improvement in AUC (14%). The accuracy of the digital twin is also improved compared to previous MRI-based predictions of breast cancer response to NAST. Both our previous study and the l-SPY study reported the best discrimination between pCR / non-pCR in TNBC using optimized in-treatment FTV, with an AUC of 0.85. Furthermore, adding post-NAST ADC to the FTV showed improved response prediction in TNBC, increasing the AUC from 0.71 to 0.81. Pharmacokinetic parameters (i.e., k ep Combining the ADC measured after one cycle of NAST with the pretreatment DCE-MRI yielded an AUC of 0.88 in breast cancer. Moreover, the predictive accuracy of digital twins is comparable to state-of-the-art machine learning (ML)-based predictions. For example, Ravichandran et al. applied a convolutional neural network (CNN) to predict pCR from pretreatment DCE-MRI, achieving an AUC of 0.77 in a total of 166 breast cancer patients. A more recent CNN-based study using both pretreatment and posttreatment DCE-MRI achieved an AUC of 0.91 in a cohort of 42 breast cancer patients.

[0215] Importantly, embodiments of the disclosed digital twin approach have several inherent advantages compared to ML algorithms. First, ML methods rely on access to large patient populations to train the algorithms, and this training dataset must contain all pathophysiological features relevant to the disease under investigation, be annotated with high quality, and be generalizable from one population to the next. In contrast, our approach does not require population-based training or annotation labels, as patient-specific data is used to calibrate a biology-based model to make patient-specific predictions. Second, ML can be difficult to interpret biologically due to the complexity of the modeling functions. In contrast, digital twins provide accurate predictions for mechanistic interpretation of tumor development during NAST, as well as pCR status at the end of NAST. For example, our modeling framework can capture early and subsequent responses to A / C for patients with very different response kinetics, as depicted in Figures 19A-D. Third, digital twin parameters quantify and elucidate observed tumor response kinetics, thus offering another potential application of predicting response to multiple candidate therapy regimens. Therefore, a digital twin built on quantitative imaging data could provide a practical way to optimize individual treatment plans and facilitate truly personalized cancer care at the early stages of NAST.

[0216] Potential Variations: First, in various embodiments, instead of utilizing DCE-MRI to inform the spatial distribution of the delivered drug, a tumor response model (i.e., Eq. (T2)) can be coupled with a drug delivery and / or tumor angiogenesis model to more accurately estimate drug distribution. Second, in various embodiments, instead of assuming that cell death is simply proportional to drug concentration (i.e., the last term on the right-hand side of Eq. (T2)), a model that takes into account detailed therapy mechanisms and pharmacokinetics can improve prediction accuracy. Third, the methodology in this "Additional Examples" section currently cannot predict lymph node or axillary invasion, which limits the accuracy of predicting pathological response or residual cancer burden (see Supplementary Section 2.3 for additional analysis). Of course, building a more comprehensive model would require more assumptions or more measurements, so a careful balance between model complexity and prediction accuracy must be sought.

[0217] Although the dataset used in this study in the "Additional Examples" section is much larger and more homogeneous than previous pilot studies, there are still some points about the cohort that may affect the analysis. Framework 1 excluded 4 patients with no visible tumor in V2. This decision was made to avoid overestimating the accuracy of the prediction, as such patients also had no visible tumor in V3, resulting in a 100% accurate prediction without actually testing the model. Additionally, there was no pathological evaluation at intermediate time points, which led to a lack of "ground truth" for Framework 1. Framework 2 excluded 9 patients due to enrollment in other trials, which enriched the cohort for pCR patients. This may lead to an overestimation of the accuracy of distinguishing pCR / non-pCR (see Supplementary Section 2.4 for additional analysis).

[0218] Processing (e.g., segmentation, registration) steps are sources of potential errors that can propagate through the modeling pipeline and lead to bias in predictions. Detailed investigations suggested that unexpected changes in ADC values ​​and distributions, as well as the appearance of necrotic regions, are potential sources of error (see Supplementary Section 2.5). Additionally, Framework 1 includes sampling of drug decay rates, which introduces uncertainty in the predicted tumor response and limits the accuracy of the final pathology prediction (see Supplementary Section 2.6). However, compared to previous attempts to simultaneously calibrate drug efficacy and decay, this procedure not only ensures a more robust model calibration, but also allows for the uncertainty in tumor dynamics to be quantified and interpreted. Another source of uncertainty is the setting of paclitaxel drug efficacy and decay rates in Framework 2, due to the lack of imaging in the paclitaxel portion of NAST. One solution would be to incorporate more measurements during NAST, especially after alternating therapies, so that the digital twin can be updated to maintain accurate predictions. Of course, it is important to note that the goal of the digital twin is not to provide a perfect reproduction of the patient's situation. Rather, a realistic goal is to provide a precise, practical formalism that provides clinically actionable insights. In the present contribution, the inventors achieve this goal in the context of predicting the response of early-stage triple-negative breast cancer to neoadjuvant systemic chemotherapy.

[0219] This study also supports the value of longitudinal MRI in cancer care. Currently, only pretreatment MRI is standard for the evaluation of breast cancer patients. Although increasing to multiple MRIs would increase the cost of imaging, the benefit of early detection of chemotherapy resistance (for example) would help avoid unnecessary toxicity and cost. Both this study and previous studies have shown that follow-up imaging after the first 1–3 cycles of NAST can help with early prediction of response. The inclusion of even one follow-up MRI (enough to allow calibration of the model) would be of great advantage.

[0220] These embodiments provide a patient-specific digital twin through the integration of longitudinal multi-parametric MRI data with biology-based mathematical modeling. This technique accurately captures the spatiotemporal response of TNBC to NAST and achieves high accuracy and specificity for predicting the final pathological status of each individual patient. The success of this approach demonstrates the potential of digital twins to shift the paradigm from response assessment to prediction and ultimately response optimization. [Table 13]

[0221] Supplementary Section 1: Supplementary Methods Supplementary Section 1.1 Determining the spatial distribution of drugs As shown in equation (T4), the j 番目 i caused by administration of 番目 The spatial distribution of the drug is C(x,τ ij ) as determined by the enhancement observed via DCE-MRI measurements (3). Specifically, for each voxel in the breast, we calculate the dynamic area under the curve (AUC) of the baseline-subtracted DCE time course (i.e., the DCE-MRI intensity time course minus the intensity of the pre-contrast frames), and then normalize this voxel-wise AUC by dividing it by the maximum AUC value within the tumor. Thus, the voxel-wise normalized AUC is expressed as a function of time τ ij Figure 1 shows maps of infusion-induced peak drug concentrations in the 100-well plateau at 100 nm. This method allows us to capture the spatially heterogeneous distribution of each drug, which is determined by the perfusion characteristics of individual patients.

[0222] For each injection, the calculations described in the last paragraph are performed on the last DCE-MRI scan image collected before the injection of the drugs. Specifically, in this study, the spatial distribution of doxorubicin and cyclophosphamide caused by the first and second cycles is calculated from the DCE-MRI data of the V1 scan. The spatial distribution of doxorubicin and cyclophosphamide caused by the third and fourth cycles is calculated from the DCE-MRI data of the V2 scan. The spatial distribution of paclitaxel caused by each cycle is calculated from the DCE-MRI data of the V3 scan. This method allows for the consideration of changes in vascular structure throughout the therapy regimen.

[0223] The above methodology is based on two assumptions: 1) the drug concentration at each location in the domain is approximately proportional to the auc of the time course of the local DCE-MRI intensity, and 2) the washout and efficacy decay of the delivered drug can be approximately represented by a voxel-wise exponential decay. We recognize that these assumptions, while realistic, may not be the most accurate. In future work, we hope to combine the tumor response model (i.e., Eq. 2 from the main text) with a rigorous description of drug delivery and tumor angiogenesis (4, 5) to achieve a more accurate spatial and temporal representation of drug distribution.

[0224] Supplementary Section 1.2: Concordance correlation coefficient The concordance correlation coefficient (CCC) (6) measures the degree of concordance between two samples (i.e., x and y) via:

number

[0225] During the ceremony,

number

number

[0226] Supplementary Section 2: Supplementary Results Supplementary Section 2.1: Uncertainty in Drug Efficacy As shown in the methodological discussion in the Additional Examples section, there are uncertainties in both Framework 1 and Framework 2 in explaining drug efficacy and attenuation.

[0227] For Framework 1, the mathematical model is calibrated using images of V1 and V2. 番目 The efficacy and decay of a drug are expressed as efficacy (α i ) and decay (β i ) rate is strongly coupled to equation (T2) (i.e.,

number

number

[0228] Fortunately, based on the results of this study, this approach introduces only minor uncertainty in the model output of predicted tumor dynamics in Framework 1 (see Figure 19). Additionally, even with uncertainty, the model does a good job of capturing the substantial differences in tumor dynamics between patients with different responses. As shown in Figures 19A and 19B of the main text, poor responders (i.e., tumor volume shrinks by 33% at the completion of A / C) have a much higher tumor cell proliferation rate than good responders (i.e., tumor volume shrinks by 100% at the completion of A / C): median k(x) = 0.05 days. -1 vs. 0.004 days -1 Additionally, poor responders showed high sensitivity to cyclophosphamide but low sensitivity to doxorubicin (median α1 = 0.24 days). -1 , α2=0.60 days -1 In contrast, good responders showed moderate sensitivity to both doxorubicin and cyclophosphamide (median α1 = 0.47 days). -1 , α2=0.34 days -1Future studies may be able to more fully address the question of how much uncertainty in these parameters can be tolerated before it begins to affect the final model predictions by combining a large sample set of drug decay rates with a Bayesian approach to quantifying error propagation from model parameters to outcomes.

[0229] For Framework 2, the mathematical model is calibrated using images from V1, V2, and V3. The data from the three time points provide a validity rate (α i ) and decay rate (β i ) can be independently calibrated. Therefore, we chose to simultaneously calibrate {α1, α2, β1, β2} without using a sampling scheme (described in the last paragraph) to avoid the uncertainty of the A / C regimen. However, another uncertainty exists for the second regimen, since no data are available to calibrate the decay and efficacy rates of paclitaxel (i.e., α3 and β3). To address this situation, the efficacy of paclitaxel is set to the literature value. α3=0.3 days -1 (6). The decay rate of paclitaxel was obtained as the average of the calibrated A / C decay rates: β3=(β1+β2) / 2.

[0230] This uncertainty may undermine potential heterogeneity in patient sensitivity to paclitaxel therapy in Framework 2. This limitation arose from a lack of monitoring data that would not allow updating patient-specific digital twins when therapy regimens were changed from A / C to paclitaxel. In general, we lack data types that would help inform each patient's sensitivity to paclitaxel. Future efforts to address this point are therefore warranted. One option is to monitor patient response longitudinally (multiple times), especially when therapy is changed, allowing the digital twin to be continually calibrated and refined to maintain accurate expectations. Another is to integrate multi-dimensional / multi-type data (e.g., genetic profiles) to inform model parameter assignment.

[0231] Supplementary Section 2.2: Interpreting uncertainties in model predictions In Framework 1, the range of model predictions represents the uncertainty in individual patient tumor response kinetics under alternative pharmacokinetic (i.e., drug decay rate) scenarios. As evident in Figures 19E and 19F, most patients in the cohort have small uncertainty in their predicted tumor response kinetics, while a few patients have relatively large uncertainty.

[0232] For example, FIG. 22 shows the predicted TTC time course (median and range) for the patient showing the largest range in FIG. 19E.

[0233] Two TTC time courses within the expected range are highlighted: one assuming fast drug decay (dashed magenta curve; sampled β = 0.7 days); -1 and β2 = 5.2 days -1 ) and the other assumes slow drug decay (dashed black curve. Sampled β1 = 0.1 days -1 , and β2 = 1.7 days -1 ) within the first A / C cycle (indicated by the green arrow), trajectories with rapid drug decay also identify a strong and rapid response and subsequent strong regrowth of the tumor. In contrast, trajectories with low drug decay identify a relatively slow response with little apparent regrowth. Such dramatically different dynamics in the sampled drug pharmacokinetics in this patient result in a "scissor-shaped" range within the first AC cycle, and also in a large uncertainty in the final prediction of TTC (similar to the case of TTV). This is very different from the patient seen in Figures 19C and 19D, which show a "spindle-shaped" range within the first AC cycle with a small uncertainty in the final prediction. Specifically, for the patient in Figure 19C, all TTC trajectories predict a rapid response and strong regrowth, and for the patient in Figure 19D, all TTC trajectories predict weak regrowth.

[0234] This observation means that by "seeing" the V1 and V2 data, our model "thinks" that there is a wide range of possibilities underlying the tumor dynamics that could lead to imaging-measured tumor changes in this particular patient. This is indeed an interesting insight, since one promising application of digital twin techniques is to guide personalized monitoring tests. Patients who show a large uncertainty in their early predictions may be recommended to undergo more frequent follow-up imaging (or other) tests, whereas patients with very small uncertainties may be recommended less frequent follow-ups. In future studies, the inventors would like to further explore the imaging characteristics and molecular profiles from intermediate biopsies to better understand which features actually drive the large uncertainty in the tumor response dynamics in these patients.

[0235] Supplementary Section 2.3: Comparison of TTV / TTC and pathological residual cancer burden In addition to predicting pCR status, we further compared the model predictions and imaging measurements with pathological residual cancer burden (RCB). Specifically, RCB was estimated from routine pathological sections of the primary breast tumor site and regional lymph nodes after completion of neoadjuvant therapy. Of the 37 patients included in Framework 2, 19 patients were in pCR, 5 patients were in RCB-I, 11 patients were in RCB-II, and 2 patients were in RCB-III. RCB-I, -II, -III classes are considered as non-pCR groups.

[0236] Figure 23 shows the predicted TTV at end of treatment (EoT), the predicted TTC at EoT, the measured TTV at V3, and the measured TTC at V3 versus RCB class. Wilcoxon tests were also performed between each pair of RCB classes for each index investigated (in Figure 23, * indicates P<0.05 and ** indicates P<0.01).

[0237] For the predicted TTV at EoT, significant differences were found between the pCR group and each non-pCR group; P = 0.003 for pCR vs RCB-I, P = 0.001 for pCR vs RCB-II, and P = 0.046 for pCR vs RCB-III, respectively. Similarly, for the predicted TTC at EoT, significant differences were found between the pCR group and each non-pCR group; P = 0.003 for pCR vs RCB-I, P = 0.001 for pCR vs RCB-II, and P = 0.034 for pCR vs RCB-III, respectively. However, no significant differences were found in predicted TTV or TTC between either pair of non-pCR groups. In contrast, for TTV measured at V3, significant differences were found only between the pCR group and the RCB-II group, P = 0.006. Similarly, in V3, the measured TTC was significantly different only between the pCR and RCB-II groups, P = 0.006. No significant differences in measured TTV or TTC were found between other pairs of RCB classes.

[0238] These results support the position that digital twin-based prediction improves the ability to distinguish pCR from non-pCR. However, neither prediction to the end of NAST nor imaging measurements at the end of A / C can provide adequate discrimination between the different non-pCR RCB classes. This is probably due to the fact that neither the imaging measurements nor the prediction scheme in the Additional Examples section take into account lymph node involvement.

[0239] Supplementary Section 2.4: Stratified analysis of the nine patients excluded from Framework 2 A stratified analysis was performed taking into account the nine patients who were excluded from Framework 2 due to enrollment in other clinical trials, and the inventors performed a stratified analysis on a subset of nine patients to assess the potential impact of this exclusion on the evaluation.

[0240] We repeated the evaluation of Framework 1 on this subset of nine patients. The CCC between the measured and expected changes in TTC was 0.83, the CCC between the measured and expected changes in TTV was 0.79, and the mean (95% Cl) difference between the measured and expected changes in tumor cell distribution was -1.33% (-41.83% to 39.18%). Thus, the predictive accuracy of Framework 1 in this subset of nine patients is lower than that of the entire cohort (shown in Table T1). This observation suggests that the exclusion of these nine patients from Framework 2 may lead to an overestimation of the accuracy of distinguishing pCR from non-pCR.

[0241] Supplementary Section 2.5: Investigation of patients with suboptimal predictive accuracy As shown in FIG. 20E, four patients (i.e., patients 5, 11, 26, and 49) were identified as having suboptimal prediction accuracy (i.e., predictions were outside the mean 95% Cl range for the cohort). Careful investigation showed that the suboptimal predictions in these patients were associated with unexpected changes in the spatial pattern of measured tumor response across treatment. For example, the panels in FIG. 24 show the measured and predicted tumor cell distributions on the central tumor slice of patient 26 in FIG. 20E.

[0242] Considering the top row, we can see that the tumor cells are distributed almost uniformly in V1 measurements, tend to accumulate in the inner posterior and outer anterior of the tumor in V2, and change to accumulate along the tumor periphery in V3. The change in cell distribution from V2 to V3 is completely unexpected based on the imaging measurements in V1 and V2. Thus, after calibration with V1 and V2, the model predicted tumor response in a spatial pattern of high proliferation around the inner posterior and outer anterior regions and high death in other regions, which led to a large difference compared to the measurements in V3. This change in the spatial pattern of tumor response across longitudinal imaging may be due to the evolution of the tumor microenvironment and clonal populations, which leads to a large difference in tumor response between early and late stages of AC therapy. Alternatively, it may be due to errors in the measured ADC. In future studies, the inventors would like to further investigate the image quality, imaging characteristics, and molecular profiles from intermediate biopsies of these patients to better understand tumor response to NAST.

[0243] Another factor that can cause unexpected changes in spatial patterns is the presence of large necrotic regions. For example, the panels in Figure 25 show the measured and predicted tumor cell distributions on the central tumor slice of patient 49 in Figure 20E. The appearance of large necrotic regions in V3 makes it difficult to predict the final tumor shape.

[0244] Supplementary Sections 2.6: 2-MRI Calibration and 3-MRI Calibration The study showed that using additional MR scans for calibration improved the accuracy of prediction and reduced the uncertainty. Consider the patient in Figure 21C as an example. Figure 26 plots the overall NAST prediction by a two-scan calibrated model (2405 and 2410) and a three-scan calibrated model (2415). The results show that the two-scan calibrated model successfully predicts the patient's early and strong response, but is unable to predict the complete shrinkage of the tumor. In contrast, including a third MRI scan in the calibration allows the model to fully capture the shrinkage of the tumor, leading to a successful prediction of pCR for this patient.

[0245] Non-limiting exemplary embodiments are provided herein. Embodiment A: A method includes acquiring, by a computing system, magnetic resonance imaging (MRI) data corresponding to a plurality of MRI scan images of an anatomical region including a tumor of a patient, the plurality of MRI scan images including a first image set obtained through a first scan performed prior to administration of a therapy to the patient and a second image set obtained through a second scan performed after administration of the therapy to the patient, the therapy including administration of a plurality of drugs; determining, by the computing system, tissue characteristics of tissue surrounding the tumor from the MRI data; registering, by the computing system, image-related data generated from the second image set to the image-related data generated from the first image set; determining, by the computing system, diffusion characteristics of the tumor based on the tissue characteristics; determining, by the computing system, growth characteristics of the tumor based on the tissue characteristics; determining, by the computing system, for each drug of the plurality of drugs, an effect of the drug on tumor cells; and generating, by the computing system, a score indicative of a predicted response of the tumor to the therapy based on the diffusion characteristics of the tumor, the growth characteristics of the tumor, and the determined effect of each drug on cells included in each voxel. Embodiment B: The method of embodiment A, further comprising performing tumor segmentation to identify a tumor region of interest (ROI) based on the MRI data prior to determining tissue characteristics. Embodiment C: The method of embodiment A or B, wherein the tissue characteristic is related to vascular structure within the tissue. Embodiment D: The method of any one of embodiments A to C, wherein the MRI data includes dynamic contrast-enhanced MRI (DCE-MRI) data, and the tissue characteristics are quantified based on the DCE-MRI data. Embodiment E: The method of any one of embodiments A to D, wherein the tissue trait is quantified based on a pharmacokinetic model. Embodiment F: The method of any one of embodiments A to E, wherein the tissue characteristics are quantified based on a fluid dynamics model. Embodiment G: The method of any one of embodiments A to F, wherein the tissue characteristic is quantified based on the Kety-Tofts model and / or a modified version of the Kety-Tofts model. Embodiment H: The method of any one of embodiments A to G, wherein the MRI data includes diffusion-weighted MRI (DW-MRI) and the method further includes generating a map of the apparent diffusion coefficient (ADC) of water. Embodiment I: A method according to any one of embodiments A to H, wherein registering the image-related data from the second set of images to the image-related data generated from the first set of images aligns the images to the map within a common domain. Embodiment J: The method of any one of embodiments AI, further comprising determining drug distribution within each voxel of the tissue. Embodiment K: The method of any one of embodiments A to J, wherein determining the drug distribution within each voxel of the tissue includes generating a normalized map of blood volume to define the initial drug distribution across the domain at the time of each administration of therapy. Embodiment L: The method of any one of embodiments A-K, wherein the diffusion characteristics correspond to diffusion of tumor cells mechanically associated with the material properties of the tissue via a physical stressor, thereby representing changes in the tumor that may cause deformation within the tissue. Embodiment M: The method of any one of embodiments A-L, wherein the growth characteristics are based on carrying capacity, which is related to the maximum number of tumor cells that can physically fit within a voxel. Embodiment N: The method of any one of embodiments A-M, wherein the growth characteristics are based on proliferation rates per voxel, and the proliferation rates are calibrated for each voxel within the patient's tumor ROI. Embodiment O: The method of any one of embodiments A-N, wherein the effect of each drug on the tumor cells is based on the concentration of the drug in the tissue. Embodiment P: The method of any one of embodiments A to O, wherein the effect of each drug on the tumor cells corresponds to the spatiotemporal distribution of each drug within the tissue. Embodiment Q: The method of any one of embodiments A-P, wherein the effect of each drug on the tumor cells is based on at least one of the drug's efficacy parameter α, the drug's washout parameter β over time after each administration, or the drug's initial concentration. Embodiment R: The method of any one of embodiments A-Q, wherein the effect of each drug on the tumor cells is based on the drug's efficacy parameter α, the drug's washout parameter β over time after each administration, or the drug's initial concentration. Embodiment S: The method of any one of embodiments A-R, wherein at least one of the efficacy parameter and the washout parameter is calibrated for a patient or drug. Embodiment T: The method of any one of embodiments A to S, wherein the calibration of a patient's washout parameters is constrained using boundaries defined from the range of the terminal elimination half-life of the drug. Embodiment U: The method of any one of embodiments A-T, wherein the MRI data includes dynamic contrast-enhanced MRI (DCE-MRI) data and the initial concentration is approximated based on the DCE-MRI data. Embodiment V: The method of any one of embodiments A to U, further comprising performing rigid or non-rigid intra-scan image registration of the MRI data in each of the first and second image sets. Embodiment W: The method of any one of embodiments A-V, further comprising determining a modified therapy based on a score indicative of the tumor's predicted response to the therapy. Embodiment X: The method of any one of embodiments A-W, further comprising administering a modified therapy based on a score indicative of the tumor's predicted response to the therapy. Embodiment Y: The method of any one of embodiments A to X, wherein the therapy and / or modification therapy is neoadjuvant therapy (NAT). Embodiment AA: A computing system comprising one or more processors and a computer readable memory having instructions configured to cause the one or more processors to: acquire MRI data corresponding to a plurality of MRI scan images of an anatomical region including a tumor, the plurality of MRI scan images including a first image set obtained through a first scan performed prior to administration of a therapy to the patient and a second image set obtained through a second scan performed after administration of the therapy to the patient, the therapy including administration of a plurality of drugs; determine from the MRI data one or more tissue attributes of tissue surrounding the tumor; register image-related data generated from the second image set to the image-related data generated from the first image set; determine diffusion characteristics of the tumor based on the tissue attributes; determine growth characteristics of the tumor based on the tissue attributes; determine, for each drug of the plurality of drugs, an effect of the drug on tumor cells; and generate a score indicative of a predicted response of the tumor to the therapy based on the diffusion characteristics of the tumor, the growth characteristics of the tumor, and the determined effect of each drug on cells included in each voxel. Embodiment BB: The system of embodiment AA, wherein the instructions are further configured to cause one or more processors to perform tumor segmentation to identify a tumor ROI based on the MRI data prior to determining tissue characteristics. Embodiment CC: A system described in embodiment AA or BB, wherein the tissue characteristic is related to vascular structures within the tissue. Embodiment DD: A system described in any one of embodiments AA to CC, wherein the MRI data includes dynamic contrast-enhanced MRI (DCE-MRI) data and tissue characteristics are quantified based on the DCE-MRI data. Embodiment EE: A system described in any one of embodiments AA-DD, wherein tissue characteristics are quantified based on a pharmacokinetic model, a fluid dynamics model, a Kety-Tofts model, and / or a variant of the Kety-Tofts model. Embodiment FF: A system described in any one of embodiments AA to EE, wherein the MRI data includes diffusion-weighted MRI (DW-MRI) and the instructions are further configured to cause one or more processors to generate a map of the apparent diffusion coefficient (ADC) of water. Embodiment GG: A system described in any one of embodiments AA to FF, wherein aligning the image-related data from the second set of images to the image-related data generated from the first set of images aligns the images to the map within a common domain. Embodiment HH: A system described in any one of embodiments AA to GG, wherein the instructions are further configured to cause one or more processors to estimate drug distribution within each voxel of the tissue. Embodiment II: A system described in any one of embodiments AA to HH, wherein determining the drug distribution within each voxel of the tissue includes generating a normalized map of blood volume to define the initial drug distribution across the domain at the time of each administration of therapy. Embodiment JJ: A system described in any one of embodiments AA to II, wherein the diffusion characteristics correspond to diffusion of tumor cells mechanically associated with the material properties of the tissue via a physical stressor, thereby representing changes in the tumor that may cause deformation within the tissue. Embodiment KK: A system described in any one of embodiments AA-JJ, wherein the growth characteristics are based on carrying capacity, which relates to the maximum number of tumor cells that can physically fit within a voxel. Embodiment LL: A system described in any one of embodiments AA to KK, wherein the growth characteristics are based on proliferation rates per voxel, and the proliferation rates are calibrated for each voxel within the patient's tumor ROI. Embodiment MM: A system described in any one of embodiments AA-LL, wherein the effect of each drug on tumor cells is based on the concentration of the drug in the tissue. Embodiment NN: A system described in any one of embodiments AA-MM, wherein the effect of each drug on tumor cells corresponds to the spatiotemporal distribution of each drug within the tissue. Embodiment OO: A system described in any one of embodiments AA to NN, wherein the effect of each drug on tumor cells is based on at least one of the drug's efficacy parameter α, the drug's washout parameter β over time after each administration, or the drug's initial concentration. Embodiment PP: A system described in any one of embodiments AA-OO, wherein the effect of each drug on the tumor cells is based on the drug's efficacy parameter α, the drug's washout parameter β over time after each administration, or the drug's initial concentration. Embodiment QQ: A system described in any one of embodiments AA to PP, wherein at least one of the efficacy parameters and the washout parameters is calibrated for a patient or drug. Embodiment RR: A system according to any one of embodiments AA to QQ, wherein the calibration of the patient's washout parameters is constrained using boundaries defined from the range of the drug's terminal elimination half-life. Embodiment SS: A system described in any one of embodiments AA to RR, wherein the MRI data includes DCE-MRI data and the initial concentration is approximated based on the DCE-MRI data. Embodiment TT: A system described in any one of embodiments AA to SS, wherein the instructions are further configured to cause the one or more processors to perform rigid and / or non-rigid intra-scan image registration of the MRI data in each of the first image set and / or the second image set. Embodiment UU: A system described in any one of embodiments AA to TT, wherein the instructions are further configured to cause one or more processors to determine modified therapy based on a score indicating the tumor's predicted response to the therapy. Embodiment VV: A system described in any one of embodiments AA-UU, wherein the therapy and / or modification therapy is neoadjuvant therapy (NAT).

[0246] As used herein, the terms "approximately," "about," "substantially," and similar terms are intended to have a broad meaning consistent with common and widely accepted usage by those of ordinary skill in the art to which the subject matter of this disclosure pertains. It should be understood by those of ordinary skill in the art who review this disclosure that these terms are intended not to limit the features described and claimed to the precise numerical ranges provided. These terms should therefore be interpreted as indicating that insubstantial or insignificant modifications or variations of the subject matter of the claims are believed to be within the scope of the invention as set forth in the appended claims.

[0247] It should be noted that the use of words "exemplary," "example," "potential," and variations thereof to describe various embodiments are intended to indicate that such embodiments are possible examples, representatives, or illustrations of possible embodiments (and such terms are not intended to imply that such embodiments are necessarily particular or best examples).

[0248] The term "coupled" and its variations as used herein means that two members are directly or indirectly coupled to each other. Such coupling may be static (e.g., permanent or fixed) or movable (e.g., removable or releasable). Such coupling may be achieved when the two members are directly coupled to each other, when the two members are coupled to each other using a separate intervening member and any additional intermediate members coupled to each other, or when the two members are coupled to each other using an intervening member integrally formed with one of the two members. When "coupled" or its variations are modified by an additional term (e.g., directly coupled), the general definition of "coupled" provided above is modified by the plain meaning of the additional term (e.g., "directly coupled" means that two members are coupled without a separate intervening member), resulting in a definition narrower than the general definition of "coupled" provided above. Such coupling may be mechanical, electrical, or fluid.

[0249] As used herein, the term "or," when used to connect a list of elements, is used in its inclusive (and not exclusive) sense, such that the term "or" refers to one, some, or all of the elements in the list. Conjunctions such as the phrase "at least one of X, Y, and Z" are understood to convey that an element may be either X, Y, Z, X and Y, X and Z, Y and Z, or X, Y and Z (i.e., any combination of X, Y, and Z), unless otherwise indicated. Thus, such conjunctive terms generally do not imply that a particular embodiment requires that at least one of X, at least one of Y, and at least one of Z are each present, unless otherwise indicated.

[0250] The terms used herein for the location of elements (e.g., "top", "bottom", "upper", "lower") are used only to describe the orientation of the various elements in the figures. It should be noted that the orientation of the various elements may vary according to other exemplary embodiments, and such variations are intended to be encompassed by the present disclosure.

[0251] The embodiments described herein have been described with reference to the drawings, which show certain details of particular embodiments implementing the systems, methods, and programs described herein, however, describing the embodiments with the drawings should not be construed as imposing limitations on the present disclosure that may be present in the drawings.

[0252] It is important to note that the construction and arrangement of the devices, assemblies, and steps as shown in the various exemplary embodiments are merely illustrative. Additionally, any element disclosed in one embodiment may also be incorporated or utilized in any other embodiment disclosed herein. Although only one example of an element from one embodiment that may be incorporated or utilized in another embodiment is described above, it should be understood that other elements of the various embodiments may be incorporated or utilized in any of the other embodiments disclosed herein.

[0253] The above description of the embodiments has been presented for purposes of illustration and description. It is not intended to be exhaustive or to limit the disclosure to the precise form disclosed, and modifications and variations are possible in light of the above teachings or may be acquired from the present disclosure. The embodiments have been chosen and described in order to explain the principles of the present disclosure and its practical application, so that those skilled in the art can utilize the various embodiments and make various modifications as appropriate for the particular use contemplated. Other substitutions, modifications, changes, and omissions can be made in the design, operating conditions, and arrangement of the embodiments without departing from the scope of the present disclosure, as set forth in the appended claims.

Claims

1. 1. A method comprising: acquiring, by a computing system, magnetic resonance imaging (MRI) data corresponding to a plurality of MRI scan images of an anatomical region including a tumor of a patient, the plurality of MRI scan images including a first image set obtained through a first scan performed prior to administration of a therapy to the patient and a second image set obtained through a second scan performed after the administration of the therapy to the patient, the therapy including administration of a plurality of drugs; determining, by the computing system, tissue characteristics of tissue surrounding the tumor from the MRI data; registering, by the computing system, image-related data generated from the second set of images to image-related data generated from the first set of images; determining, by the computing system, a diffusion characteristic of the tumor based on the tissue characteristics; determining, by the computing system, a growth characteristic of the tumor based on the tissue characteristics; determining, by the computing system, for each drug of the plurality of drugs, an effect of the drug on tumor cells; generating, by the computing system, a score indicative of the tumor's predicted response to the therapy based on the diffusion characteristics of the tumor, the growth characteristics of the tumor, and the determined effects of each drug on cells contained in each voxel.

2. The method of claim 1 , further comprising performing tumor segmentation to identify a tumor region of interest (ROI) based on the MRI data prior to determining the tissue characteristics.

3. The method of claim 1 , wherein the tissue characteristic is related to vascular structure within the tissue.

4. The method of claim 1 , wherein the MRI data includes dynamic contrast enhanced MRI (DCE-MRI) data, and the tissue characteristics are quantified based on the DCE-MRI data.

5. The method of claim 1 , wherein the MRI data includes diffusion-weighted MRI (DW-MRI), and the method further includes generating a map of the apparent diffusion coefficient (ADC) of water.

6. The method of claim 1 , further comprising determining drug distribution within each voxel of tissue.

7. 2. The method of claim 1, wherein the diffusion characteristics correspond to diffusion of the tumor cells mechanically linked to material properties of the tissue via a physical stressor, thereby representing changes in the tumor that can cause deformation within the tissue.

8. The method of claim 1 , wherein the growth characteristics are based on carrying capacity, which is related to the maximum number of tumor cells that can physically fit within a voxel.

9. The method of claim 1 , wherein the growth characteristics are based on a voxel-by-voxel growth rate, the growth rate being calibrated for each voxel within the tumor ROI of the patient.

10. The method of claim 1 , wherein the effect of each drug on tumor cells corresponds to the spatiotemporal distribution of each drug within the tissue.

11. 2. The method of claim 1, wherein the effect of each drug on tumor cells is based on at least one of an efficacy parameter α of the drug, a washout parameter β of the drug over time after each administration, or an initial concentration of the drug.

12. 10. The method of claim 1, further comprising determining a modified therapy based on the score indicative of the predicted response of the tumor to the therapy.

13. 10. The method of claim 1, further comprising administering a modified therapy based on the score indicative of the predicted response of the tumor to the therapy.

14. 1. A computing system comprising one or more processors and a computer readable memory having instructions, the instructions causing the one or more processors to: acquiring MRI data corresponding to a plurality of MRI scan images of an anatomical region including a tumor, the plurality of MRI scan images including a first image set obtained through a first scan performed prior to administration of a therapy to a patient and a second image set obtained through a second scan performed after the administration of the therapy to the patient, the therapy including administration of a plurality of drugs; determining one or more tissue characteristics of tissue surrounding the tumor from the MRI data; registering image-related data generated from the second set of images to image-related data generated from the first set of images; determining a diffusion characteristic of the tumor based on the tissue characteristics; and determining a growth characteristic of the tumor based on the tissue characteristics; and determining, for each drug of the plurality of drugs, an effect of the drug on tumor cells; and generating a score indicative of a predicted response of the tumor to the therapy based on the diffusion characteristics of the tumor, the growth characteristics of the tumor, and the determined effect of each drug on cells contained in each voxel.

15. 15. The computing system of claim 14, wherein the instructions are further configured to cause the one or more processors to perform tumor segmentation to identify a tumor ROI based on the MRI data prior to determining the tissue characteristics.

16. The computing system of claim 14 , wherein the tissue characteristic relates to a vascular structure within the tissue.

17. 27. The computing system of claim 26, wherein the MRI data includes diffusion weighted MRI (DW-MRI), and the instructions are further configured to cause the one or more processors to generate a map of an apparent diffusion coefficient (ADC) of water.

18. The computing system of claim 14 , wherein the growth characteristics are based on a carrying capacity related to the maximum number of tumor cells that can physically fit within a voxel.

19. 15. The computing system of claim 14, wherein the effect of each drug on tumor cells is based on at least one of an efficacy parameter α of the drug, a washout parameter β of the drug over time after each administration, or an initial concentration of the drug.

20. The computing system of claim 14 , further comprising determining a modified therapy based on the score indicative of the predicted response of the tumor to the therapy.

Images