# ReadMe.md This ReadMe is for the data collection "1D Flow in Simulated Pulmonary Resection". The generation and analysis of this data is discussed in detail in our article "A Computational Framework for Pulmonary Assessing Wave Intensity Following Simulated Lung Resection" (https://doi.org/10.64898/2026.03.16.712097). We analyse pre-segmented patient-specific human pulmonary arterial trees in order to generate 1D computational domains in which to simulate fluid flow. To these computaional domains, we make systematic geometric changes by sequentially removing vessels to capture the change in pulmonary arterial geometry following lung resection. We simulate flow in all these derived altered and unaltered geometries. The flow simulations are saved and analysed. This data collection contains the original segmented pulmonary surfaces (in TreeData.zip) plus the derived centrelines for the vessels and surfaces. When unzipped TreeData contains 48 sub-folders most of which contain 5 files: the original stl surface, two centreline files (csv, vtp), and two surface files (csv, vtp). The centrelines are analysed using the scripts in TreeScripts to generate the 1D computational domains as described below (and in TreeScripts/TreeScriptsDocs.pdf). The 1D model code (1DPulmCode.zip) is called upon all of the generated computational domains (stored in SummaryTables.zip) to generate the simulated flows (FlowSimulationResults.zip). These are analysed using the contents of AnalysisScripts.zip to make the tables seen in the article. ## Tree Scripts 1. `GetInletIndices.m`: Run this to open all the STL files. While each is open, write down the index that appears on the MPA. These are needed in the next step. 2. `STL2Centerlines.m`: For each tree, create the centrelines that run between the inlet specified in the previous step. Each tree has its own folder. These are saved as `centrelines.vtp`. Use `DrawAllCentrelines.m` to draw the collection of centrelines to see if any look bad. 3. `STL2GroupSurface.m`: For each tree, create the group surfaces that constitute each branch. These are saved as `GroupSurface.vtp`. 4. `Centreline2Branches.m`: Open the centreline files with Centreline2Branches.m to make branches (from the inlet to each outlet for each network). Networks 10, 25, 35 do not work. Exclude these. Branches are saved in branches.mat. Each cell is a structure with attributes: - name, a string - branches, a collection of cells which contain matrices with columns: - radius, edge array 1, edge array 2, EdgePCoordArray, Points:0, Points:1, Points:2 We only need the first of these and the last 3 since the vertices appear in order. Branches are saved in `branches.mat`. From the branches, we can make segments using 5. `Branch2Segment.m`: Create segments from the branches. This has `MakeSegs.m` as a dependency. It creates segments given a collection of edges `E`, a list of terminal nodes `T`, a list of body nodes `B`, a list of junction nodes `J`, and an inlet `I`. 6. `reorient.m`: ensure the correct orientation of the segment away from the inlet. This has `NetworkStructure.m` as a dependency. Self loops and twins are removed in this step. There is now a cell structure saved as `RawSegments.mat` 7. For each tree in `RawSegments.mat` remove self-loops. We can look at the graphs of all the trees by calling `InspectGraphs(data)` (with `ChooseSize` as a dependency). These segments produce pictures as seen in Fig. 6(a) of the paper: the bifurcations are in the wrong place and the vessels are too short. 1. Call `AnalyseSegments.m` to analyse the segments by making any series joins, extending the segments and extrapolating the radius data, linearly interpolating the spatial trace of each segment to get a consistent step size between segments, and shifting junction points. Save the analysed segments as `AnalysedSegments.mat` 2. Visualise all of the trees to ensure that the MPA is segment{1}, LPA is segment{2}, and the RPA is segment{3}. This is currently done manually. It is important to do now so determining left and right vs operative and non-operative sides is easier later. 3. Now make alterations to the trees using `RemoveDownstreamVessels.m`. This has dependencies - `CategoriseJunctions.m` which determines the maximal junction order (bif-, trif-, …, n -furcation) in each tree. - `removeAdditionalVessels.m` which determines the vessels to remove for each altered tree. This outputs a cell structure with one cell per original tree. Inside each of these is a further collection of cells which give the vessels to be removed. The first cell is always empty: for the unaltered tree, no vessels are removed; cell 2 contains the indices of all the left vessels, and cell 3 contains the indices of all the right vessles. - `FindMainTrunk.m` - `MakeST.m` which makes summary tables Essentially, loop through the trees for which you have working boundary conditions, remove vessels, and create summary tables for each original and altered tree. The summary tables are found in the (compressed) folder called SummaryTables. Unzip before use. 4. Now call `RunInNetworks.m` on the collection of summary tables. Useful functions: `CylAroundSeg.m` draws the net of the vasculature and the collection of cylinders. `ChooseSize.m` gives a suggestion for the best layout for a subfigure layout. `InspectGraphs.m` can be called on a specific tree or all trees to draw the graphs. ## 1D Code The 1D model code constructs a network of 1D elastic walled blood vessels given a summary table. While the tree is being constructed, structured trees are generated if a vessel has `init` value 4 (column 4 of the summary table). The inlet vessel has `init` value of 1 which allows the inlet flow boundary condition to be applied. For each time step, we solve the model equations at all points in space in the computational domain given the boundary and initial conditions. The solver is a two-step Lax-Wendroff scheme. Junctions are matched by the conditions given in junctions.c. After the first period has run, we can make use of the convolution integral for the structured tree boundary conditions. The programme runs an integer number of periods until the solution in two successive periods converges or the stopping condition (based on number of iterations) is met. The code makes one `.2d` file per summary table in `.tsv` it was given. The names of these files match until the extension. The 2d files have columns: - vessel index - temporal point - spatial point - pressure - flow - area as per `PrintPxt` in `arteries.C`. ## Analysis scripts We have run the 1D computational model in all of the networks available. There is one output file per network file. The names match, so the output from simulating flow in the tree `human_001.7.tsv` is saved in `human_001.7.2d`. This tree is an alteration to the original tree `human_001.0.tsv` and the vessels listed in `to_rm{1}{8}` are removed. This cell structure is generated by calling `to_rm = removeAdditionalVessels(data);`. The offset from between the naming convention and the cell structure index arises because MATLAB indexes from 1 and C indexes from 0. The output 2d files have columns: - vessel index - temporal point - spatial point - pressure - flow - area as per `PrintPxt` in `arteries.C` For simplicity, we want to read all output files into MATLAB and make a cell structure. We do this with `WrangleOutput.m`. This script loops through the outputs, gets the tree and the alteration from the file name, first three vessels, and three unknowns ($p$, $q$, $A$) to make a nested cell strucutre called `input`. This has one cell per original tree. The name of each tree is saved in the cell in an attribute called `name`. The other attribute of each cell is the data. The length of `input{i}.data` for tree `input{i}.name` is `n+1` where `n` is the number of alterations to the original tree, and the `+1` is the original tree itself. Each cell of `input{i}.data` is itself a cell. The numbering from `to_rm` from above is preserved, so `input{i}.data{j}` is the data from tree simulation in which branches `to_rm{i}{j}` were removed from `to_rm{i}{1}`. There is one cell per alteration to the tree in addition to the original tree. Each alteration specific cell contains 3 cells. These are for $p$, $q$, $A$ in the MPA, LPA, and RPA, respectively. These data are from the spatial mid-point of the vessel. The data are saved in `SimulationOutput.mat` as `input` as this is the *input* for all following work. In order for the Fourier space boundary condition to work, the number of time-steps in the numerical scheme must be a factor of 2. The number of time points at which we save the simulation data must be some divisor of a power of 2 which is also a power of 2. The period of our simulations of $T = 0.7$s and we are aiming for a time step of 1ms or 701 timesteps per period. The closest power of 2 to 701 is 512. The data coming from the simulation have 512 time steps. The data are linearly interpolated and the Savitzky-Golay filter is applied in `smoothSimulatedData.m`. The data are interpolated first, and then the Savitzky-Golay filter is applied to three successive periods of flow $q$. The (odd valued) window width is chosen such that the SSE between the original signal and the filtered signal is less than 1 in the central period. This scheme preserves periodicity away from the end points. The same window size is applied to $p$ and $A$. This is repeated for all simulations. The central period is saved and has the desired resolution. Window size for each simulation is saved as $W$ in `WindowSize.mat` . Tables 1 & 2 are generated using the script `PreOpTable.m` which loads the data and calls `removeAdditionalVessels.m` to get the alterations. Wave intensity is computed for the MPA, LPA, and RPA in each tree (where present). This is acheived by calling `ComputeWaveIntensity` on the whole collection of simulation results. This has dependency `wave_intensity_calculation` which computes wave intensity for a single simulation.The time at which the period starts is shifted to align with the method described by Glass by calling `shift_t0`. The four wave components are computed by calling `GetWaveComponents.m` and from this, the parameters of interest are computed using `ParamsOfInterest.m`. Collect the pre-op values for each of the parameters of interest for each of the PAs. Generate the two tables. Table 4 is generated by running `LPAvsRPA.m` which works in essentially the same way as `PreOpTable.m`. Table 5 is generated by running PostOpTable.m. It, again, basically works the same way as the previous scripts: load the 1D flow simulation data, compute wave intensities and wave components to get parameters of interest, then slice the data to draw conclusions. Highlight the cells manually to show agreement with Glass.