Solving ODEs vs VoronoiMesh

I’m using VoronoiMesh to model cells in an epithelium. The concentration of certain proteins within the cells are determined by solving a system of ODEs per cell. My goal is to model emerging patterns based on the concentrations of such proteins in each cell, where the colour of each cell is determined by the concentration of a specific protein.

I have two questions regarding the implementation of this simulation:

  1. Is it possible to change the colour of each cell in a Voronoi mesh, given, for example, an array of the solutions of such concentrations? I’m aware of the MeshPrimitives command, but I’m not sure how to efficiently use it.

  2. Also, these concentrations depend on the protein concentrations of neighbouring cells (these are cell-cell signalling dynamics). Is it possible to detect neighbouring cells so that I can incorporate their protein concentrations in the ODE system for a specific cell? Say, something like NeighbourQ? Maybe changing the whole thing into a graph could help, but I’d appreciate any suggestions.

Any ideas? Thank you.