Ecological Release and Patch Geometry Can Cause Nonlinear Density-Area Relationships

Hi, everyone. Welcome to today's Wolfram guest speaker webinar. This is Jamie Peterson with Wolfram U. I'm pleased to welcome Professor Jerome Goddard II today to share with us about his research and recent paper on ecological release and patch geometry. Professor Goddard is a principal investigator of a cross disciplinary team at Auburn University Montgomery. He'll share presentation slides Wolfram notebook showing examples from his research. So, thank you again. Thank you, Jerome, for being here and thanks to everyone in the audience. I will hand things over to you.

All right, I appreciate it. Let me just start by saying thank you for this opportunity to speak. I'm an avid Mathematica user. I've been using it since early days of graduate school. It's been an integral part of my research. And so I was just really excited and honored to be invited to speak today. So let me turn the camera off and let's jump over to my slides. All right. First, let me mention this work today is a collaborative effort. My research team is funded by the National Science foundation, so we're very thankful for their support. On the math side, we have a math modeling group headed by Shibaji at uncg, by myself at aum. We also work with an ecologist, Jim Cronin, and his group out at lsu. Here's a couple of pictures of Jim in the field showing some of our math students how he collects insects and bugs and things like that. So this is again a collaborative effort. Let me start by defining what is a density area relationship. So basically that's just the relationship between area of a patch and population density. Why do we care? So DAR is critical to conservation if we only have a finite amount of area, amount of land for conservation. We got two questions. Do we put it all together in one big patch or maybe do we break it up into smaller patches? DAR plays a role here. If you have a negative relationship between density and area, that would mean larger areas have smaller density. If that's the case, we want to go with a lot of small patches. However, if we had a positive Darwin, we'd want to put all our resources into a single large patch. So understanding something about the DAR is important for conservation as well as other areas. DAR has been well described for many animal taxa. Interestingly, empirical studies report positive, negative as well as neutral DAR within species and even within the same studies. I mean, there have been studies where one particular species exhibited all three of these. We don't really have a good understanding of the mechanistic causes behind these differences in dar. The theory most often used to Explain. It would be the equilibrium theory of island biogeography, which should imply a neutral DAR resource concentration hypothesis, which should give you a positive dar, or density composition compensation hypothesis, which should give you a negative Darwin. But really there's no good understanding of why we would have intermediate levels of dar. Maybe somewhere in the middle is the right place. So let me introduce a little bit about the framework. We're going to use a theoretical population model to try to understand some of the mechanisms behind this DAR that we see in the field. So as ecologists, what we still don't fully understand is the relationship between persistence of a population in a patch, the patch size, the patch shape and the patch matrix. Now, as a mathematician myself, this matrix is an unfortunate word that's used here. To an ecologist, the matrix is the area around a patch. To a mathematician, the matrix is a rectangular array of numbers. But for this talk, anytime I say matrix, I'm gonna mean the area surrounding a patch. Well, this is way too much variability. Let's fix something. Let's freeze the shape. So how are we going to do that? Let's take a patch omega. Let's just say it's a two dimensional patch for now that has unit area. So it says say for example, 1km squared area. We take that patch omega. I'm then going to create another patch omega naught by scaling the patch by lx. So each of the spatial coordinates gets scaled by L. L is like patch. So now as I increase the patch size, I'm stretching the patch, or if I decrease L, I'm shrinking the patch, but I'm maintaining the same shape. What this does is it allows me to fix patch geometry and just explore patch size versus matrix hostility. Wonderful. So for today we'll be talking about two different kinds of models. The first is a single species model. Imagine we have an island surrounded by a hostile matrix. So you could think something like an island surrounded by the ocean at every time T and location X in the patch. We're going to model density U according to a standard reaction diffusion model. Notice here I have a absorbing boundary condition. No animals are able to live on the boundary. They're immediately killed or taken off and never to return to this particular island. In this model, D1 is the patch diffusion rate. R1 is the patch intrinsic growth rate. U naught is the location and number of animals at time zero. It's the initial density profile. R1F tilde is the per capita growth rate in the patch. For today's talk, I'm going to just look Exclusively at logistic growth. But in the paper we explored some other growth types. Here K1 is the patch carrying capacity given in terms of density. For the two species competition model, we still have an island, it's still surrounded by a hostile matrix. But now at every time T in location X, I'm modeling the density of two competing species, U and V. We follow a standard Lacta volterra diffusive competition model. Oh, that's a mouthful. So in this model V, I represent patch diffusion rate, ri, patch intrinsic growth rate, ki patch carrying capacity CI. Here is the interspecific competition strength. What does that mean? The C1, for example, measures how much competition U is sensing from V. In other words, the higher the level of V, the higher the density of V, the more it starves U for resources. So that's measured in this parameter C throughout here I will be 1, meaning u and 2 will mean v. Well, as a mathematician there's way too many parameters. We're gonna do a non dimensionalization. Don't worry so much about the math here. The beauty is it takes me to a model that's much more simple. So now my patch has been scaled down to have unit area. The patch area itself has been moved inside the differential equation in terms of a parameter called lambda. So what is lambda? It's got in it the growth rate, the diffusion rate and the area of the patch. So if you see lambda, just think area. I could also think of it as the patch area times a parameter R1, which is the growth to diffusion ratio. It's the ratio of growth in the patch versus diffusion in the patch as well. For the competition system we have the same lambda that appears, but now we have a few extra parameters. D naught's the ratio of the diffusion rates. R Naught's a ratio of growth rates. And the competition coefficients CI have been scaled by bi. Don't worry so much about this math. We're going to get to some more important things. But I got to have this in here to make the talk legal as math talk. Now, what are the steps for computing Dar? I'm going to give you some results here that we proved analytically. Some theorems we were able to prove about dar, but also want to show some computational results where we involve Wolfram and Mathematica to do the computations. So what are the steps in that case? Well, we would fix the model parameters. We would choose an average island geometry. What does that mean in nature? If you had a series of islands, they'd all have different shapes. So we can't really expand over all those different shapes. So we pick one shape that's kind of the average island geometry. We call that Omega. We assume the population density is near an equilibrium. We then vary the island area and compute the average density in our patch as predicted by the model. So the way we would do this is take the single species case. We define a function D here to represent density. Average density given an area. We simply integrate the density over the entire patch. We then divide by the area of the patch. That would give me the average density. But if you remember, our patch has unit area, so that goes away to one. We do a similar thing for a competition system. So for U, you have DU of A defined in this way and have a similar formula for V. Remember, U and V are in competition. We can then readily take these formulas and use Mathematica and come up with computational results for this density area relationship. Wonderful. Let me talk very briefly about the mathematics behind these equilibrium solutions. So when I go to the steady state model, my single species model becomes 5. This particular boundary value problem has a long history. This is really not a math talk per se, so I'm going to gloss over some of history. Just give a couple of references here. The key result that we need, that I'll recall from the literature, is if I have a logistic growth model, then 5 has a unique globally asymptotically stable positive solution. For all lambdas greater than this value, E1D. What the heck is E1D? E1D is a principal eigenvalue of Laplace's equation with absorbing boundary conditions. E1D is a number that quantifies the loss an animal would experience when they encounter that hostile matrix. Depending on the shape or the geometry of the patch, E1D is a crucial number in this analysis. Lambda basically represents patch size or area. So if the area is big enough, you'll have a positive steady state. You'll have a population. If that area is too small, the patch can't support a population. That's basically what the theorem says. For the competition model, when we go to the steady state, we end up with an extra parameter R. Here, R is the ratio of the growth to diffusion ratios, ratios floating around everywhere. Remember, R1 measures the ratio between the growth rate of U versus its diffusion rate. Similarly for R2 and V. Wonderful. In the literature, most of the time a model like 6 has been studied by varying the diffusion rates. But that's not really our interest today. Today we want to fix a patch size and vary the Area, I want to know how the density changes with respect to area or lambda. In our problem, I want to vary lambda. Interestingly, if you fix the diffusion rates and just vary lambda, you get some substantial differences compared to just varying diffusion rates. And so in a recent paper by my research group and nonlinear analysis, we showed lots of different things. But I'll summarize these in this one table. Don't get scared too much by the math. Let's just dive in and see what happens. What does this mean? Lambda, remember, is like patch area. If the area of an island is too small, coexistence between these two species is never possible. If the area of an island is big enough and the competition is weak, that means the inter specific competition less than the intra specific competition. If that's the case, coexistence always happens in the middle. If I have a semi strong case where the competition is really, really strong towards one species compared to the other, then I need a trade off. I need for example, V to have a smaller minimum patch size requirement and U to have less competition. That trade off between those two will allow coexistence. Similarly, I could reverse these two between the GD ratio and competition and get coexistence again. Alright, let's get into some results for the single species case with logistic growth. We were able to prove in this paper that if you have logistic growth, the single species model will always predict a strictly positive and continuous dar as long as the A, the area of the island is bigger than this crucial threshold where the E1D was that special principal eigenvalue and the R1 is the ratio of the growth rate to diffusion rate. Now of course we want to use Mathematica to get some hands on experience here to see what this DAR relationship would really look like. So let's take a convex patch, a disk, and let's jump over to Mathematica and run some code and see what the DAR actually looks like. So let me switch screens, I'm gonna switch over to code and this code has been shared with y'. All. Let me make sure this loads. Okay, good. There's a setup cell here that you'll run. I'll point out here that to calculate E1D, I'm using Mathematica's built in eigensystemsolver for PDEs. It's a wonderful function call here that lets us numerically calculate this E1D given any patch geometry that that we can represent in Mathematica. So we'll run that cell. All right, great. Now, for the disk, if I just use the disk command in Mathematica, I get a disk that has unit radius. But remember, I need unit area. The area of this disk is PI. So we're gonna rescale. Using the command region resize, I'm gonna rescale the disk so that the area becomes one. So that doesn't take much. Just a little bit of math here to figure out what you need to scale the radius by. And I'll end up with one. Now, with this disk, with unit area, I can compute the E1D value, which is 18.16. Now, where that really comes into play here is when you define the R1. Remember, that was the ratio of growth to diffusion. So if I assume the growth and diffusion rates are the same, the minimum patch size for an island would be 18 units. So 18 units of area could be kilometers squared, could be meters squared. Whatever your units are. If I let the intrinsic growth rate be higher than the diffusion rate, say 2, that lowers the minimum area required for the species to exist. If the diffusion rate is much higher. In other words, animals are moving out, encountering the hostile matrix and dying at a higher rate. I'm required to have a much larger patch size in order for that population to survive in the island. Next, we're going to use the command parametric ndsolve to solve my single species equation with dark lay boundary conditions. Remember, absorbing boundary. Because of the method we're using, we're going to specify a seeding value for the method. This seeding value is important because the zero or trivial solution is always there. And. And if the seeding value is too small, we might pick up that solution. So we're gonna seed the solver here with a value bigger than zero to try to stay away from the trivial solution. Let's run that command. Now, let's see visually what the density would look like in my patch. So I'm gonna run this code here in the next cell if this works. What we see here is the area of the patch. We see the patch as a disk. Right now, we have no density. There's no animals able to colonize the patch. Notice the area is 16, which is below the minimum patch size. As I get close to the minimum patch size and in fact cross it, even just barely, say 18.3, all of a sudden, almost miraculously, the patch comes alive and we're able to maintain a population. Now, the z axis here is density, and it's scaled according to carrying capacity. So this at most on the z axis here you get a one which will represent 100% of carrying capacity. Because of the non dimensionalization, we scale that parameter out so 0.01 would be like 1% of carrying capacity. It's very, very small. As I increase area, what's happening is you can't see it. The disk is not changing, but I'm really stretching the size of the patch. It's a bigger disk, it's got more and more area. As I increase the area, you'll notice the scale here is going up. We are able to maintain the population at higher and higher density levels. And in fact you can prove this and you can see it as well computationally. As the area goes to infinity, the density inside that patch approaches 100% of carrying capacity. So you can see that a little bit computationally here. As we increase the slider, you'll see that this density inside is approaching 100% of carrying capacity. So what we want to do is we want to create a curve where you give me an area, I take the density predicted by the model, I calculate the average density here and then plot that as a point area, comma, average density. So let's do that. I'm going to delete this manipulate for now so I don't mess up code later on. So here, just like in the slides, we have a function single species averagedensity ssad that will integrate the solution, the numerical solution of my problem over the patch. So let's go ahead and run that. For example, if I want the average density for a patch of Area 20, I would end up with this small percentage. It's a percentage here. So like five or five and a half percent of carrying capacity. Now the next code is going to generate our density area relationship. So I've got a list here, udar points that's going to hold points of the form area and average density of the patch. Remember I told you I needed to set the seed value a little bit higher to make sure that I don't end up with the trivial solution. For this example, we're going to let that growth to diffusion ratio be one. I'm going to do two passes. So I'm going to start with 0.1 for area, go up to 100.1 and do 50 steps. I'm also going to do a pass that goes from point one to a thousand with 50 steps. So 100 total steps. We're gonna calculate or mesh the area interval there in two different groups as one, that's the areas in step one as two, that's areas in Step two, join those together so that as is all the areas I wanna calculate the average density for. We then feed that set into the command parallel map which will evaluate slot here will pull out one area at a time from the list and check to see if we're below the minimum patch size. Remember we proved already from the literature if the area is below this level, there is no population. You can go ahead and zero it out. There's no need to calculate anything. If we're above this critical threshold, then we'll use our command here that we defined earlier to generate the average density. So we're going to run this command. While it's running, it takes roughly about 45 seconds on my laptop. My laptop's older generation, maybe an eighth generation intel with six cores. So it's going to take probably about 45 seconds. While it's running, I'll go ahead and talk through the rest of the code. I'm going to prepend to this list of points, our bifurcation point. So the critical patch size, area here and 0. Well then use the sort command to sort them by the first component, which is area. So as it's comparing each of those points in the list, it will pull out the first component and sort them according to area. We'll finally plot those points using a list plot command. You can prove that the curve we're gonna get is continuous. So I can justify joining the points. There's nothing wrong with that mathematically. And then I'm gonna give a nice label here for patch area. It's running a little slower than what I thought. So let's see. We're up to a minute. Let's see if we can let it finish here. Okay, so we got patch area versus average density. All right, let's see. This says just a few more seconds. By the way, I love, love, love this parallel map shows kind of some built in time tracking of how much longer you have to wait. It's a game changer to me in terms of Mathematica. Love that feature. And then we get it. A beautiful result. Look at this. Oh man, this is so cool. So you've got area on the X axis, average density as a percentage of carrying capacity on the Y axis. You'll notice for areas below the critical threshold, nothing can survive. Once you pass that critical threshold, then we see a strictly positive continuous DAR relationship, just as was predicted by our theorem. Wonderful. Let me switch back over to the presentation. In fact, while I'm. No, let me just go ahead and switch, switch back to here. So we get back to our presentation. Let me make sure you all can see that screen. Yes, good. And we can pull up a nice density area relationship. Wonderful. What about non convex geometry? So I kind of have an arbitrary looking patch here where it's got three subcores, little peninsulas, two peninsulas and kind of a core area here. I'm calling them subcore 1, 2, and 3. Let's jump back over to Mathematica and let's see what would happen if we chose this kind of a non convex. So I'll switch my screen again. We're only gonna do this a couple of times, switching back and forth. Then I'll stick with the presentation. So let's go back to Mathematica. I'll share this again. I'm gonna pause for a second to make sure y' all can see that. Good. So I've contrived this region using region union and rectangles. I've done it in such a way that the area of this patch is one. Remember, I need a unit area to be able to stretch the patch bigger and smaller. We calculate E1D again. It's much bigger this time because the smaller peninsula, like things allow animals to enter into the hostile matrix at a higher rate than, say, a disk would. So that makes perfect sense that E1D would be higher. This command is the same as before. I've just changed the patch to now be omega 2. Let's just very quickly look at the visualization of this particular patch with our PDE numerically solved and put on top. So notice the minimum patch size here is 73. We're sitting at that level. Nothing's in the patch. As I increase the area, you'll notice almost immediately that bigger core is able to maintain a population, but the peninsulas cannot, even though the bigger core is subsidizing them. Because these peninsulas are so skinny, there's so many animals entering the matrix and dying or leaving, that these peninsulas are not able to maintain a noticeable population. But as I increase area, we will cross a threshold here where this left peninsula will all of a sudden be able to maintain a population. Notice the smaller peninsula still is somewhat dead. As I increase some more, we will cross another critical threshold for area and we will end up back with that peninsula being able to maintain a positive density. And as area goes to infinity, the patch will approach 100% of carrying capacity throughout. At least away from the boundary. You'll get a boundary layer that forms around the boundary, but in the interior you'll approach 100% of carrying capacity. Wonderful. Let me delete this. Manipulate. We've got a similar function to calculate single species average density. Only thing that's changed is now we have Omega 2. It's the exact same code as what I just ran above, except I'm using SSAD2, so I'm not going to go through the code again. I'm just going to run it. Hopefully it won't take quite as long to run this time. Let's see. We've got some kind of error here. What did we miss? All right, let me run this again. Let's make sure we have that and we have this. Okay, we'll delete this, run this one again, and try it again. That's interesting. Never seen that error message before, but we'll see if it will compute something or not. Well, that's live for you. You practice something, you never get an error message. When you try to do it in real time, you end up with things happening against Murphy's Law. If it can go wrong, it will. That's all right. Got some pictures in the PowerPoint if we have to get back over there. Oh, it produced it. I don't know what's going on with this message today, but that's okay. Oh, look at this. We get a result. This is beautiful. Notice it may be hard to see, but at those critical thresholds where the different little peninsulas came online or activated, there's a small dip in the DAR curve. And if you zoom into it close enough, you will see that the DAR slope is actually discontinuous there. So something very interesting is happening at that point. So let me jump back over to my PowerPoint. Switch back to that. I've got a little bit better picture on the PowerPoint and make sure y' all can see the screen good. So we got this. I did a little bit more detail here and actually gave you the density in each of those subcores. So green here represents the density just in subcore 1. Red represents the density in subcore 2. Blue represents the density in subcore 3. What we saw from Mathematica is once we passed that critical A1 threshold, we had the center or the bigger core part activated, but the others had no density, relatively no density in them. Once we got past the second critical threshold, the red core or sub S2 activated, causing a sharp increase in our overall DAR curve. Finally, as we cross that final threshold, the blue subcore activated, giving another sharp increase in the overall DAR curve, leading to DAR slope discontinuity. And so that's an interesting result there. Finally, one more thing. I'll show you. We wanted to know, could we really take this beautiful framework that Wolfram has been developing, this nice technology stack? Could we take a real world patch and do the same kind of calculation? The answer we found is yes. So we found a paper, Sozio and others in 2016, I literally found this picture in their paper and copied and pasted it into Mathematica. So let's jump over to Mathematica one more time and I want to show you that we can really apply this theory, this code, to a real world patch. So let's get back to the code. Hang with me one second. Wonderful. Everybody can see.

Good.

Here's the patch literally copied and pasted from their paper. We're going to use the command colornegate to negate the color and then image mesh to create a mesh out of that shape. We now have the shape of a real world patch. I believe this is a patch of forest located in Italy. Notice right now it has an area of 2400. To work in our framework, you've got to resize that so that the area is one you could use find root or something. I think I just played around with this scaling parameter here and got the area close to one. It doesn't have to be exact, so we're close there. You can calculate that E1D value. Again, I'm gonna skip a little bit of this code and I'll let y' all evaluate the manipulate on your own. You can play with the slider yourself. I wanna go down to the code. I've already pre run this code because it takes at least on my machine two or three minutes to do and we get a beautiful result here. What you'll notice is there is definitely a big time discontinuity in the DAR slope. And this is where one of those big peninsulas kind of activates. There's a couple more hidden here that's harder to see, but if you zoom in, you can find them. So yes, we can even use a real world patch from the literature. Put it in the Wolfram technology stack and output a DAR curve. This is so cool. All right, let me switch Back to the PowerPoint. We'll stay there. I won't switch back again unless there's questions later on. So let's go to the PowerPoint. I will pause for a second to make sure y' all are back with me. Yes, good. So again, DAR slope discontinuities arise. All right, let me talk a little bit about competing species. Dar by the way the code for the competing species is almost identical, you can use the same code. You just have to have your ND solve needs to work for the Lotka Volterra system, so you make a few minor modifications. The same code of work. I'll leave that as an exercise to the interested reader. All right. For competition, in the absence of a competitor, our model predicts that U and V have strictly positive continuous dar. What does this mean? If you have competition and you don't have the competitor there, you're back to the single species. So I'd expect it to be very much like it was before. For you. In blue, for any patch bigger than its critical area, you'll have a strictly positive continuous dar. For V in red, again, for area larger than this critical value, you'll have a strictly increasing positive dar. For the rest of the talk, I'm going to assume R2 is bigger than R1. That'll just make the results easier to display. Now, if you go to competition, let's take the weak case first. Whichever competitor has the larger minimum patch size, they always are guaranteed to have an overall positive dar. What do I mean overall? I can't guarantee that the curve is monotone. It could fluctuate. But overall, if I took a bunch of discrete points and I put a linear regression on it, I get a positive slope. So that's what I mean by overall. For the case of the competitor with the smaller minimum patch size, some interesting things happen. So let me set this up for you. First, in the literature, we know that the coexistent state UV will approach B1 Star and B2 Star, which is given just solely in terms of the competition strength as the area goes to infinity. We're going to make use of that in just a second. So here's the result. I can use Mathematica on a convex patch. I can see this beautiful outcome here. Area vs average density. For the competitor U who has the larger minimum hatch size. The green here is in the coexistent state. So this is U in the presence of V. From its initial kind of critical threshold up to the asymptotic value of B1 star, there's an overall increasing relationship. Again, it may wiggle. I can't guarantee it's monotone, but overall it's going up for V. If the average density at this value where U comes online, the coexistence first begins. If that average density is less than the asymptotic value for V in the coexistent state, then you will have a strictly Positive followed by overall positive. Dar, dar, give you one more case here and then I'll move along for sake of time. If you have a situation for V where the critical density, so the density of V where coexistence first starts, if that value is almost the same as B2 star, then you have strictly positive followed by overall neutral, which is very interesting. So that as the area increases, the model predicts that overall the density of V in the presence of its competitor doesn't change much. What's happening here is even though the area is getting bigger, it's benefiting not only V but also its competitor because of the competitive strength, somehow balancing out the benefit of the area increase and giving you a flat dar. I don't have time to show you the third case. The third case here. Let me jump to strong competition. In the case of strong competition we end up with something very interesting here. We have for V positive initially changing over to negative after that or a hump shaped dar. This has really interesting implications for conservation. Under this scenario, the V species would be better optimized density wise for an intermediate patch area, not one small, not one large, but somewhere in the middle. So that's a very interesting result that we got from this. I'll also mention that discontinuities can actually arise in DAR itself because of the competition. But let me not talk too much about that at this point. Now what about the non convex case? Let's come back to weak competition. For a disk we saw overall positive for B turning into neutral. If you go to our strange three peninsula looking thing, it qualitatively, not just quantitatively qualitatively changes the dar so that I go from overall neutral to one hump, two humps. So it changed even the qualitative shape of the DAR curve, just changing, keeping everything else the same, only changing patch geometry. What about the strong competition case? For a convex patch we saw V would have an overall hump shape. For the non convex patch we get not one, not two, but three humps and a bonus, a dead area in the middle. So you could maximize the density at three different island areas in the interior large patches you wouldn't be able to exist. In terms of V. Small patches you wouldn't be able to exist. Interestingly, we get three humps so that that convexity versus non convexity plays a crucial role. The final thing I want to do in my remaining 14 minutes and I'm going to be very respectful of your time here and in with a few minutes to answer questions, the final thing I Want to do is show you how this connects to some real world studies. So we got really excited as a research group when we took these theoretical results and found some examples in the literature where DAR was non standard and we tried to use our framework to explain what was going on there. Look, in my first example here, we had two ecologists, Buckley and Rough Garden in a 2006 paper studied the bronze anole and the grenad Islands. There's 13 islands here that these lizards were found in. These lizards are in competition with several insectivorous birds. Here's one in particular, so we'll call the anoles b the birds ewe in their system. At least from the paper we found that birds couldn't survive on smaller islands. This is what's reported in their paper. Birds are superior competitors. So we infer the semi strong competition based on this wonderful paper from these two ecologists. So how did they measure density? So they would take an island like Union island, which by the way looks beautiful here. As a mathematician, I've never had the benefit of doing math research in such a beautiful place. They were really lucky here. I'm sure it was tough walking through these woods, but man, what a beautiful island. So they would walk 100 meter transects as they went along. They counted the number of lizards that they saw. They turned that into a density. So that density on the y axis and the area of these 13 islands on the x axis. And so they were able to create estimates for island density. And so as we were studying this as a research group, we said, well, what's going on? This is a positive relationship with dar. Oh, here's a negative relationship with dar. That's hump shaped. That doesn't fit any of the standard theory. Oh, but we saw hump shape from our model where you could have this kind of a situation. Let's see if their study fits our modeling framework assumptions. Well, they have true islands. They're literally surrounded by ocean. A wave could wash a lizard out to sea if it gets too close to the shore, never to return. Same thing for birds. So absorbing boundary conditions makes sense. Birds can't survive on smaller islands, so in fact they have a larger minimum hatch size. In their paper, these two ecologists produced this particular graph where they showed for islands larger than area 1km squared, birds were present. For islands smaller than that crucial value, somewhere in 1km squared, no birds were present. So we chose an R value appropriately based on that, based on their research, birds are the superior competitors. So we infer the semi strong competition we assumed an average shape of a disk. Every one of these islands is a little bit different in shape. So what we're going to give is a prediction that's an average over the various insular geometries. So here's their hump shape observed dar. Here's what our model produced. Oh man, this is so exciting. Look at this. Around one you have a change from a positive DAR which matches what they found empirically, changing over to a negative DAR which again matches what they found empirically. But our model allows us to go a little bit deeper. We can talk a little bit about the mechanism that may be behind this hump shaped dar. So for lizards, red represents a lizard in an island by itself, no competitor. Purple represents a lizard in an island with competitors. Similarly for the birds, green represents present with lizards, blue would be birds by themselves. Typically as area decreases, we expect lizard density to decrease. Typically a single species, at least from our model should produce some kind of a positive relationship here. But what happens is because the patch area is decreasing, it's also affecting their competitor, but at a steeper rate. If you'll notice, the competitor dies out somewhere around 1. That releases anoles from the competition and actually creates a negative relationship. That's what's really behind this negative relationship is an ecological release. This is not an exact fit to their data. But what we can say is qualitatively our model does a really good job of explaining what's happening. So it's a possible mechanism, it's a plausible mechanism for this nonlinear and non standard dar. Let me give one more connection and then I'll pause from this and take some questions from y'. All. Anybody who has a question. In the second case we found a paper again by two ecologists, Sapilsky and Benkman, who in 2005 studied Red Cross bills. These birds live in patches of lodge pole pine trees. So this is an example of a lodgepole pine tree. These patches, by the way, these islands are huge. Just to give you a reference, like this is the state of Montana, this is Canada at the top, Idaho down here. I believe this large patch here would be Yellowstone National Park. So their islands are huge, they're forested islands surrounded by grasslands. This cross bill bird has evolved to only eat the seeds of that one pine tree. And in fact at different islands the birds beaks have evolved into different shapes to better fit the pine cones represented on those islands. The crossbill in these islands is in competition with a squirrel known as the pine squirrel, who also only eats the seeds of this one Pine tree. So if a squirrel or a bird gets very far away from one of their islands, there's nothing to eat. It becomes hostile. They will end up dying in their system. Just from the paper, from the research done in this wonderful paper, squirrels were absent from some of the islands present in others. And squirrels were the stronger competitors. But based on what we read, it seemed to be a weak competition between the species. How did they measure density? So they would take an island like the South Hills island patch, if you will. This is an actual picture of the South Hills forested island. Those are the pine trees surrounded by grassland. What they did is every 250 meters in this patch went through and gridded it out. They would stop for 10 minutes, count the number of crossbill calls. Using those calls, they could then get an estimate for crossbill density. That's what's represented on the Y axis. Island or forest patch areas here on the X axis. So what we see is for cross bill. In islands where the squirrel was absent, what kind of shape do we see? Oh, look at this. This is a positive dar. That's exactly what our model had suggested for a competition system. When the competitor is alone, we expect a positive at least overall dar and that's what we see in their result. Now, if you'll notice, the black dots here represent islands where the squirrel or the competitor is present. For those where the bird and the squirrel are together, what do we see? Well, the density almost doesn't change. It's somewhat neutral. Again, our model has predicted in certain situations you can see a neutral dart. So let's see, does their system, empirical system satisfy the modeling assumptions of our model. So they have forested islands where the matrix, the matrix rather is lethal to both. No food in the matrix. So that justifies absorbing boundary conditions. The crossbills can survive on smaller islands. Just reading the paper we pick up that squirrels have a larger minimum path size. So we pick an R accordingly. Competition is cited as the main driver of dynamics. Squirrels here are cited to be the stronger competitor. But overall the competition is weak. Inter specific competition is less than intra specific competition. Finally, again, we used an average island shape. So our predictions are average over all the different geometries are would be an average prediction. We again use Mathematica. Here's what they observed. Just to remind you, in islands where it was only bird, you had strictly positive DAR bird plus squirrel together neutral. Running the code on Mathematica, we end up with birds by themselves strictly increasing positive Darwin birds with squirrels together neutral. This is so cool. Qualitatively Our model is able to capture exactly what they observe in the real world. Now we weren't able to fit their data. This is not a quantitative type fit, but qualitatively our model fits theirs, qualitatively speaking. So what we can say is there's at least plausible mechanism here. And in fact we can go a little deeper and try to understand what's going on. If we look at the birds by themselves in red in competition in purple, squirrels in competition in green, by themselves in an island in blue, even though the area is decreasing when birds are present with their competitor, their density was somewhat neutral. Mechanistically, in our model you could see what's happening here is an ecological release. As the area is decreasing, not only is the bird density going down, but squirrel density, their competitors going down. And that somehow offsets the losses to less area and allows a near constant or neutral density area relationship. So once again what we see is ecological release here is a plausible mechanism for this nonlinear neutral dar. Let me end by summarizing. So the DAR form has been shown to be very inconsistent across taxa and even in within the same experiment to date, there's really been no theory to explain why DAR would be higher for intermediate areas. For example, in the lizard and bird scenario they really had a peak at an intermediate level of area. So there's to date had been no theory to explain that. So we wanted to see could a model explain what's going on, at least give a qualitative fit for our single species model? We saw that strictly positive DAR in logistically growing populations is always predicted in a competition setting. You could get overall positive positive than neutral or hump shape depending on the competition scenario, depending on the parameter ranges. Discontinuous DAR is present and possible anytime competition is involved. So you could end up with that scenario. When we compared our model to a couple of these empirical studies, we saw a good qualitative fit. What that means is ecological release is a plausible mechanism for these non standard DAR forms that don't fit anything in the literature today. Going forward. If I can leave you with two things. Number one, future empirical studies should be really careful about the scale of their study. So we saw totally different DAR structures depending on where you were in terms of patch area. So the scale of the study is going to play a role. We also saw, thanks to the beautiful work with Mathematica, the code that's written there by Wolfram Research, we're able to explore these really non convex geometries and see what kind of a role that plays and so that really should be taken into account for future study. Finally, if we're talking about species conservation, if you want to use dar, you need to be careful because you could get totally different conclusions based on is a competitor present or not? Do we have logistic growth, allele effect growth? There may be no single solution for any of this. So that kind of just be a closing remark. So with that, I will pause here and ask for questions.

Thank you, Professor Goddard. There are, if you want to take a peek at that Q and A pane in the new questions area. There are a couple of questions relevant to your presentation.

All right, so let me see. Okay, so I have from the first question Henrik is asking, please repeat, how DAR area divided by density can be negative. So density, that's a great question. Density, area is. DAR is a relationship. So when we say that the relationship is negative, we refer to referring to the slope. So if I were to take some islands with different increasing areas and I plot the corresponding densities, a negative relationship would mean as the area gets bigger, the density goes down. So that would represent the relationship there. Yeah. Great question. Second question from Kit. Could we generalize these results for more general boundary conditions, such as a robin or a robin or even nonlinear boundary conditions? That's a great question. In fact, you can. We've done some of this in my research group. You lose mathematically some of the ability to describe what's going on. There's some technical details involved on the mathematical side that we're still working on. But yes, that's a future direction for us to extend these results to more realistic boundary conditions. Yeah, wonderful question.

Well, you have done great on time. Thank you so much. I'm not seeing any new questions come into that Q and A pane. If you do have any, go ahead and type them in now. It's your opportunity. It was great to hear some your excitement on some of the newer functionality. So thank you for sharing that.

Absolutely. And if I can just make a comment if there is any questions while they're typing. I mean, really, from my perspective, I'm not a computational person. Say 10 years ago, Mathematica didn't have the functionality to do this kind of PDE solving over generic regions, maybe even five minutes ago, five years ago. I mean, you could do it back then, but you would have to really spend a lot of time with coding and get into libraries and all this. As a mathematician who's not a computational person, I really appreciate the technology stack y' all are building here. Really, really helps project like this go forward?

Well, I see that some of our developers are even joining us here online so they're hearing your feedback and that's, that's very constructive. Thank you. I see one more question coming in.

We have a question from Carlos. He says this model treats population as a continuous value. That's correct. What are my thoughts on using a discrete population with the same type of analysis be possible? I would say yes. I, I don't typically work with discrete population models so I don't want to give a definite answer. I would guess that you should be able to do something like this if for example, you had a discrete time model, continuous space or maybe discrete space, continuous time, there probably would be an analog to that. But again that's outside my research area so I don't want to speak too much there. But yes, we are considering the population to be continuous. Great question.