Dynamic analysis of the individual patterns of intakes, voids, and bladder sensations reported in bladder diaries collected in the LURN study

Overview of the multistep analysis of BDs

First, we performed a traditional time-averaged analysis of the BDs. We examined the distributions and inter-subject correlations of demographic and time-averaged BD variables across the cohort of 197 participants with LUTS. Then, we introduced a dynamic approach to the analysis of BDs by examining intra-subject correlations between BD variables across the 3-day duration of the diaries, investigated cross-correlation functions and delays between these variables, and finally developed and evaluated dynamic models for individual BDs (Figure 1).

• Download figure
• Open in new tab

Figure 1. Flowchart of the multistep analysis of BDs.

Analysis steps are in rectangles, data and results are in the hexagons. Steps 1-2 constitute time-averaged analysis; steps 3-8 constitute a dynamic approach. First, we examined the distributions of BD variables across 197 participants. Second, we performed the time-averaged analysis of BDs by examining the matrix of inter-subject correlation coefficients of intake and voiding variables across 197 BDs. Third, we visualized individual BD by presenting graphs of drinking and voiding events and introduced dynamic variables describing drinking and voiding patterns. Fourth, we examined the intra­subject correlations between the introduced dynamic variables. Fifth, we calculated cross­correlation functions between the introduced dynamic variables to estimate both their similarities and time lags. Sixth, we developed and evaluated the individuals’ grey-box models predicting BFR using data on intake profiles and composition of the drinks. Seventh, we developed multivariable linear regression models for each individual urge growth rate (UrgR). Finally, we proposed behavioral modifications, i.e., restrictions of intake volume, caffeine, and alcohol consumption and simulated the effect of these modifications on UrgR and urge levels at voids.

Justification of dynamic approach

A dynamic approach to understanding reality, i.e., description of relationships between changes in variables in time rather than relationships between the values of the variables, proves to be productive in quantitative sciences, starting with Newton’s laws of motion, where acceleration, i.e., change in velocity of the object, was shown to be proportional to the force acting on the object. We propose to use a dynamic approach to the analysis of BDs, by exploring how changes in blood volume, osmolality, and concentrations of caffeine and alcohol, caused by fluid consumption affect (change) the rate of urine production by the kidneys. We propose to investigate how the sensation of urinary urge changes in time and how it is affected by the bladder filling rate and urine volume. Using this dynamic approach, we introduce urge growth rate (UrgR) as defined in the next section. We believe that UrgR often (with the exception of social or convenience voids) drives voiding behavior. For example, time interval between voids is determined by the time required to reach certain urge levels, and therefore is proportional to the reciprocal of the UrgR. In a sense, UrgR is more important than urge level at void. Any healthy person deprived of the opportunity to void will eventually get to a high level of urge, e.g., urge=3; however, it will likely take a healthy person much longer to get to this level than a person with overactive bladder (OAB). This makes UrgR a useful metric to evaluate the severity of OAB.

While implementing this dynamic approach to the analysis and modeling of BDs and to the simulation of potential behavioral modifications, we kept in mind that information available for modeling, as well as our understanding of the involved physiological processes, are incomplete. For instance, were have no information on physical activities or changes in heart rate that might affect the urine production rate by the kidney. We do not have food diaries and therefore do not have information on the fluids consumed in food (e.g., watermelon or soup). We do not know how much salt is consumed with food and when, which influences water retention. We do not know about the presence of psychological distractions and triggers that might influence the sensation of urge. All this incompleteness of information makes our task of predicting urgency episodes less like Newtonian mechanics calculation of the position of planets and more like the task of the captain determining the course of the ship having an approximate knowledge of the strength of the winds and currents, and therefore regularly measuring the position of the ship relative to the sun and the stars or using global positioning system (GPS) signals. Similarly, we try to minimize the propagation of the errors along the modeling process and, whenever possible, use the measured values of variables recorded in the BDs.

Dynamic variables

BDs report the values of intake variables (volumes and types of drinks) at times of intakes and voiding variables (volumes and urge levels) at times of voids, which do not coincide. Dynamic analysis requires the knowledge of intake variables and voiding variables at the same and preferably equidistant time points. To meet this requirement and to analyze dynamics of drinking and voiding patterns and bladder sensations, we defined several dynamic (time­dependent) variables across the whole duration of the 3-day BDs. We calculated rate of intake, bladder filling rate (BFR), urge growth rate (UrgR), and time-dependent frequency of voiding (F), as shown below. Stepwise approximations were used to derive these variables from the BD data. They were considered constant during the time intervals (durations of fluid intake and time intervals between voids), with changes occurring only at the boundaries of these intervals.

Fluid intake rate (IR(t)) was approximated as i^th intake volume (IV_i) divided by the estimated duration of i^th intake (Di) and assumed zero everywhere except during reported intake, where i=1,..M, M-number of reported intakes. No fluid intake was assumed or attributed to the food the participant ate; only fluids recorded on the BD were included in the analysis.

Frequency of voiding (F(t)) was defined as: where ΔT_j = T_j – T_j–1, time interval from previous void in minutes, T_j, j=1,..N – time of each of N voids reported by an individual. For easier interpretation, the 24 × 60 multiplier was introduced to make frequency defined by eq. 1 comparable with the averaged daily frequency, i.e., number of voids per day. Note that F defined by eq. 1 depends on time, unlike the frequency of voiding averaged across the day or several days typically used in the urologic literature [23,24]. For example, consider an individual who voided at 8:00 am, 9:00 am, 9:20 am, 9:40 am, 12:00 noon, 2:00 pm, 5:00 pm, 7:00 pm, 8:00 pm, and 11:00 pm. They would typically be described as voiding 10 times a day, or voiding every 90 minutes; however, there was a time during the day when they voided every 20 minutes, which is not reflected by the time-averaged definition of frequency, but is reflected by the peak (F(t=9:20 am)=72) of the time-dependent frequency (eq. 1). Note that the presence of such high-frequency voids at inopportune times could be much more bothersome for the individual than voiding 11 versus 10 times a day; therefore, it is desirable to have tools to predict such high-frequency episodes, as well as high-urgency episodes.

Urge was assumed to be zero immediately after each void, so UrgR(t) was defined as urge (U_j) reported at the j^th void divided by ΔT_j, assuming a linear growth of urge with time, which may not always be correct but is close to sigmoidal growth observed in several studies measuring real-time bladder sensations in individuals with and without OAB [25,26]. Even when the assumption of linear growth of urge is not valid, the variable UrgR(t)) defined as U_j/ΔT_j is of importance since it reflects how much time it takes to reach the certain level of urge. The time required to reach urgency is even more important than level of urge. For instance, an individual without OAB driving on the highway (with limited access to a toilet) might get a sensation of urge U_j =3 and need to exit the highway. However, an individual with OAB will get to urge U_j =3 much faster and would need to exit more often. Variable U_j/ΔT_j is different from frequency defined by eq. 1, since some of the voids could be social or convenience ones and occur at the low level of U_j.

The BFR, sometimes referred to as the rate of diuresis [27,28], was defined as: BFR(t) = VV_j/ΔT_j, where VV_j is the j^th voided volume (VV). Here, we assumed that the variation of PVR from void to void is small relative to the VV (justified in participants with PVR<50mL selected for this analysis); therefore, the VV can approximate the change of the urine stored in the bladder between voids.

Figure 2 illustrates the 3-day BD data (intake and VVs, and urge levels reported at voids) and its relationship with the dynamic variables derived from these data (BFR, time-dependent (F), and UrgR) for a typical study participant A. Supplemental Figures S1-S9 provide similar information for nine other representative participants (B-J). More detailed discussion of the figures is provided in the Results section; however, we think it is beneficial to introduce these figures here for better illustration of our dynamic approach to the BD analysis.

• Download figure
• Open in new tab

Figure 2. Intake and voiding profiles recorded in the 3-day BD of patient A and dynamic variables derived from the profiles.

Figure 2A: fluid intake volumes. Figure 2B: voided volumes (VV). Figure 2C: urge levels reported at voids. Data in Figures 2A-2C are presented in the form of stems instead of dots to illustrate the discrete nature of drinking and voiding events and to emphasize the impossibility of connecting the dots, e.g., intake volumes are equal to zero between intakes, and void volumes are equal to zero between voids. Figure 2D: bladder filling rate (BFR(t)). Figure 2E: time-dependent frequency of voiding (F(t)). Figure 2F: urge growth rate (UrgR(t)). Caffeine containing intakes are shown in Figure 2A in red. Note that peaks of BFR(t), F(t), and UrgR(t) are collocated with intakes of caffeine and with frequent and/or high-volume drinks. Some of these peaks are labeled to emphasize their collocation (X=time in hours from midnight of the first day of the diary; Y=values of the plotted variables). Note high (several fold) time variability of BFR.

Inter-subject and intra-subject correlations

Both inter-subject and intra-subject correlation coefficients of BD variables were calculated. The inter-subject correlations analysis was included in the time-averaged approach (step 2 in Figure 1). It was performed by first averaging the dynamic BD variables over time (3 days of the BD) and then calculating the matrix of the Pearson correlation coefficients of the time-averaged variables across 197 BDs. The intra-subject correlation analysis was part of dynamic approach (step 4 of Figure 1). It was performed by calculating matrices of Pearson correlation coefficients of dynamic variables for each of the individuals and then averaging the individual correlation matrices across all 197 individuals. Matrices of Pearson correlation coefficients were calculated using the MATLAB function corrcoef.m (MathWorks, MA).

Both for inter-subject and intra-subject correlations, we also calculated matrices of partial correlation coefficients, while controlling for certain variables. This type of analysis helps when there are three or more codependent variables (e.g., A, B, C), and we want to know how variable A responds to the changes in variable B, which are not caused by the changes in variable C. It is similar to controlling for variables in linear regression models [29] and was implemented using the MATLAB function partialcorr.m.

Cross-correlation function

Reactions of organisms to changes in external factors are not instantaneous but take time, which may be short or long depending on the physiological processes involved and may differ across individuals. For instance, it may take different time in different individuals for the BFR to react to the increase in the intake due to consumed fluid. Intra-subject correlations do not allow for capturing similarities in the presence of unknown time delays between the variables. To evaluate similarities and time lags between the functions representing the dynamic variables, we used cross-correlation function, which is routinely used to compare signals in signal processing [30] and is defined as: where f(t) and g(t) are two dynamic variables under comparison. The maximum possible value of cross-correlation function is X(τ_max) = 1, which happens when g(t + τ_max) = f(t). Therefore, the value of X(τ_max) represents the level of similarity, while τ_max represents the time lag between the two functions. We used cross-correlation function to compare the dynamic variables derived from the BD data. We also used cross-correlation function to compare BFR profiles predicted by the developed models and derived from the BDs. Cross-correlation function was implemented as MATLAB function xcorr.m.

Grey-box model of the urine formation rate

We created mathematical models for each individual participant, describing the dynamics of the BFR, which is equal to the rate of urine formation by the kidneys. In doing so, we followed a nonlinear grey-box modeling approach. Grey-box model is a type of model where the equations describing the processes are assumed known, while the parameters/coefficients of the equations have physical or physiological meaning and are determined by minimizing differences between observed and modeled processes. We described the model by a set of nonlinear ordinary differential equations (ODEs), with free parameters determined through minimization of the difference between urine formation rate profile predicted by the model and BFR profile derived from the BD. We assumed that the structure of the model was the same for all participants, while the parameters of the models were participant-specific. The core of our model is the modified model of Bighamian et al [31] describing the dynamics of the distribution of added water between blood and interstitial fluid, which is essential for determination of time lags between intakes and voids. We generalized this model by adding two more compartments (stomach and intestine) and by describing transport, kinetics, and redistribution of electrolytes, caffeine, and alcohol between blood and interstitial fluid. We also allowed for urine formation rate dependency on blood plasma volume, osmolality, and caffeine and alcohol concentrations. The model was implemented with System Identification Toolbox (MATLAB 2021a), using function idnlgrey.m. Details on the assumptions, equations, and implementation of the grey-box models of urine formation rate are presented in the Supplemental Material.

Comparing modeled and observed dynamic variables by evaluating locations and amplitudes of the peaks

The traditional metric to compare modeled and observed variables is called “fit percent” and is calculated as: where yobs(t_i) and ymod(t_i) are the values of observed and modeled dynamic variable y at all points t_i, where this variable was observed and modeled. Note that this metric compares differences in the observed and modeled values of the variable (numerator) with the overall time-variability of the observed variable (denominator), which makes it, in general, suitable for evaluating dynamic models. However, given the presence of high and narrow peaks in the variables of interest in the BDs (see BFR(t), F(t) and UrgR(t) in Figure 2 and Supplemental Figures S1-S9), which could be caused by high-volume or/and high caffeine content intakes, this metric might be unsatisfactory. For instance, it could happen that narrow peaks in yobs(t_{i obs}) and ymod(t_{i mod}) are closely collocated but not overlapping. In this case, the Fit Percent value would be lower than if the peak of interest was completely absent in the modeled function, which would misrepresent properties of the model. To avoid this situation, we introduced a new metric named Peak Fit, similar to Fit Percent, but allowing for some lags in the peaks of the modeled function; i.e., we identified ten highest observed and modeled peaks and searched for the modeled and observed maxima in the specified vicinity (30 minutes) of these peaks. The below equation for Peak Fit is similar to equation for Fit Percent but uses only the peak apex points compared with the maxima in the vicinity of the peaks. This way, it allows for limited lags in the modeled function relative to the observed function and punishes both the cases of missing peaks and of false modeled peaks, which are absent in the observed function. where t_{j obs peak} and t_{k mod peak} are positions of the apexes of the peaks in the observed and modeled functions, while t_{j mod max} and t_{k obs max} are the maximum values of modeled and observed functions in the vicinity of the above peaks.

To evaluate the average lag between observed and modeled peaks, we introduced another metric: Note that this metric is similar to the τ_max-lag provided by the cross-correlation function; however, the latter is dominated by the lag of the highest peak, while the former provides lag averaged across multiple peaks.

Evaluation and selection of grey-box models

An important quality of any model is its ability to describe/predict observations in the time range outside of the time range used for parameter fitting. To estimate this ability, we derived not only the models using 3 days of BD data (3-day model), but also the models using only first 2 days of BD data (2-day model). Then, we used the parameters of the 2-day models and 3-day intake data to predict BFR(t) during all 3 days. Details of the comparison of 2-day and 3-day models are provided in the Results section and are briefly summarized below. We evaluated the developed 197 individuals’ grey-box models of urine production rate by using the Peak Fit metric defined by eq. 4. We selected for further analysis and use only the individuals with Peak Fit>0.9 both for 3-day and 2-day grey box models of urine production rate. Models for 145 individuals satisfied these criteria. The mean value of Peak Fit for these models of BFR(t) was 0.94, and the mean value of Peak Lag was 6 minutes. We also compared the “full models”, including information on caffeine and alcohol content of the intakes, with the “no caffeine no alcohol” models ignoring this information. It was shown that “full models” provide substantially better fit.

Regression models of UrgR

Physiology of bladder sensations is less clear than that of urine formation by the kidneys; therefore, creation of the grey-box model of UrgR does not seem feasible. It is known, however, that feeling of urinary urge is initiated by stretching of the bladder wall [32]. We assumed that the reaction of the nervous system to wall stretching was nearly instantaneous, or at least happened at a much shorter time scale (seconds) than the time scale of the BD records (time intervals between drinks and voids – minutes) and the typical time required to establish equilibrium distribution of water between blood and interstitial fluid (minutes or hours). Given the assumption of the nearly instantaneous reaction of the nervous system to the input variables, we used the multiple linear regression approach to create personalized models of urinary UrgR for each participant. Analysis of the intra-subject correlations performed as described in the “Inter-subject and intra-subject correlations” subsection above and in the Results section below provided evidence of different levels of correlations of UrgR with BFR and with VV, which correspond to sensitivity, not only to the level of bladder wall stretch, but also to the velocity of stretching.

Given these observations, the inputs of the model included BFR (x₁=BFR(t)) and volume (x₂=V(t)) of the urine in the bladder (calculated as integral of the BFR from the time of previous void to the given moment plus PVR volume The last input (x₃) was binary and contained information on whether the participant was awake or asleep at the given moment, which allowed for possible bladder sensation differences in these two states. The formula for model specification can be presented as: where y = UrgR(t), while Scaling of the variables is not necessary in multivariable linear regression but was implemented for further comparison of the relative role of different input variables in predicting UrgR. For instance, intercept in this formulation represents the value of UrgR in the hypothetical case when all the input variables are kept constant during the 3-day BD and equal to the median values of these variables during these 3 days. The values of coefficients β₁, β₂ represent the increase in the UrgR in response to the change in the values of input variables equal to their medians. The above specification (eq 6) means that the model allows not only for linear dependence of UrgR(t) on the input terms, but also on their interactions. The values of UrgR(t) and BFR(t) were determined from the BDs and assumed constant during the time intervals between voids; asleep/awake data was also taken directly from the BDs. This method was chosen to describe the inputs of the model as accurately as possible. BFR (t) can be determined from the grey-box model, however, less accurately than from the BDs since the model does not take into account the changes in heart rate or the consumption of salted food or diuretics, which could contribute to changes in BFR observed from BDs.

Multivariable linear regression was performed using MATLAB function stepwiseglm.m, which uses forward and backward stepwise regression to determine the values of coefficients β_i and β_¡j. At each step, the function adds or removes terms to the model, based on the value of the ‘criterion’ chosen. In this situation, the default value of ‘criterion’ was used, i.e., the p-value for an F-test of the change in the deviance that results from adding or removing the term. The quality of the individual final models was evaluated by their adjusted R squared value (Adj. R²). Only the models with Adj. R² >0.3 were selected for further analysis and simulation.

Calculating urge level using the predicted UrgR: Predicting urgency episodes

Timings of voids are not completely determined by the intake information, since some of the voids are “social” or “convenience” voids, which occur without urge or at the low level of urge and cannot be predicted by renal physiology-based models. Urge sensation is supposed to be zero at the moments right after voids; therefore, timings of voids are necessary as independent data to determine the level of urge with UrgR predicted by the developed regression models. Therefore, we used BD data on the timings of voids (t_i) together with the predicted by the regression models UrgR (pUrgR(t)) to predict the values of urge levels at the times of voids: Note that pUrgR(t) predicted by the model (unlike UrgR estimated from the BD) is not a necessary constant in the interval between the adjacent voids; therefore, integration rather than multiplication by the time interval is warranted. We then used eq. 7 to calculate urge levels at the time of each void and evaluated percentage of true positive and false positive predictions of urge episodes by comparing with urge levels reported in BD. Note that, while observed/reported values of urge at the moments of void have discrete integer values 0,1,2,3,4, the predicted values of urge pU(t_i) can have any positive value. Given the previously-outlined definitions of urge levels [19], we counted all voids where U(t_i)=2,3, or 4 as urgency episodes. In comparison, for the predicted urgency episodes, we implemented a “personalized urgency threshold” approach by allowing the threshold level for each patient to change from 1.05 to 2 and selecting for each patient the threshold that minimized the sum of percentages of false negatives and false positives. This is equivalent to adding one more parameter to the individual’s model of sensation of urge.

Simulation of “what if’ scenarios and personalized behavioral modifications

We used the developed models to simulate personalized behavioral modifications aimed to avoid or minimize urgency episodes. Simulation is different from model development since it does not involve parameter fitting. Parameters of each individual model were already determined at the stage of model development. In simulation, parameters of the models are fixed, inputs of the models are modified, and the changes of the outputs caused by the modifications of the inputs are reported. Simulating modifications of intake patterns involve running with the modified intake data the developed individual grey-box models of urine production and the regression models of UrgR. In particular, some of the simulated modifications excluded use of alcohol and/or caffeinated drinks, substituting them with the alcohol-free and decaffeinated drinks with the same osmolality and volume. The BFR(t) profiles with proposed modified inputs and with the initial true inputs are calculated using the grey-box models; the difference in two profiles is determined and added to the BFR(t) profile determined from the BD. The result constitutes our best guess of the BFR(t) profile and serves as the first input variable for the developed multivariable linear regression model of UrgR(t), while the second input variable is determined by the time integral of the first one. Note that these simulations are different from “no caffeine no alcohol” models discussed above. While the comparison of “no caffeine no alcohol models” with “full models” served to demonstrate the importance of considering the effect of caffeine and alcohol on BFR for achieving adequate fit of modeled and observed profiles, the simulations with excluded caffeine and alcohol containing beverages served to examine if abstinence of these beverages could help to avoid urgency episodes.

Distributions of the BD variables across 197 participants

Our cohort of 197 participants with BDs suitable for dynamic analysis is quite heterogeneous. There is substantial variability across the cohort by age, body mass index (BMI), 3-day intake volumes, mean intake volumes, types of drinks (osmolality, caffeine and alcohol content), voided volumes (VVs), bladder capacity (estimated as maximum VV), urge levels, and numbers of leaks, as illustrated by Table 1 and histograms in Figure 3.

• Download figure
• Open in new tab

Figure 3. Distributions of the BD data and the derived dynamic variables across the 197 participants.

Rows 1-2: intake properties. First row: 3-day intake volume (mL), mean intake volume (mL), standard deviation of intake volume (mL), ratio of maximum hourly intake volume, and total intake volume. Second row: mean osmolality of consumed drinks (mOsm), standard deviation of osmolality of the drinks, total consumed caffeine (mg), total consumed alcohol (mL). Rows 3-4: voids properties. Third row: 3-day void volume (mL), mean void volume (mL), maximum void volume, i.e., bladder capacity (mL), number of voids per day. Fourth row: log2 (max BFR(t)/median (BFR(t))) – measure of BFR versus time variability, mean urge, maximum urge, 3-day number of leaks. Note that our cohort is quite heterogeneous, with 3-day intake volume as low as 2L for some patients and as high as 14L for others, osmolality of the drinks varying from 0 to 1200 mOsm, both within the BD for a given participant and across participants. Bladder capacity varied from 200mL to 1400 mL, number of voids per day from 3 to 19, ratio of maximum to median BFR(t) from 1.6 to 128, mean urge from 0 to 4, and number of leaks during 3 days from zero to 33.

Intake volume variability during the day differed substantially across participants, with the ratio of maximum hourly intake volume to total 3-day intake volume ranging from 7% to 30%. Time-variability of void volumes during the day, measured as the ratio of mean to maximum void volume ranges from 0.18 to 0.88, with cohort median equal 0.48 and 75 percentile equal 0.58, supporting the observation [33] that the majority of people void when their bladder is at most 50% full. The ratio of maximum to median BFR for a given participant varies from a rather low 1.6 to extremely high value 128. Specifically, max/median BFR >4 in 70% of participants and >8 in 36% of participants, demonstrating that the high variability of BFR observed in Figure 2 and Supplemental Figures S1-S9 is a general feature of the cohort.

Table 2 provides values of significant correlation coefficients between the variables from Table1, sorted from stronger to weaker, and their partial correlation coefficients calculated with the rest of the variables kept constant.

Most of the inter-subject correlations are either trivial (e.g., correlations between 3-day intake and void volumes, number of leaks and mean urgency during 3 days, or mean daily frequency and 3-day void volume), or weak, or become weak when controlled for other variables, therefore not allowing for identification of definitive predictors of frequency, urgency, and incontinence episodes in our heterogeneous cohort and demonstrating the need for a dynamic approach to the analysis of individual BDs.

More interesting information is provided by the analysis of the distribution of VVs divided by the individual bladder capacity (estimated as the maximum VV during the 3-day BD). First, we performed analysis similar to [33] by creating histograms of the number of voids with the given ratio of VV to the individual bladder capacity. Then, we focused the analysis by looking only at the “very urgent” voids, i.e., voids with U_j >=3. Finally, we investigated how many individuals had low-volume voids and low-volume “very urgent” voids. Results of this analysis are presented in Figure 4. Low-volume voids appeared to be quite prevalent: 134 of 197 (68%) individuals had low-volume voids, while 52 of 197 (26%) individuals had low-volume “very urgent” voids, i.e., voided with urge U_j >=3, when only 20% or less of their bladder volume were filled with urine. While for the majority of individuals, low-volume “very urgent” voids were infrequent (20% or less of their voids), for some individuals, such voids happened as often as in 50% of cases. This analysis indicates that, at least for some of the individuals, the fraction of bladder volume filled with urine is not the only or even the main determinant of urinary urge. Next, we performed dynamic analysis of individual BDs in search of such determinants.

• Download figure
• Open in new tab

Figure 4. Distributions of the ratio of the voided volumes (VVs) to individual bladder capacity.

Figure 4A: number of voids with the given ratio of VV to the individual bladder capacity in the cohort of 197 individuals. Figure 4B: number of individuals with low-volume voids (ratio<=20%) versus the number of such voids in the 3-day BD. Figure 4C: number of individuals with low-volume voids (ratio<=20%) versus the fraction of such voids in the 3-day BD. Figure 4D: number of “very urgent” voids (U_j >=3) with the given ratio of VV to individual bladder capacity in the cohort of 197 individuals. Figure 4E: number of individuals with low-volume “very urgent” voids (ratio<=20%) versus the number of such voids in the 3-day BD.

Figure 4F: number of individuals with low-volume “very urgent” voids (ratio<=20%) versus the fraction of such voids in the 3-day BD.

Intake and voiding profiles

Three main observations can be made from visualization of the individual intake and voiding profiles and dynamic variables derived from the BDs. Here, we demonstrate it for a typical participant A (Figure 2) and for nine other representative participants B to J (Supplemental Figures S1-S9). First, dynamic variables, i.e., bladder-filling rate (BFR), instantaneous time-dependent frequency of voiding (F), and urge growth rate (UrgR) can vary several-fold (e.g., 12-fold [Figure 2], 16-fold [Supplemental Figures S2, S3, S5, S7, S9]) during the day. Second, peaks of these functions occur at time points collocated with fluid intakes. The highest peaks of BFR(t), F(t), and UrgR(t) are closely following consumption of caffeine and/or alcohol, e.g., in participant A (Figure 2), and in participants B, C, D, and I (Supplemental Figures S2, S3, S4, S9). Other peaks are collocated with frequent and/or high-volume drinks. Third, periods of low BFR(t), F(t), and UrgR(t) are collocated with time intervals with no fluid consumption, usually at nighttime.

Next steps serve to formalize and quantify these observations by examining intra-subject correlations and cross-correlation functions of the above dynamic BD variables, and by developing personalized models of BDs.

Intra-subject correlations of dynamic BD variables

In this section, we will show how correlation analysis of time-dependent (dynamic) variables reveals a powerful correlation of BFR with UrgR and F. We will also show, in contrast, a lack of significant correlation of urge level at void with either BFR or VV. BD data for a typical participant A (Figure 2) can be presented in a different way by plotting UrgR versus BFR (Figure 5A) and versus VV (Figure 5B). Similarly, F can be plotted versus BFR (Figure 5C) and versus VV (Figure 5D). In this presentation, information regarding the sequence of events is lost; however, presence or absence of correlations between variables becomes more visible. For instance, for participant A, BFR is strongly and significantly correlated with UrgR (R=0.88, p= 3 · 10^-11) and with F (R=0.85, p= 2 · 10^-10), while correlation of these variables with VV is weak and unsignificant (R=-0.19, p=0.29; R=-0.30, p=0.08). Correlations of urge level at the time of void (URG) is weak and unsignificant with VV (R=0.32, p=0.06) and BFR (R=0.31, p=0.08) (Figures 5E-5F). As shown in Figure 6, these observations hold for the majority of 197 patients. While correlations of UrgR and F with BFR are almost always positive and are above 0.5 for the majority of participants, correlations of these variables with VV can be both positive and negative and are almost all weak (|R|<0.5). Figure 6C presents the correlations of URG with VV and BFR. Unlike UrgR and F, where correlation with BFR was a strong positive for almost all participants, there are participants with positive, zero, and negative correlations of URG with BFR. Similarly, there are participants with positive, zero, and negative correlations of URG with VV. Four extreme cases with the strongest positive and negative correlations are presented in Table 3. Note that the strongest negative correlation of URG with VV is observed for participant Y, consuming high osmolality drinks as well as caffeinated and alcohol-containing drinks. Negative correlations with BFR and VV are observed in participant Z, consuming the highest amount of caffeine. Likely, for these participants, urge level is determined not by the urine volume but by the composition of the drinks. The strongest positive correlation of URG with VV is observed for participant W, who consumed low osmolality drinks without caffeine and alcohol. The strongest positive correlation of URG with BFR was observed for participant X, with the highest time variability of BFR (max/median BFR=95.8). Observed differences in the correlations of URG with BFR and VV across participants suggest the differences in mechanism of urge, which could be potentially used for subtyping of patients with LUTS.

• Download figure
• Open in new tab

Figure 5. Example for a typical participant A of the correlations between the dynamic BD variables.

Figure 5A: UrgR versus BFR (R=0.88, p= 3 · 10^-11). Figure 5B: UrgR versus VV (R=-0.19, p=0.29). Figure 5C: F versus BFR (R=0.85, p= 2 · 10^-10). Figure 5D: F versus VV (R=-0.30, p=0.08). Figure 5E: URG versus BFR (R=0.31, p=0.08). Figure 5F: URG versus VV (R=0.32, p=0.06).

• Download figure
• Open in new tab

Figure 6. Scatter plots representing the values of correlation coefficients between dynamic variables for 197 participants.

Figure 6A: correlation of UrgR with BFR (mean R=0.73) versus correlation of UrgR with VV (mean R=-0.07). Figure 6B: correlation of F with BFR (mean R=0.75) versus correlation of F with VV (mean R=-0.17). Figure 6C: correlation of URG with BFR (mean R=0.12) versus correlation of URG with VV (mean R=0.36). Four cases (W, X, Y, and Z with the highest and lowest values of R) are labeled and described in Table 3.

Table 4 presents the values of intra-subject correlation coefficients between dynamic BD variables calculated for each participant and then averaged across 197 participants. As seen, BFR, F, and UrgR are strongly and significantly correlated (correlations are preserved when controlled for other variables; see values of partial correlation coefficients); correlations of these variables with VV are weak, indicating that BFR rather than VV is the main determinant of urinary frequency and UrgR. This conclusion corroborates the observation of numerous low-volume “urgent” voids illustrated by Figure 4. There are two aspects of the correlation of UrgR with BFR. The first is obvious: A higher BFR results in faster bladder-filling to a level that produces a sensation of urge and a greater frequency of urination. The second aspect is revealed when we compare the full correlation coefficients with the partial correlation coefficients calculated when certain variables are fixed. For instance, partial correlation of BFR with UrgR (R=0.758) with fixed VV is even higher than full correlation of BFR with UrgR (R=0.733), while partial correlation of BFR with UrgR with fixed urinary frequency is relatively weak (R=0.3, p-value<0.001), indicating that the role of high BFR cannot be reduced to the time necessary to fill the bladder to a certain level. We further investigate the role of BFR and urine volume in the regression models of UrgR.

Investigating similarities between BD variables using cross-correlation functions

Intra-subject correlation analysis in the previous section helped to reveal important relationships between variables in individual BDs; however, in such analysis, information about the sequence of intake and voiding events is lost. To preserve this information and analyze similarities in profiles of the variables, including co-occurrence of events, we used cross­correlation functions, calculated as explained in the Methods section.

Cross-correlation functions of dynamic BD variables for a typical participant A are presented in Figure 7 and for nine other typical participants (B-J) in Supplemental Figure S10. The maxima of the cross-correlation functions of BFR(t), F(t), and UrgR(t) for participant A occurred at zero lag and have high values of over 0.96 (highest possible value of cross­correlation function is 1), demonstrating the high level of similarity of the profiles and absence of time delay between these functions. This result corroborates the observation from Figure 2 that peaks of BFR(t), F(t), and UrgR(t) are co-occurring (e.g., labeled peaks at X=30, Figure 2). It also corroborates the result of the previous section that urinary frequency and UrgR are highly correlated with the BFR. The presence of two smaller maxima of these correlation functions at ±24 hours indicate the similarity of voiding patterns in days 1,2,3 of the BD, i.e., stability of the daily voiding pattern for this participant. Note that peaks of cross-correlation function at ±24 hours occur for some (B, D, G, H, I, J) but not for all participants in Supplemental Figure S10, indicating that similarity of voiding patterns of days 1,2,3 of the BD is not a universal, but a patient-specific property.

• Download figure
• Open in new tab

Figure 7. Cross-correlation functions between dynamic BD variables of typical participant A.

Cross-correlation functions of BFR(t), UrgR(t), and F(t) have high maxima (Xcorrmax=0.96) at zero lag (X=0), indicating high similarity and absence of time delay between these functions. Two more maxima of these cross-correlation functions (X=±24 hour) indicate similarity of voiding patterns in days 1,2,3 of BD. Lower values of maxima of cross-correlation functions of frequency and UrgR with intake rate IR(t) (X=-0.05, Xcorrmax=0.34) indicate that peaks in intake volumes alone cannot predict peaks in voiding variables.

Cross-correlation functions of urinary frequency and UrgR with IR(t) demonstrate much lower values with maxima of about 0.34, indicating that peaks in intake volume alone cannot predict peaks in voiding variables. This corroborates with observations in Figure 2, where the highest peak (at X=30) co-occurred not with the highest intake volume but with the intake containing caffeine. Similar cross-correlation functions can be observed in Supplemental Figure S10 for the nine representative participants (B-J). For all of these participants, the maxima of cross-correlations of F and UrgR with BFR are in the range of [0.88, 0.98], which is three- to four-fold higher than their correlations with the IR(t). This observation makes sense, since IR(t) captures information only of the volumes and durations of the drinks, while urine production rate and therefore BFR may be affected by the composition of the drinks, e.g., osmolality, caffeine, and alcohol content. Also, importantly, BFR is affected not only by the volume of fluid consumed at the most recent intake but at the previous intakes as well, as seen for the peak of BFR (at X=68.3) following several peaks of intake (around X=66.8) in Figure 2, or similarly, for the peaks of BFR and UrgR (at X=14.23) following a series of intakes (at X=13.5) in Supplemental Figure S2. Therefore, analysis with cross-correlation functions confirmed conclusions of the analysis of intra-subject correlations on high correlations and co-occurrence of peaks of urinary frequency and UrgR with the BFR, but also demonstrated that more sophisticated methods of analysis are necessary to predict these peaks from the intake profiles.

Predicting BFR profile with the grey-box model of urine formation rate

As described in the Methods section and in Supplemental Material, we developed renal physiology-based grey-box models of urine formation rate for 197 individuals by fitting parameters of the model equations to minimize the difference of the urine formation rate predicted by the models and BFR calculated from the BDs of these individuals.

Figure 8 presents the comparison of BFR profiles derived from the BD with results of the grey-box models of urine formation rate for participants A (Figure 8A), B (Figure 8B), D (Figure 8C), and J (Figure 8D). Two versions of the model are presented, i.e., the “full model” and the model ignoring caffeine and alcohol contents of drinks (“no caf no alc”). As seen, the “full model” is capable of better describing the peaks of BFR especially collocated with intakes containing caffeine. The “no caf no alc” model predicts peaks collocated with intakes; however, these peaks are lower and broader than observed and predicted by the full model, as seen for participants A and C (Figures 7A, 8C). Figures 8B and 8D demonstrate that models ignoring caffeine and alcohol content of the drinks not only have problems in accurately predicting high peaks collocated with caffeine and alcohol consumption, but also sometimes create high false peaks collocated with drinks without caffeine and alcohol (e.g., around 7 pm on the first and third days in Figure 8B and on multiple locations in Figure 8D). This is the consequence of “no caf no alc” models not differentiating between different types of drinks and therefore attributing higher effects to the volume of the drinks. Although not ideal, “full models” are capable of better prediction of BFR profiles, with fewer peaks missed and fewer false peaks created. When comparing predicted and derived (from the BD) BFR profiles, it is important to remember that the later profiles are approximations as well, based on the assumption of BFR being constant during the interval between the voids, which is not necessarily correct and might be affected by multiple factors not recorded in the BD, e.g., heart rate, level of physical activity, consumption of salted food and/or diuretics.

• Download figure
• Open in new tab

Figure 8. Comparison of BFR profiles derived from the BDs with the profiles predicted by the grey-box models of urine formation rate.

Figure 8A: participant A. Figure 8B: participant B. Figure 8C: participant D. Figure 8D: participant J. Both “full models” (in red) and models ignoring caffeine and alcohol content of the drinks (in green) are presented.

Evaluation and selection of the BFR models

We evaluated the developed grey-box models by comparing BFR profiles predicted by the models and derived from the BD. We performed not only visual comparison (Figure 8) but also quantified similarities by using cross-correlation function (eq 1) and functions Peak Fit and Peak Lag (eqs. 3-4). As noted in the Methods section, an important quality of any model is its ability to describe/predict observations in the time range outside of the time range used for parameter fitting. To estimate this ability, we developed models using 3 days of BD data (3-day model), as well as models using only the first 2 days of BD data (2-day model). Then we used the parameters of the 2-day model and 3-day intake data to predict BFR(t) during all 3 days.

The two rows in Figure 9 serve to compare the similarity of the observed BFR(t) and the 3-days (Figure 9A-9C) versus 2-day models (Figure 9D-9F). As seen, the histograms are rather similar, indicating that 2-day model predictions are only slightly less accurate than 3-day models. For the majority of the 3-day and 2-day models, maximum of the cross-correlation function (max X_corr) is above 0.8 (Figures 9A, 9D). For the majority of the models, Peak Fit>0.9 (Figures 9B, 9E) and Peak Lag <6.5 minutes (Figures 9C, 9F). We used the Peak Fit metric to select satisfactory models for further analysis. Specifically, we selected BDs for whom both 3-day and 2-day models have similar accuracy, i.e., Peak Fit >0.9 for 3-day and 2-day models. We found 145 such BDs. Table 5 provides the comparison of the 3-day and 2-day models for these selected BDs by showing mean values of maximum and lag of cross-correlation function, as well as mean values of Peak Fit and Peak Lag metrics. As seen, both 3-day and 2-day models provide satisfactory prediction of peak heights and positions, and the 2-day models are only slightly less accurate than the 3-day models.

• Download figure
• Open in new tab

Figure 9. Distributions of metrics of similarity between modeled and observed (derived from BD) profiles of BFR for 197 participants.

Figures 9A-9C: 3-day models. Figures 9D-9F: 2-day models. Figures 9A, 9D: distribution of maximum values of cross-correlation function. Figures 9B, 9E: distribution of Peak Fit values. Figures 9C, 9F: distribution of Peak Lag values.

Predictions of linear regression models of urinary UrgR

As described in the Methods section, we developed multivariable linear regression models of UrgR for each individual with the following three inputs: BFR, volume of urine in the bladder, and awake (yes/no or 1/0). Figure 10 provides the comparison of profiles of UrgR(t) derived from the BDs with those predicted by the developed multivariable linear regression models for four typical participants (A, B, E, and F).

• Download figure
• Open in new tab

Figure 10. Comparison of profiles of UrgR predicted by the models and derived from the BDs for four typical participants.

Figure 10A: participant A (Adj.R² =0.88). Figure 10B: participant B (Adj.R² =0.87). Figure 10C: participant E (Adj.R² =0.86). Figure 10 D: participant F (Adj.R² =0.77). Note that intake and voiding profiles for these participants are presented in Figure 2 and Supplemental Figures S1, S5, and S6.

The quality of the multivariable linear regression models (the proportion of variance explained by the model) is estimated by Adj.R². The distribution of Adj.R² for 145 participants is presented in Figure 11A (0.1< Adj.R²<0.95). There is no universal rule to judge the quality of the model by the value of Adj.R². It is recognized that higher values of Adj.R² are expected in physics and chemistry and lower in biology, medicine, and sociology. Typically, models of the processes involving human behavior are considered “good enough” if Adj.R²>0.3 [34]. We applied this threshold and selected for further analysis the 124 participants for whom Adj.R²>0.3 in the multivariable linear regression model of UrgR. These 124 participants are a subset of 145 participants for whom the satisfactory grey-box models (Peak Fit>0.9) of urine production rate were created. We used the predicted profiles of UrgR(t) to calculate predicted urge scores at the times of void by using eq. 7 and compared the urge scores predicted by the models and reported in BD, as described in Methods, to determine percentages of true positive and false positive predictions of urgency episodes during individual 3-day BDs. The distributions of true positives and false positives are presented in Figures 11B and 11C. For 87 (70%) of individuals, models correctly predicted more than 90% of urgency episodes. Similarly, models predicted less than 10% false positive episodes in 82 (66%) individuals. The mean percentage of true positive predictions across 124 participants is 84.6%, while the mean false positive is 10.9%.

• Download figure
• Open in new tab

Figure 11. Quality measures of the multivariable linear regression model of UrgR and urge scores at voids.

Figure 11A: Distribution of adjusted R squared values for 145 individuals’ UrgR model (mean Adj.R² =0.67). Figure 11B: Distribution of percentage of true positive predictions of urgency episodes in subset of 124 participants with Adj.R² >0.3. Figure 11C: Distribution of percentage of false positive predictions of urgency episodes in the same subset of 124 participants. There is a rather weak positive correlation between Adj.R² of the model and percentage of true positive predictions of urgency episodes (R=0.3, p=0.006), and very weak negative correlation between Adj.R² and false positive predictions (R=-0.17, p=0.054).

One goal of creating multivariable linear regression models was to investigate how BFR and urine volume affect the UrgR in various individuals. We compared the group of participants (n=106) experiencing at least one episode of urgency (URG≥2) with those without urgency (n=18). Table 6 provides the mean and standard deviation values of the linear model coefficients for these two groups.

Intercept is significantly higher in the group with urgency. Similarly, β₁ is positive in both groups but is significantly higher in the group with urgency, indicating that BFR is a strong driver of urgency. Mean of ß₂ is negative in both groups but has two-fold larger absolute value in the group with urgency. Mean negative value of the coefficient for ß₂ is in agreement with the negative mean intra-subject correlation of UrgR with urine volume (Table 4). The negative value of ß₂ might seem counterintuitive; one would expect urge to be higher with the higher urine volume when everything else is kept equal. However, it is important to remember that our model

(eq 6) predicts urge growth rate not urge level, meaning that with negative β₂, urge is growing with the increase of urine volume but slower than linearly, suggesting possible saturation of the signal produced by the stretched bladder wall. Mean of β₃ is negative in both groups but has nearly three-fold larger absolute value in the group with urgency, indicating that, when everything else is kept equal (i.e., BFR and urine volume), on average, UrgR is higher when participants are asleep than when they are awake in both groups, but to a larger extent in participants with urgency.

Distributions of the intercept and model coefficients for participants with and without urgency in Figure 12 corroborates this observation. It also demonstrates that distributions of intercept and all coefficients are broader for participants with urgency relative to those without urgency. Rephrasing Leo Tolstoy, one can say that all people without urgency are “non-urgent alike”, while each person with urgency has urgency in their own way. Note that, while intercept and β₁ are positive for all 124 participants (both with and without urgency), the values of intercept and β₁ are higher for participants with urgency; β₂ is negative for all participants without urgency and for the majority (107, or 92%) of participants with urgency. Coefficient β₃ is negative for the majority of participants (85 with urgency and 12 without urgency), but is positive for some participants (31 with urgency and 6 without urgency), indicating that being awake could increase or decrease urgency growth rate in participants from both groups.

• Download figure
• Open in new tab

Figure 12. Distributions of the intercept and 6 coefficients of multi linear regression models of UrgR across 124 participants.

Blue: participants with urgency (n=106). Red: participants without urgency (n=18). Bottom-right panel compares the distributions of Adj.R² values for participants with and without urgency (no significant difference).

There are several metrics to estimate the severity of urgency in a patient; one possibility is to count the number of urgent voids. Table 7 provides correlations between the number of urgent voids with coefficients of the linear regression models across 124 participants.

Number of urgent voids is strongly, significantly, and positively correlated with the intercept and coefficient ß1. Correlation with the intercept is obvious and trivial (it just means that participants with higher UrgR have higher number of urgent voids). Positive correlation with the coefficient β₁ is more interesting and indicates that, not only the BFR is driving urge, but the participants with more frequent urgent voids are more sensitive to the increased BFR. Number of urgent voids is also significantly negatively correlated with β₂, suggesting that participants with a high number of urgent voids are less sensitive to the increased volume of urine in the bladder relative to the increased BFR.

It is of interest to examine if the slopes (sensitivities to different factors) are correlated, e.g., are people more sensitive to the increase of BFR also more sensitive to the increase of urine volume or to being awake. Scatter plots of intercept and coefficients β₁ and β₂ in Figure 13 demonstrate that intercept and β₁ are positively correlated (R=0.47, p<0.001), meaning that the larger the mean UrgR for the participant, the higher their sensitivity to the increase in BFR. On the contrary, correlation between β₁ and β₂ is negative (R=-0.33, p<0.001), indicating that participants more sensitive to increased BFR are less sensitive to increased urine volume. Note that correlation between β₁ and β₂ is different from the interaction term β₁₂, which deals with interactions of BFR and urine volume within a model for each individual. Negative ß₁₂ indicates slower than linear growth of UrgR with BFR. Correlation between β₁ and β₂ are across the individuals. Negative correlation indicates that participants more sensitive to increased BFR are less sensitive to increased urine volume.

• Download figure
• Open in new tab

Figure 13. Scatter plots of intercept and two model coefficients for 124 multivariable linear regression models of UrgR.

Figure 13A: intercept versus β₁ (R=0.47, p<0.001). Figure 13B: ß2 versus ß1 (R=-0.33, p<0.001).

Of 21 correlations between 7 coefficients of the models, 12 are significant (Table 8) but the majority of them are weak to moderate. Of main interest are positive correlations between intercept and sensitivity to BFR (β₁), and negative correlation between intercept and ß₂ and β_3, indicating that participants with higher UrgR are more sensitive to BFR, less sensitive to urine volume, and have larger difference between day and night UrgRs. Correlations between β_i and β_ik, and β_ij and β_ki are stronger; for instance, strong negative correlation of ß₂ and ß₂₃ indicates that people more sensitive to urine volume are also more sensitive to it during the night than during the day. Overall, the complex structure of correlations between coefficients of the multivariable linear regression models of UrgR indicates the possibility of different mechanisms of urge in different people. It is not clear (due to a limited sample size, n=124) if it is possible to identify distinct subtypes of patients with urgency based on these data, or if it is a continuum of different levels of sensitivities to various factors. Subtyping of patients with urgency would benefit from inclusion of data outside of BDs (patient-reported symptoms, physical examination, etc.) and would require consideration of the correlations across these variables [35]. Such study is of interest but is beyond the scope of the current paper.

Comparison of participants with satisfactory and unsatisfactory models

As described above, models of urine production rate and UrgR for 124 (63% of 197 participants) satisfied the selection criteria and were used for further analysis and simulation. It is of interest to evaluate how these 124 participants differ from 73 participants for whom satisfactory models were not created. Table 9 presents the comparison of the mean values and ranges of the BD variables for these two groups of participants. Variables are sorted starting with those with the most significant differences between the groups. Participants for whom satisfactory models were not created demonstrated significantly higher maximum across 3-day values of time-dependent instantaneous urinary frequency F(t) (Max F). They consumed more alcohol and high osmolality drinks. This demonstrates that our models have some difficulties in predicting extremely high peaks of F and in predicting the effects of alcohol and high-osmolality drink consumption. On the other hand, there are no differences in total intake and void volumes, maximum urge, number of voids and leaks, and sex of the participants. Differences between the groups in caffeine content of the drinks and mean urge levels are borderline significant.