Logic in Activating DCs

[drbergman]

To me, there appears to be an implied relationship between is_Dirichlet, dirichlet_activation_vector, and dirichlet_activation_vectors. It seems dirichlet_activation_vectors is an $n_\text{voxels}\times n_\text{substrates}$ array and the others are (supposed to be) the result of any along the two dimensions.

If this is correct, then current implementations in BioFVM_microenvironment.cpp do not maintain this consistency. It can lead to unexpected behaviors in which a command such as update_dirichlet_node targeting a single substrate actually affects other substrates as well.

[drbergman]

As an example of unexpected behavior:

The implementation of update_dirichlet_node, which seems meant for being the simplest way of setting a DC for a single substrate at a single voxel, can lead to unexpected behavior. To see this, make the template project. Add DCs to substrate with value 38. Then, add a second substrate with DCs off. Finally, use the following in custom.cpp for setup_microenvironment and run.

void setup_microenvironment( void )
{
    initialize_microenvironment();
    microenvironment.update_dirichlet_node( 175 , 1 , 1 ); // note that this supposedly leaves "substrate" (with index 0) alone
    return; 
}

the first substrate will end up with a DC of value 0 also at voxel 175:

[MathCancer]

HI, I wanted to touch base with you. Is this still unresolved?

[drbergman]

Yeah, it still is unresolved when using the latest PhysiCell release.

[MathCancer]

Thanks, Daniel. Sorry I didn't see the response.

Maybe we should have a quick call about it in the next week or so, so we can examine and plan expected behavior for a fix in 1.11.1.

It might be a bug, or could be poor documentation, so I'd like to work it super carefully. Thank you!

[drbergman]

That would work for me!

[rheiland]

I created a unit test of sorts which seems to do what I expect. If nothing else, we have a concrete test for reference. I will likely do a PR to PhysiCell dev eventually; never too many tests. Rf the Readme.txt, the custom.cpp setup_microenvironment, and the 2 .png images. https://github.com/rheiland/PhysiCell/tree/development/unit_tests/custom_DCs_2substrates

[rheiland]

And I just created a custom_v2.cpp which is more similar to your example in that it only changes a single (center) voxel value of just one substrate (oxygen) and leaves the other alone. See the updated files and Readme.txt. I also fixed the config file so the PMB/Studio (v2.19.2) would be happy. See *v2.png files too.

[rheiland]

To more closely mimic Daniel's original test, I committed a new config file with the glucose BCs=99:

        <variable name="glucose" units="dimensionless" ID="1">
            <physical_parameter_set>
                <diffusion_coefficient units="micron^2/min">0</diffusion_coefficient>
                <decay_rate units="1/min">0</decay_rate>
            </physical_parameter_set>
            <initial_condition units="particles/micron^3">0</initial_condition>
            <Dirichlet_boundary_condition units="particles/micron^3" enabled="True">99</Dirichlet_boundary_condition>
            <Dirichlet_options>
                <boundary_value ID="xmin" enabled="True">99</boundary_value>
                <boundary_value ID="xmax" enabled="True">99</boundary_value>
                <boundary_value ID="ymin" enabled="True">99</boundary_value>
                <boundary_value ID="ymax" enabled="True">99</boundary_value>
                <boundary_value ID="zmin" enabled="False">0</boundary_value>
                <boundary_value ID="zmax" enabled="False">0</boundary_value>
            </Dirichlet_options>
        </variable>

Running the test using the simpler custom_v2.cpp where we only do:

void setup_microenvironment( void )
{
	initialize_microenvironment(); 	

    int idx_oxygen = 0;  // just hard-coded; better to get index by name.
    double oxy_value = 40.0;
    microenvironment.update_dirichlet_node( 4, idx_oxygen, oxy_value);
}

it seems like we still have the expected results? Perhaps I'm still missing something though. See https://github.com/rheiland/PhysiCell/blob/development/unit_tests/custom_DCs_2substrates/oxygen_v2.png and https://github.com/rheiland/PhysiCell/blob/development/unit_tests/custom_DCs_2substrates/glucose_v2.png

[rheiland]

I finally got back to this and think I've fixed the bug, taking an approach that Paul agreed with, whereby we assume DC activation is false everywhere (instead of true, as it was previously), then just set it true when needed.

Please test if you get a chance, from this branch: https://github.com/MathCancer/PhysiCell/tree/dev-randy-DCs

• see diffs in BioFVM_microenvironment.h and .cpp
  • changed true to false in many places
  • added a new method get_substrate_dirichlet_value
  • commented out a loop at the very end setting activation vectors (rf. April 2023 comment)
• see unit_tests/custom_DCs_2substrates (and Readme.txt there)

Not sure if we also want to fix a method with confusing args, rf. TODO? fix confusing swapped usage of args comment in BioFVM_microenvironment.cpp

[rheiland]

I just pushed this same fix into my (rheiland) fork (dev branch) for 1.12.0 (https://github.com/rheiland/PhysiCell/tree/development) and will do a PR into the MathCancer dev branch after some more testing. Please feel free to test your use case(s). I created a new, very simple one which also illustrated the bug before the fix:

• create 2 substrates: subA and subB
• specify DC only on xmax boundary of subA
• specify DC only on ymax boundary of subB
• run and plot results

Bug: Both DCs on both boundaries appear on both substrates.

[MathCancer]

Hello! Do we feel this is now resolved in 1.13.0? Thanks for your great work on this!

[drbergman]

Yes! I believe it is resolved!

[MathCancer]

@drbergman is a hero and this is long since resolved. thank you!